Probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory

Abstract

The probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory includes step 1: constructing a probabilistic model of wind turbine, photovoltaic and load output power; step 2: adding a set of estimation points on the basis of the point estimation method to form a new estimation point matrix and weight; step 3: performing deterministic calculation of harmonic power flow on the basis of the probabilistic model constructed in step 1 to solve the harmonic voltage of each node; step 4: establishing constraints on harmonic voltage, and constructing the optimal probabilistic distribution path of harmonic voltage through maximum entropy theory. The method of the present invention takes into account the uncertainty of wind turbines and photovoltaics and the characteristics of load volatility, obtains the harmonic distribution of the distribution network, and then better studies the level of random harmonics in the distribution network.

Description

Technical Field

[0001] The invention relates to the field of probabilistic harmonic power flow of distribution network, and in particular to a probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory. Background Art

[0002] In the context of low-carbon transformation of energy structure, building a new power system with new energy as the main body will become an important means to achieve low-carbon goals, which makes the power system transform from a deterministic system to a highly uncertain system. In these new power systems, the control characteristics of power electronic devices represented by synchronous machine rotor motion and electromagnetic transient equations have the problem of high power electronics, which will bring problems such as low inertia, weak anti-interference, harmonic pollution and broadband oscillation to the system, affecting the normal operation of the power grid.

[0003] Among them, harmonics are considered to be a public nuisance that threatens the safety of the power grid. Since harmonic currents are essentially generated by the conversion of fundamental currents in nonlinear elements, distributed power sources with nonlinear inverter devices connected to the power grid are receiving more and more attention due to their randomness and volatility as well as the harmonic problems they bring.

[0004] At present, the difficulty of considering probabilistic harmonic power flow calculation lies in how to take into account the following requirements:

[0005] 1) Calculate the statistical characteristics and probability density of the output random variables to further analyze the system harmonic distribution;

[0006] 2) It is necessary to ensure both calculation accuracy and improve calculation efficiency;

[0007] 3) Be able to handle nonlinear correlation problems between random variables.

[0008] Therefore, when calculating the harmonic flow of a power grid containing renewable energy, it is necessary to consider the objective factor of the uncertainty of renewable energy and reasonably select a probability model to maximize the power quality while ensuring the accuracy and efficiency of harmonic flow calculation.

[0009] Traditional point estimation is often combined with Gram-Charlier series expansion to obtain the probability density of input random variables. However, in practical problems, when the input variables have strong randomness and volatility, the output variables often have large high-order cumulants. When the third-order or fourth-order cumulants of the input random variables exceed a certain range, the Gram-Charlier series method will have a negative probability density function, resulting in large errors and complicated calculations. Summary of the invention

[0010] Based on the shortcomings of obtaining statistical characteristics of harmonic power flow calculation and low fitting accuracy of probability density function in the prior art, the present invention provides a probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory, which takes into account the uncertainty of wind power, photovoltaic and load output, so that the obtained probabilistic harmonic power flow calculation method is applicable to all working conditions; and while ensuring the calculation accuracy, the calculation efficiency is greatly improved.

[0011] The technical solution adopted by the present invention is:

[0012] The probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory includes the following steps:

[0013] Step 1: Construct a probability model of wind turbine, photovoltaic, and load output power;

[0014] Step 2: Add a set of estimation points based on the point estimation method to form a new estimation point matrix and weights;

[0015] Step 3: Based on the probability model constructed in step 1, perform deterministic calculation of harmonic power flow to solve the harmonic voltage of each node;

[0016] Step 4: Establish constraints on harmonic voltage and construct the optimal probability distribution path of harmonic voltage through maximum entropy theory.

[0017] In the step 1,

[0018] (1): The probability model of the wind turbine's output power is constructed as follows:

[0019]

[0020] In formula (1), f(v) is the probability density function of wind speed; v is the wind speed; is the average wind speed; σ v is the standard deviation of wind speed distribution; e is the base of the natural logarithm function.

[0021] The output power of the fan is shown in formula (2):

[0022]

[0023] In formula (2), P(v) is the output power of the fan; P N is the rated capacity of large-scale wind power; v is the wind speed;

[0024] v in 、v N 、v out They are respectively the cut-in wind speed, rated wind speed and cut-out wind speed (i.e. the output power at the rated wind speed) of the wind turbine.

[0025] When the probability density of the average wind speed in one hour and the wind turbine output power function are known, that is, when formula (1) and formula (2) are known, the average value of the probability density of the wind turbine output in that hour can be calculated using formula (3) and formula (4): With standard deviation

[0026] The average value of the probability density of wind turbine output It is expressed by formula (3):

[0027]

[0028] In formula (3), is the average value of the probability density of wind turbine output; E(P(v)) is the first-order moment of wind turbine output power, that is, the expected value of wind turbine output power; P(v) is the output power of the wind turbine; f(v) is the probability density function of wind speed; +∞ is an infinite positive number; ∞ is an infinite number; v is the wind speed;

[0029] v in 、v N 、v out They are respectively the cut-in wind speed, rated wind speed and cut-out wind speed (i.e. the output power at the rated wind speed) of the wind turbine.

[0030] Standard deviation of wind turbine output probability density It is expressed by formula (4):

[0031]

[0032] In formula (4), is the average value of the probability density of wind turbine output; is the standard deviation of the probability density of wind turbine output; P(v) is the output power of the wind turbine; f(v) is the probability density function of wind speed; It is the second-order origin moment of the difference between the wind turbine output power and the average value of the wind turbine output probability density; +∞ is an infinite positive number; v is the wind speed;

[0033] (2): The probability model of photovoltaic output power is constructed as shown in formula (5):

[0034]

[0035] In formula (5), f(r) is the probability density function of the photovoltaic output power; Γ(·) is the Gamma function, which is a type of function that extends the factorial function on real and complex numbers; r and r maxare the actual light intensity at a certain moment and the maximum light intensity in that period respectively; α and β are the shape parameter and scale parameter of Beta distribution, which are determined by the expected value and variance of the light intensity in each period.

[0036] Ignoring the nonlinear part of the photovoltaic cell output characteristics, the output power of the photovoltaic array is approximately considered to be proportional to the light intensity, as shown in formula (6).

[0037] P PV (r) = rA PV η PV (6);

[0038] In formula (6), P PV (r) is the photovoltaic output power; r is the light intensity; A PV is the total area of ​​the photovoltaic array; η PV is the photoelectric conversion efficiency of the photovoltaic array.

[0039] (3): The probability model of the load output power is constructed as shown in formula (7):

[0040]

[0041] In formula (7), f(P LD ) is the probability density function of the active power of the load; P LD are the active power of the load; μ P is the expected value of active power absorbed by the load; σ P is the standard deviation of the active power absorbed by the load; e is the base of the natural logarithm function.

[0042] The probability density function of load reactive power is shown in formula (8):

[0043]

[0044] In formula (8), f(Q LD ) is the probability density function of load reactive power; Q LD are the reactive power of the load; μ Q is the expected value of reactive power absorbed by the load; σ Q is the standard deviation of reactive power absorbed by the load; e is the base of the natural logarithm function.

[0045] The step 2 includes the following steps:

[0046] S2.1: Based on formulas (9) to (16), the point estimation theory is constructed as follows:

[0047] The estimated point x generated by connecting wind turbines, photovoltaics, and loads as input random variables to the system i,j, as shown in formula (9):

[0048] x i,j =μ i +ξ i,j σ i i=1,2,...,n,j=1,2,...,m (9);

[0049] In formula (9), x i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; n is the total number of connected wind turbines, photovoltaics and loads as input random variables; m is the total number of estimated points; i is the i-th random variable among the n connected wind turbines, photovoltaics and loads as input random variables; j is the j-th estimated point among the m estimated points; μ i and σ i are the expectation and standard deviation of each wind turbine, photovoltaic and load connected as input random variables.

[0050] Weight w i,j It can be expressed as formula (10):

[0051]

[0052] In formula (10), n is the total number of input random variables generated by uncertain wind turbines, photovoltaics and loads; m is the total number of estimation points; i is the i-th random variable among the n random variables generated by uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; w i,j The weight coefficient for the random variables generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken; It is the sum of the weight coefficients of the random variables generated by all uncertain wind turbines, photovoltaics and loads at all estimation points.

[0053] ξ i,j and w i,j The calculation of satisfies formula (11):

[0054]

[0055] In formula (11), w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; n is the number of wind turbines, photovoltaic and load connected as input random variables; ξ i,j For the estimated point x i,j The position coefficient of i,k is the normalized k-th order central moment, i.e., the k-th order central moment and standard deviation σ of the input random variables of wind turbines, photovoltaics and loads i The ratio of to the kth power.

[0056] Normalized k-th order central moment λ i,k As shown in formula (12):

[0057]

[0058] In formula (12), λ i,k is the normalized k-th order central moment; M k (X i ) is the input random variable X, which is connected to wind turbines, photovoltaics and loads i The k-th central moment of i The input random variables X are connected to wind turbines, photovoltaics and loads. i The expectation and standard deviation of ; k is the order of the central moment.

[0059] Connect wind turbines, photovoltaics and loads as input random variables X i The k-th order central moment M k (X i ), as shown in formula (13).

[0060]

[0061] In formula (13), M k (X i ) is the input random variable X, which is connected to wind turbines, photovoltaics and loads i The k-th central moment of i The input random variables X are connected to wind turbines, photovoltaics and loads. i Expectations; x i The wind turbine, photovoltaic and load are connected as one of the random variables in the input random variables; +∞ is an infinite positive number; -∞ is an infinite negative number; f(x i ) is the probability density function of the random variable; k is the order of the central moment.

[0062] The calculation process of point estimation is shown in formula (14):

[0063]

[0064] In formula (14), X nm is the mth estimated point of the nth random variable.

[0065] Assume that there are n random variables in the system, among which there are m probability sets. If each random variable can construct m estimation points, the total number of estimation points is m×n. Then the first 2m-1 moments of the random variable to be determined in the function can be obtained, thereby determining its probability distribution. One estimation point counts as n calculations, two estimation points require 2n calculations, and the accuracy is close to a third-order polynomial, and so on.

[0066] When each connected wind turbine, photovoltaic and load is used as the input random variable X i , take the jth estimated point x i,j When the remaining n-1 input random variables all take their expected values, the digital characteristics of the harmonic voltage function to be determined can be obtained by formula (15):

[0067]

[0068] In formula (15), X(i,j) is the input random variable of wind turbine, photovoltaic and load, and takes its jth estimated point x i,j When , the remaining n-1 random input variables are taken from the estimated point vectors at each mathematical expectation; μ n is the expected value of the input random variables of wind turbine, photovoltaic and load; i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; Z is the harmonic voltage of the random variable to be determined; w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; h(·) is the function value corresponding to the mathematical expectation of the wind turbine, photovoltaic and load as input random variables, which reflects the nonlinear relationship between the input variable and the output variable; l is the origin moment order of the variable to be determined, the harmonic voltage; E(Z l ) is the l-order origin moment of the random variable harmonic voltage to be determined.

[0069] According to the first two order origin moments, the mathematical expectation μ of the variable harmonic voltage Z can be obtained: Z and standard deviation σ Z As shown in formula (16):

[0070]

[0071] In formula (16), Z is the random variable harmonic voltage to be determined; E(Z 2 ) is the second-order origin moment of the random variable harmonic voltage to be determined; μ Z is the expected value of the random variable harmonic voltage to be determined; σ Z is the variance of the random variable harmonic voltage to be determined.

[0072] Using the Taylor series expansion of multivariate functions combined with formula (11), the output variable harmonic voltage Z is used as the input random variable X in wind turbines, photovoltaics and loads. i The expected value of is expanded to obtain the position coefficient ξ i,j And each estimated point x i,j The corresponding weight w i,j .

[0073] Each connected wind turbine, photovoltaic and load is used as the input random variable X i Position coefficient ξ i,j The calculation is shown in formula (17).

[0074]

[0075] In formula (17), ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; k is the order of the central moment; λ i,k is the kth order central moment normalized by the i-th random variable; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0076] Weight w i,j The calculation of is shown in formula (18).

[0077]

[0078] In formula (18), ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; k is the order of the central moment; λ i,k is the kth order central moment normalized by the i-th random variable; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0079] When the estimated point m = 3, according to formula (11), the position coefficient ξ can be obtained: i,j And each estimated point x i,j The corresponding weight w i,j , as shown in formula (19) and formula (20).

[0080]

[0081] In formula (19), λ i,4 is the normalized 4th-order central moment, also the input random variable X i The skewness coefficient can be obtained based on the first four origin moments of the random variable; ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; λ i,3 is the normalized third-order central moment; ξ i,3is the location coefficient of the third estimation point; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0082]

[0083] In formula (20), w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; w i,3 is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the third estimation point is taken; i,4 is the normalized 4th-order central moment, also the input random variable X i The skewness coefficient can be obtained based on the first four origin moments of the random variable; ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; λ i,3 is the normalized third-order central moment; ξ i,3 is the location coefficient of the third estimation point; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0084] When m=3, there are 3n estimation points in total. However, as shown in formula (11), each random variable X i There is an estimated point in the mathematical expectation μ i So there are n points corresponding to the same estimated vector. These n identical estimated point vectors only need to calculate the function value once, and the sum of their weight coefficients can be obtained, as shown in formula (21):

[0085]

[0086] In formula (21), w0 is the estimated point in the mathematical expectation μ i The weight at the point; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; w i,3 is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the third estimation point is taken; i,4 is the normalized fourth-order central moment; λ i,3 is the normalized third-order central moment.

[0087] Therefore, the statistical characteristics of the output variable corresponding to the three-point estimation method can be expressed as formula (22):

[0088]

[0089] In formula (22), X(i,j) is the input random variable of wind turbine, photovoltaic and load, and its jth estimated point x is taken. i,j , the remaining n-1 random input variables are all taken from the estimated point vectors at each mathematical expectation time; i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimated point among the m estimated points; μ n is the expected value of the input random variables of wind turbine, photovoltaic and load; i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; Z is the harmonic voltage of the random variable to be determined; w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; w0 is the weight coefficient when the random variables generated by the i-th uncertain wind turbine, photovoltaic and load are taken as the expected value; h(·) is the function value corresponding to the mathematical expectation of the connected wind turbine, photovoltaic and load as input random variables, which reflects the nonlinear relationship between the input variable and the output variable; l is the origin moment order of the variable to be determined, the harmonic voltage; E(Z l ) is the l-order origin moment of the random variable harmonic voltage to be determined.

[0090] S2.2: According to formulas (23) to (25), a set of estimation points is added to form a new estimation point matrix and weights, as follows:

[0091]

[0092] In formula (23), i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; ξ' i,j is the new position coefficient corresponding to each node; w' i,j is the new weight corresponding to each node; n is the number of wind turbines, photovoltaics and loads connected as input random variables.

[0093] Substitute ξ' in formula (23) i,1 , ξ' i,2 Substitute it back into formula (9) to get a new set of estimated points x' i,j , as shown in formula (24):

[0094]

[0095] In formula (24), x' i,j is the newly added estimation point; μ iis the expected value of the i-th connected wind turbine, photovoltaic and load as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among m estimation points; n is the number of connected wind turbines, photovoltaics and loads as input random variables.

[0096] Then, formula (22) estimates the origin moment of the desired output variable harmonic voltage, which can be updated to formula (25).

[0097]

[0098] In formula (25), X(i,j) is the random input variable generated by uncertain wind turbines, photovoltaics and loads when its jth estimated point x i,j , the remaining n-1 random input variables are all taken from the estimated point vectors at each mathematical expectation time; i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimated point among the m estimated points; μ n is the expectation of the random variables generated by uncertain wind turbines, photovoltaics and loads; i,j is the estimated point generated by the random input variable generated by each uncertain wind turbine, photovoltaic and load connected to the system; Z is the random variable harmonic voltage to be calculated; E'(Z l ) is the l-order origin moment of the newly added random variable harmonic voltage; w i,j is the weight coefficient of the random variable generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken; w0 is the weight coefficient of the random input variable generated by the uncertain wind turbine, photovoltaic and load when the estimation point is the expected value; h(·) is the function value corresponding to the mathematical expectation of the random input variable generated by the uncertain wind turbine, photovoltaic and load, which reflects the nonlinear relationship between the input variable and the output variable; X μ is the estimated point vector at the expected value; l is the origin moment order of the variable harmonic voltage to be determined; w i ' ,j The newly added weight coefficient is the random variable generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken.

[0099] In step three, the harmonic voltage of each node is solved by the harmonic network equation as shown in equation (26).

[0100] I h =Y h ·U h (26);

[0101] In formula (26), I h is the hth harmonic current; Y h is the hth harmonic admittance matrix; Uh is the hth harmonic voltage;

[0102] In step 4, a constraint model for harmonic voltage is established according to formulas (27) to (29), and the optimal probability distribution path of harmonic voltage is constructed by maximum entropy theory, as follows:

[0103] The maximum entropy of the harmonic voltage function is taken as the objective function, and the first four-order origin moments E of the harmonic voltage statistical characteristics are used as the Z,l (l=0,1,...,4) is the constraint condition to construct the model, as shown in formula (27).

[0104]

[0105] In formula (27), H is the entropy function; max is the largest set element, that is, the maximum value; ∞ is an infinite number; -∞ is an infinite negative number; f Z (z) is the probability density function of harmonic voltage; E Z,l is the origin moment of the harmonic voltage; l is the order of the origin moment.

[0106] The constraints of the maximum entropy model are shown in formula (28).

[0107]

[0108] In formula (28), H is the entropy function; st is the constraint condition; f Z (z) is the probability density function of harmonic voltage; ∞ is an infinite number; -∞ is an infinite negative number; E Z,l is the origin moment of the harmonic voltage; l is the order of the origin moment.

[0109] The Lagrange multiplier method is used for the entropy function H to construct a correction function, as shown in formula (29).

[0110]

[0111] In formula (29), L is the correction equation; H is the entropy function; f Z (z) is the probability density function of harmonic voltage; E Z,l is the origin moment of the harmonic voltage; +∞ is an infinite positive number; -∞ is an infinite negative number; l is the order of the origin moment; λ0, λ1, ..., λ4 are unknown constants.

[0112] By solving the nonlinear equation, the Lagrange coefficients λ0, λ1, ..., λ4 can be obtained, and the maximum entropy probability density function is shown in formula (30):

[0113]

[0114] In formula (30), f Z (z) is the probability density function of harmonic voltage; Z is the variable harmonic voltage to be determined; e is the base of the natural logarithm function; l is the order of the origin moment; λ0, λ1, ..., λ4 are the coefficients to be determined by the Lagrange multiplier method.

[0115] The present invention provides a probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory, and the technical effects are as follows:

[0116] 1) Considering that the traditional three-point estimation and series expansion have large errors in the value of input random variables during the calculation process, the present invention will improve them from two aspects: statistical feature extraction and probability density function fitting, and propose a probabilistic harmonic power flow algorithm based on improved point estimation combined with maximum entropy, so as to analyze and evaluate the system status.

[0117] 2) In general, the accuracy of point estimation is positively correlated with the number of points taken. However, in engineering practice, if a large number of estimation points are taken, the amount of system calculation will be greatly increased, resulting in reduced calculation efficiency. On the other hand, as can be seen from formula (11), the increase in the number of estimation points will introduce high-order standard central moments, which will cause a large number of high-order moments without actual physical meaning to be mixed in the subsequent complex calculations. In this regard, the present invention adds a set of estimation points x' on the basis of the three-point estimation. i,j Multiple sampling points are constructed to achieve accurate solution of the statistical characteristics of random variables. This not only ensures the calculation accuracy, but also avoids the introduction of high-order moments, making three-point sampling more practical.

[0118] 3) Based on the traditional three-point estimation method, the present invention changes the original weights by adding a group of estimation points, and reconstructs the uncertainty probability model including wind turbines, photovoltaics and loads, thereby improving the accuracy of statistical characteristics in probabilistic harmonic currents. The planning scheme can be applied to all working conditions.

[0119] 4) In response to the need to describe the statistical characteristics of the output results, the present invention introduces the maximum entropy theory to improve the statistical characteristics, which can accurately estimate all the statistical information of the output random variables, thereby obtaining a relatively complete harmonic power flow calculation theory based on point estimation. The present invention introduces the maximum entropy theory on the basis of the probability model constructed in (3) above, combines the improved three-point estimation with the maximum entropy theory, reconstructs the probability density function of the harmonic voltage, and accurately fits its probability distribution characteristics. BRIEF DESCRIPTION OF THE DRAWINGS

[0120] Figure 1 Flow chart of the method of the present invention.

[0121] Figure 2 The wiring diagram for the power distribution network of the system.

[0122] Figure 3 Schematic diagram of harmonic sources.

[0123] Figure 4 The following is the probability density function curve of the photovoltaic node under the working condition.

[0124] Figure 5 The probability density function curve of the fan node under the following working conditions.

[0125] Figure 6 It is the probability density function curve of the photovoltaic node under working condition 2.

[0126] Figure 7 It is the probability density function curve of the fan node under working condition 2.

[0127] Figure 8 It is the probability density function curve of photovoltaic nodes under three working conditions.

[0128] Fig. 9 It is the probability density function curve of the wind turbine node under three working conditions.

[0129] Fig.10 It is the probability density function curve of photovoltaic nodes under four working conditions.

[0130] Fig.11 It is the probability density function curve of the wind turbine node under four working conditions.

[0131] Fig.12 The following is a comparison chart of the voltage distortion rate distribution of each harmonic voltage. DETAILED DESCRIPTION

[0132] The probabilistic harmonic power flow calculation method based on the improved three-point estimation and maximum entropy theory changes the original weights by adding a group of estimation points on the basis of the traditional three-point estimation method, reconstructs the uncertainty probability model including wind turbines, photovoltaics and loads, and forms an estimation point matrix. The maximum entropy theory is introduced on the basis of the constructed probability model, and the improved three-point estimation is combined with the maximum entropy theory to reconstruct the probability density function of the harmonic voltage, so as to accurately fit its probability distribution characteristics. The method of the present invention takes into account the uncertainty of wind turbines and photovoltaics and the characteristics of load volatility, and obtains the harmonic distribution of the distribution network, thereby better studying the level of random harmonics in the distribution network. The obtained calculation method is applicable to all working conditions, and the fitting effect of the probability density function is closer to engineering practice, which verifies the feasibility and accuracy of the proposed method in dealing with probabilistic harmonic power flow problems. Specifically, it includes the following steps:

[0133] Step 1: Construct the probability model of wind turbine, photovoltaic and load according to formula (1) to formula (8):

[0134]

[0135] In formula (1), f(v) is the probability density function of wind speed; v is the wind speed; is the average wind speed; σ v is the standard deviation of wind speed distribution; e is the base of the natural logarithm function.

[0136]

[0137] In formula (2), P(v) is the output power of the fan; P N is the rated capacity of large-scale wind power; v is the wind speed; v in 、v N 、v out They are respectively the cut-in wind speed, rated wind speed and cut-out wind speed (i.e. the output power at the rated wind speed) of the wind turbine.

[0138] When the probability density of the average wind speed in one hour and the wind turbine output power function are known, that is, when formula (1) and formula (2) are known, the average value of the probability density of the wind turbine output in that hour can be calculated using formula (3) and formula (4): With standard deviation

[0139] The average value of the probability density of wind turbine output It is expressed by formula (3).

[0140]

[0141] In formula (3), is the average value of the probability density of wind turbine output; E(P(v)) is the first-order moment of wind turbine output power, that is, the expected value of wind turbine output power; P(v) is the output power of the wind turbine; f(v) is the probability density function of wind speed; +∞ is an infinite positive number; ∞ is an infinite number; v is the wind speed; v in 、v N 、v out They are respectively the cut-in wind speed, rated wind speed and cut-out wind speed (i.e. the output power at the rated wind speed) of the wind turbine.

[0142] Standard deviation of wind turbine output probability density It is expressed by formula (4).

[0143]

[0144] In formula (4), is the average value of the probability density of wind turbine output; is the standard deviation of the probability density of wind turbine output; P(v) is the output power of the wind turbine; f(v) is the probability density function of wind speed; is the second-order origin moment of the difference between the wind turbine output power and the average value of the wind turbine output probability density; +∞ is an infinite positive number; v is the wind speed; v in 、v N 、v out They are respectively the cut-in wind speed, rated wind speed and cut-out wind speed (i.e. the output power at the rated wind speed) of the wind turbine.

[0145] The probability density function of photovoltaic output power is shown in formula (5).

[0146]

[0147] In formula (5), f(r) is the probability density function of the photovoltaic output power; Γ() is the Gamma function, which is a type of function that extends the factorial function on real and complex numbers; r and r max are the actual light intensity at a certain moment and the maximum light intensity in that period respectively; α and β are the shape parameter and scale parameter of Beta distribution, which are determined by the expected value and variance of the light intensity in each period.

[0148] Ignoring the nonlinear part of the photovoltaic cell output characteristics, the output power of the photovoltaic array is approximately considered to be proportional to the light intensity, as shown in formula (6).

[0149] P PV (r) = rA PV η PV (6);

[0150] In formula (6), P PV (r) is the photovoltaic output power; r is the light intensity; A PV is the total area of ​​the photovoltaic array; η PV is the photoelectric conversion efficiency of the photovoltaic array.

[0151] The probability density function of the active power of the load is f(P LD ) is shown in formula (7).

[0152]

[0153] In formula (7), f(P LD ) is the probability density function of the active power of the load; P LD are the active power of the load; μ P is the expected value of active power absorbed by the load; σ Pis the standard deviation of the active power absorbed by the load; e is the base of the natural logarithm function.

[0154] The probability density function of load reactive power f(Q LD ) is shown in formula (8).

[0155]

[0156] In formula (8), f(Q LD ) is the probability density function of load reactive power; Q LD are the reactive power of the load; μ Q is the expected value of reactive power absorbed by the load; σ Q is the standard deviation of reactive power absorbed by the load; e is the base of the natural logarithm function. Step 2: Construct a point estimation theory with reference to formulas (9) to (16), and add a set of estimation points based on it according to formulas (23) to (25) to form a new estimation point matrix and weights.

[0157] The estimated point x generated by connecting wind turbines, photovoltaics and loads as input random variables to the system i,j As shown in formula (9).

[0158] x i,j =μ i +ξ i,j σ i i=1,2,...,n,j=1,2,...,m (9);

[0159] In formula (9), x i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; n is the total number of connected wind turbines, photovoltaics and loads as input random variables; m is the total number of estimated points; i is the i-th random variable among the n connected wind turbines, photovoltaics and loads as input random variables; j is the j-th estimated point among the m estimated points; μ i and σ i are the expectation and standard deviation of each wind turbine, photovoltaic and load connected as input random variables.

[0160] Weight w i,j It can be expressed as formula (10):

[0161]

[0162] In formula (10), n is the total number of input random variables generated by uncertain wind turbines, photovoltaics and loads; m is the total number of estimation points; i is the i-th random variable among the n random variables generated by uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; w i,jThe weight coefficient for the random variables generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken; It is the sum of the weight coefficients of the random variables generated by all uncertain wind turbines, photovoltaics and loads at all estimation points.

[0163] ξ i,j and w i,j The calculation of satisfies formula (11):

[0164]

[0165] In formula (11), w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; n is the number of wind turbines, photovoltaic and load connected as input random variables; ξ i,j For the estimated point x i,j The position coefficient of i,k is the normalized k-th order central moment, i.e., the k-th order central moment and standard deviation σ of the input random variables of wind turbines, photovoltaics and loads i The ratio of to the kth power.

[0166] Normalized k-th order central moment λ i,k As shown in formula (12).

[0167]

[0168] In formula (12), λ i,k is the normalized k-th order central moment; M k (X i ) is the input random variable X, which is connected to wind turbines, photovoltaics and loads i The k-th central moment of i The input random variables X are connected to wind turbines, photovoltaics and loads. i The expectation and standard deviation of ; k is the order of the central moment.

[0169] Connect wind turbines, photovoltaics and loads as input random variables X i The k-th order central moment M k (X i ) is shown in formula (13).

[0170]

[0171] In formula (13), M k (X i ) is the input random variable X, which is connected to wind turbines, photovoltaics and loads i The k-th central moment of i The input random variables X are connected to wind turbines, photovoltaics and loads. iExpectations; x i The wind turbine, photovoltaic and load are connected as one of the random variables in the input random variables; +∞ is an infinite positive number; -∞ is an infinite negative number; f(x i ) is the probability density function of the random variable; k is the order of the central moment.

[0172] The calculation process of the point estimate is shown in formula (14).

[0173] n random variables

[0174] In formula (14), X nm is the mth estimated point of the nth random variable.

[0175] Assume that there are n random variables in the system, among which there are m probability sets. If each random variable can construct m estimation points, the total number of estimation points is m×n. Then the first 2m-1 moments of the random variable to be determined in the function can be obtained, thereby determining its probability distribution. One estimation point counts as n calculations, two estimation points require 2n calculations, and the accuracy is close to a third-order polynomial, and so on.

[0176] When each connected wind turbine, photovoltaic and load is used as the input random variable X i Take the jth estimated point x i,j When the remaining n-1 input random variables all take their expected values, the digital characteristics of the harmonic voltage function to be determined can be obtained through formula (15).

[0177]

[0178] In formula (15), X(i,j) is the input random variable of wind turbine, photovoltaic and load, and takes its jth estimated point x i,j When , the remaining n-1 random input variables are taken from the estimated point vectors at each mathematical expectation; μ n is the expected value of the input random variables of wind turbine, photovoltaic and load; i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; Z is the harmonic voltage of the random variable to be determined; w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; h(·) is the function value corresponding to the mathematical expectation of the wind turbine, photovoltaic and load as input random variables, which reflects the nonlinear relationship between the input variable and the output variable; l is the origin moment order of the variable to be determined, the harmonic voltage; E(Z l ) is the l-order origin moment of the random variable harmonic voltage to be determined.

[0179] According to the first two order origin moments, the mathematical expectation μ of the variable harmonic voltage Z can be obtained: Z and standard deviation σ Z As shown in formula (16).

[0180]

[0181] In formula (16), Z is the random variable harmonic voltage to be determined; E(Z 2 ) is the second-order origin moment of the random variable harmonic voltage to be determined; μ Z is the expected value of the random variable harmonic voltage to be determined; σ Z is the variance of the random variable harmonic voltage to be determined.

[0182] Using the Taylor series expansion of multivariate functions combined with formula (11), the output variable harmonic voltage Z is used as the input random variable X in wind turbines, photovoltaics and loads. i The expected value of is expanded to obtain the position coefficient ξ i,j And each estimated point x i,j The corresponding weight w i,j .

[0183] Each connected wind turbine, photovoltaic and load is used as the input random variable X i Position coefficient ξ i,j As shown in formula (17):

[0184]

[0185] In formula (17), ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; k is the order of the central moment; λ i,k is the kth order central moment normalized by the i-th random variable; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0186] The weight is shown in formula (18):

[0187]

[0188] In formula (18), ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; k is the order of the central moment; λ i,kis the kth order central moment normalized by the i-th random variable; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0189] When the estimated point m = 3, according to formula (11), the position coefficient ξ can be obtained: i,j And each estimated point x i,j The corresponding weight w i,j , as shown in formula (19) and formula (20).

[0190]

[0191] In formula (19), λ i,4 is the normalized 4th-order central moment, also the input random variable X i The skewness coefficient can be obtained based on the first four origin moments of the random variable; ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; λ i,3 is the normalized third-order central moment; ξ i,3 is the location coefficient of the third estimation point; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0192]

[0193] In formula (20), w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; w i,3 is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the third estimation point is taken; i,4 is the normalized 4th-order central moment, also the input random variable X i The skewness coefficient can be obtained based on the first four origin moments of the random variable; ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; λ i,3 is the normalized third-order central moment; ξ i,3 is the location coefficient of the third estimation point; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points.

[0194] When m=3, there are 3n estimation points in total. However, as shown in formula (11), each random variable X i There is an estimated point in the mathematical expectation μ i So there are n points corresponding to the same estimated vector. These n identical estimated point vectors only need to calculate the function value once, and the sum of their weight coefficients can be obtained as formula (21):

[0195]

[0196] In formula (21), w0 is the estimated point in the mathematical expectation μ i The weight at the point; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; w i,3 is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the third estimation point is taken; i,4 is the normalized fourth-order central moment; λ i,3 is the normalized third-order central moment.

[0197] Therefore, the statistical characteristics of the output variable corresponding to the three-point estimation method can be expressed as formula (22):

[0198]

[0199] In formula (22), X(i,j) is the input random variable of wind turbine, photovoltaic and load, and its jth estimated point x is taken. i,j , the remaining n-1 random input variables are all taken from the estimated point vectors at each mathematical expectation time; i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimated point among the m estimated points; μ n is the expected value of the input random variables of wind turbine, photovoltaic and load; i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; Z is the harmonic voltage of the random variable to be determined; w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; w0 is the weight coefficient when the random variables generated by the i-th uncertain wind turbine, photovoltaic and load are taken as the expected value; h(·) is the function value corresponding to the mathematical expectation of the connected wind turbine, photovoltaic and load as input random variables, which reflects the nonlinear relationship between the input variable and the output variable; l is the origin moment order of the variable to be determined, the harmonic voltage; E(Z l ) is the l-order origin moment of the random variable harmonic voltage to be determined.

[0200] The construction of the newly added estimation points and weights is specifically shown in formula (23):

[0201]

[0202] In formula (23), i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; ξ' i,j is the new position coefficient corresponding to each node; w' i,j is the new weight corresponding to each node; n is the number of wind turbines, photovoltaics and loads connected as input random variables.

[0203] Substitute ξ' in formula (23) i,1 , ξ' i,2 Substitute it back into formula (9) to get a new set of estimated points x' i,j As shown in formula (24).

[0204]

[0205] In formula (24), x' i,j is the newly added estimation point; μ i is the expected value of the i-th connected wind turbine, photovoltaic and load as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among m estimation points; n is the number of connected wind turbines, photovoltaics and loads as input random variables.

[0206] Then, the origin moment of the output variable harmonic voltage estimated by formula (22) can be updated to formula (25).

[0207]

[0208] In formula (25), X(i,j) is the random input variable generated by uncertain wind turbines, photovoltaics and loads when its jth estimated point x i,j , the remaining n-1 random input variables are all taken from the estimated point vectors at each mathematical expectation time; i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimated point among the m estimated points; μ n is the expectation of the random variables generated by uncertain wind turbines, photovoltaics and loads; i,j is the estimated point generated by the random input variable generated by each uncertain wind turbine, photovoltaic and load connected to the system; Z is the random variable harmonic voltage to be calculated; E'(Z l ) is the l-order origin moment of the newly added random variable harmonic voltage; w i,jis the weight coefficient of the random variable generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken; w0 is the weight coefficient of the random input variable generated by the uncertain wind turbine, photovoltaic and load when the estimation point is the expected value; h(·) is the function value corresponding to the mathematical expectation of the random input variable generated by the uncertain wind turbine, photovoltaic and load, which reflects the nonlinear relationship between the input variable and the output variable; X μ is the estimated point vector at the expected value; l is the origin moment order of the variable harmonic voltage to be determined; w' i,j The newly added weight coefficient is the random variable generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken.

[0209] Step 3: Based on the probability model constructed in step 1, the harmonic power flow of the example is deterministically calculated, and the harmonic voltage of each node is solved by the harmonic network equation as shown in formula (26).

[0210] I h =Y h ·U h (26);

[0211] In formula (26), I h is the hth harmonic current; Y h is the hth harmonic admittance matrix; U h is the hth harmonic voltage;

[0212] Step 4: Establish a constraint model for harmonic voltage by referring to formulas (27) to (29), and construct the optimal probability distribution path of harmonic voltage through maximum entropy theory, as shown in formula (30).

[0213] The maximum entropy of the harmonic voltage function is taken as the objective function, and the first four-order origin moments E of the harmonic voltage statistical characteristics are used as the Z,l (l=0,1,...,4) is the constraint condition to construct the model, as shown in formula (27).

[0214]

[0215] In formula (27), H is the entropy function; max is the largest set element, that is, the maximum value; ∞ is an infinite number; -∞ is an infinite negative number; f Z (z) is the probability density function of harmonic voltage; E Z,l is the origin moment of the harmonic voltage; l is the order of the origin moment.

[0216] The constraints of the maximum entropy model are shown in formula (28).

[0217]

[0218] In formula (28), H is the entropy function; st is the constraint condition; f Z (z) is the probability density function of harmonic voltage; ∞ is an infinite number; -∞ is an infinite negative number; E Z,l is the origin moment of the harmonic voltage; l is the order of the origin moment.

[0219] The Lagrange multiplier method is used for the entropy function H to construct a correction function, as shown in formula (29).

[0220]

[0221] In formula (29), L is the correction equation; H is the entropy function; f Z (z) is the probability density function of harmonic voltage; E Z,l is the origin moment of the harmonic voltage; +∞ is an infinite positive number; -∞ is an infinite negative number; l is the order of the origin moment; λ0, λ1, ..., λ4 are unknown constants.

[0222] By solving the nonlinear equation, the Lagrange coefficients λ0, λ1, ..., λ4 can be obtained, and the maximum entropy probability density function is finally obtained as shown in formula (30).

[0223]

[0224] In formula (30), f Z (z) is the probability density function of harmonic voltage; Z is the variable harmonic voltage to be determined; e is the base of the natural logarithm function; l is the order of the origin moment; λ0, λ1, ..., λ4 are the coefficients to be determined by the Lagrange multiplier method.

[0225] Figure 1 A probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory is given. The specific calculation process is as follows:

[0226] 1) Calculate the probability and output size of wind power, photovoltaic power and load output power;

[0227] 2) Using the traditional three-point estimation theory, the estimation points and weights of the random variables generated by each uncertain wind turbine, photovoltaic and load are calculated, and the estimation point matrix is ​​constructed;

[0228] 3) Add a set of estimation points based on the three-point estimation method, obtain the corresponding weights, and construct a new estimation point matrix;

[0229] 4) Bring the reconstructed new estimated point matrix and weights into the example to perform deterministic power flow calculation and collect the power flow information of the system;

[0230] 5) Extract the fundamental voltage result in the fundamental power flow calculation and calculate the fundamental injection current;

[0231] 6) Assign values ​​to harmonic currents at wind turbine and photovoltaic nodes according to the harmonic content rate;

[0232] 7) Establish a harmonic power flow model, import the example initialization data, and calculate the harmonic admittance matrix through iteration;

[0233] 8) Obtain information such as harmonic power flow based on the obtained harmonic current and harmonic admittance matrix;

[0234] 9) Statistical probability information of harmonic voltage at distribution network nodes as output random variables, including expectation, variance, etc.;

[0235] 10) By reconstructing the probability density of harmonic voltage using the maximum entropy theory, we can more accurately estimate all the statistical information of the output random variables, thereby obtaining a relatively complete harmonic power flow calculation theory based on improved three-point estimation.

[0236] Verification example:

[0237] The parameters mentioned in the present invention are explained as follows:

[0238] (1): The example uses the IEEE33 node system for analysis. The system's power distribution network wiring diagram is shown in Figure 2 shown.

[0239] (2): When the wind turbine and photovoltaic are independent of each other, they are connected to nodes 16 and 30 respectively and used as harmonic source nodes. The schematic diagram of the harmonic source is as follows: Figure 3 As shown. The example takes the 5th and 7th harmonics generated by the harmonic source as an example. When analyzing the fundamental wave power flow, the wind turbine and photovoltaic nodes are simplified to PQ nodes. The harmonic content of the wind turbine and photovoltaic is taken from the literature (Chicco G, Schlabbach J, Spertino F. Characterization and assessment of the harmonic emission of grid-connected photovoltaic systems [C] / / 2005 IEEE Russia Power Tech. IEEE, 2005: 1-7. DOI: 10.1109 / IECON.2009.5415045) and:

[0240] Literature (Liu Bo, Yang Xu, Kong Fanlin, Ye Haizhong, Yu Hong. Control strategy of three-phase photovoltaic grid-connected inverter [J]. Transactions of China Electrotechnical Society, 2012, 27(08): 64-70. DOI: 10.19595 / j.cnki.1000-6753.tces.2012.08.009. ISSN: 1000-6753, CN: 11-2188 / TM), as shown in Table 1.

[0241] Table 1 Harmonic content rate

[0242]

[0243] (3): Based on the above, uncertain factors such as wind turbine, photovoltaic output and load are modeled. The expected value of the load node is the load prediction value, and its power factor is considered to remain constant during load and output fluctuations, with a coefficient of variation of 5%. The cut-in wind speed, rated wind speed and cut-off wind speed of the injection wind turbine are set to 3m / s, 14m / s and 25m / s respectively.

[0244] (4): Considering that wind speed fluctuation has a significant impact on the output of connected wind turbines and affects the harmonic distribution of the system, two scenarios with wind speeds of 7m / s and 13.5m / s were simulated for probabilistic harmonic power flow analysis.

[0245] (5): The distribution of connected wind turbines and photovoltaics is shown in Table 2. The active rated power of the wind turbine is P r And the maximum light intensity (1KW / m 2 ) when the maximum output P max All are set to 0.5MW.

[0246] Table 2 Distribution of wind turbines and photovoltaics

[0247]

[0248] (6) To prove the accuracy and feasibility of the proposed method, four operating conditions are simulated under the fifth harmonic and seventh harmonic conditions to analyze and calculate the probabilistic harmonic power flow, as shown below:

[0249] Working condition 1: 5th harmonic, wind speed 7 m / s;

[0250] Working condition 2: 5th harmonic, wind speed 13.5 m / s;

[0251] Working condition 3: 7th harmonic, wind speed 7 m / s;

[0252] Working condition 4: 7th harmonic, wind speed 13.5 m / s.

[0253] Based on the above-mentioned model, the statistical moments of Monte Carlo simulation (MCS), Point Estimate and Gram-Charlier (PEM&GC) and the improved three-point estimation method and maximum entropy theory (ITPEM&ME) proposed in this invention are compared. The statistical characteristics of harmonic voltage obtained by the system under the 5th harmonic using different algorithms are shown in Tables 3 and 4.

[0254] According to the analysis of the calculation results, when there are n wind turbines, photovoltaics and fluctuating loads connected to the system, the algorithm proposed by the present invention is used to generate a total of 2n+1 harmonic power flow calculations, which can quickly obtain statistical information such as the expectation and variance of the harmonic voltage, and the larger the number of samples, the more obvious the reduction in the amount of calculation. In obtaining the expected value of the 5th harmonic voltage, whether the wind speed is 7m / s or 13.5m / s, the calculation results of the three algorithms are basically the same, and are also relatively close to the calculation standard of the Monte Carlo simulation method. However, with the increase of the statistical moment order, the statistical characteristics of the three algorithms gradually differ.

[0255] Table 3 The first four order moments of harmonic voltage under different algorithms under the working conditions

[0256]

[0257] Table 4 The first four order moments of harmonic voltage of different algorithms under working condition 2

[0258]

[0259] In order to further verify the accuracy of the proposed algorithm, the present invention uses the probability harmonic voltage results obtained by 10,000 Monte Carlo simulations as a reference, takes the expected value and standard deviation error indicators as evaluation indicators of the algorithm accuracy, and analyzes the accuracy of the obtained probability harmonic power flow.

[0260] The relative error is defined as shown in formula (31):

[0261]

[0262] In formula (31), ε μ is the relative error of the expected value between algorithms; μ MCS represents the expected value of the output variable harmonic voltage obtained by MCS; μ PEM It is the expected value of the output variable harmonic voltage obtained by different point estimation methods.

[0263] The relative error is defined as shown in formula (32):

[0264]

[0265] In formula (32), ε σ is the relative error of the standard deviation between algorithms; σ MCS Represents the standard deviation of the output variable harmonic voltage obtained by MCS; σ PEM is the standard deviation of the output variable harmonic voltage obtained by different point estimation methods. According to formulas (31)~(32), the relative error results are shown in Table 5.

[0266] Table 5 Relative error of harmonic voltage statistical characteristics under different algorithms

[0267]

[0268] It can be clearly seen from Table 5 that the relative error of the ITPEM&ME method is smaller when solving the harmonic voltage of the two harmonic source nodes of wind and solar. The relative error of the two methods in the expected value is not much different, and both are less than 0.0001%, which fully meets the accuracy requirements of the distribution system operation analysis. In addition, in terms of the relative error of the standard deviation, the relative error obtained by the ITPEM&ME algorithm at 16 nodes is 0.0164%, 0.0164%, 0.0166% and 0.0163% lower than that of the PEM&GC method under four working conditions, and the relative error obtained at 30 nodes is 0.0366%, 0.0341%, 0.0374% and 0.0358% lower than that of the PEM&GC method under four working conditions, and the highest relative error reaches nearly 0.04%. This is because the standard deviation is calculated through the 1st and 2nd order origin moments, and the expected value is also the first order origin moment of the random variable. Since there are errors in the estimation of each order origin moment of the random variable, the errors of each order statistical moment gradually increase due to error accumulation, but they are all within the allowable range of the project.

[0269] Table 6 shows the average relative error of all nodes under four working conditions. It can be seen from Table 6 that in terms of the overall accuracy of the system, the expectations of harmonic voltages are basically consistent. For the standard deviation of harmonic voltages, the average relative errors of ITPEM&ME under four working conditions are 0.0233%, 0.0229%, 0.0236% and 0.0232% lower than those of PEM&GC, respectively. It can be seen that its average error is the smallest, which is significantly improved over the unoptimized PEM&GC. This also further illustrates that the ITPEM&ME method proposed in the present invention is more accurate, because when sampling the estimation points, ITPEM&ME further improves the tail characteristics of the sampling points, while PEM&GC ignores this point, resulting in the existence of errors. Therefore, when dealing with the problem of uncertain input variables, the ITPEM&ME method proposed in the present invention can obtain accurate system output responses using a small number of sampling points, which effectively improves the efficiency of the traditional model and has obvious advantages over other methods.

[0270] Table 6 Average relative error of harmonic voltage of different algorithms under four working conditions

[0271]

[0272] Under the four working conditions of the proposed parameter (6), the differences in probability density function reconstruction between the proposed method (ITPEM&ME) and PEM&GC and Monte Carlo are compared and analyzed. The comparison results are as follows: Figures 4 to 11 shown. Figures 4 to 11 They respectively represent the probability density function of the photovoltaic nodes under working condition one, the probability density function of the wind turbine nodes under working condition one, the probability density function of the photovoltaic nodes under working condition two, the probability density function of the wind turbine nodes under working condition two, the probability density function of the photovoltaic nodes under working condition three, the probability density function of the wind turbine nodes under working condition three, the probability density function of the photovoltaic nodes under working condition four, and the probability density function of the wind turbine nodes under working condition four.

[0273] from Figures 4 to 11It can be seen that the probability density estimation curves of the ITPEM&ME method proposed in the present invention are closest to the actual situation simulated by MCS under four different working conditions compared with the traditional PEM&GC method. The reason for this phenomenon is that when the system is connected to the wind turbine, due to the strong volatility of wind speed, the output of the wind turbine node generally has a high high-order cumulative amount, resulting in the maximum entropy theory only using the first four-order origin moments of the input harmonic voltage as constraints to obtain the probability density function through the Lagrange multiplier method when fitting the probability density, while PEM&GC obtains the first seven-order origin moments of the harmonic voltage and the semi-invariant pair series to expand and fit, so the accumulation causes errors. This proves that it is more reasonable to use the statistical moments of each order obtained from statistical data as constraints to approximate the probability density function of the random variable to be calculated, and further verifies that the proposed algorithm ITPEM&ME has a high degree of accuracy while also effectively ensuring the accuracy of the probability density estimation.

[0274] According to the obtained statistical data, the 5th and 7th harmonic voltages are calculated at each node, and the harmonic voltages under the 11th and 13th harmonics are calculated in the same way. The ratio of each harmonic voltage content to the fundamental voltage can be obtained to obtain the voltage distortion rate of each node. The distribution of the voltage distortion rate of each harmonic voltage is obtained through simulation, such as Fig.12 shown.

[0275] The impact of each harmonic of the comprehensive node on the system is Fig.12 It can be seen that the total voltage harmonic distortion rate generated by the ITPEM&ME algorithm is much lower than 4%, which meets the harmonic voltage standard of the public power grid.

[0276] The method of the present invention is run on the MATLAB platform of an Intel Core i5-8500 CPU computer, and the running calculation time under various working conditions is shown in Table 7.

[0277] Table 7 Comparison of calculation time under each algorithm

[0278]

[0279] It can be seen from Table 7 that the ITPEM&ME proposed in the present invention is slightly slower than PEM&GC in calculation speed, but much faster than the Monte Carlo method, and can quickly complete the harmonic power flow calculation in the distribution network. In summary, the method of the present invention has obvious high efficiency in terms of calculation efficiency and probability density fitting.

Claims

1. A probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory is characterized by The following steps are involved: Step 1: Construct a probability model of wind turbine, photovoltaic, and load output power; Step 2: Add a set of estimation points based on the point estimation method to form a new estimation point matrix and weights; Step 3: Based on the probability model constructed in step 1, perform deterministic calculation of harmonic power flow to solve the harmonic voltage of each node; Step 4: Establish constraints on harmonic voltage and construct the optimal probability distribution path of harmonic voltage through maximum entropy theory; In the step 3, the harmonic voltage of each node is solved by the harmonic network equation as shown in formula (26); I h =Y h ·U h (26); In formula (26), I h is the hth harmonic current; Y h is the hth harmonic admittance matrix; U h is the hth harmonic voltage; In step 4, a constraint model for harmonic voltage is established according to formulas (27) to (29), and the optimal probability distribution path of harmonic voltage is constructed by maximum entropy theory, as follows: The maximum entropy of the harmonic voltage function is taken as the objective function, and the first four-order origin moments E of the harmonic voltage statistical characteristics are used as the Z,l (l=0,1,...,4) is the constraint condition to build the model, as shown in formula (27); In formula (27), H is the entropy function; max is the largest set element, that is, the maximum value; ∞ is an infinite number; -∞ is an infinite negative number; f Z (z) is the probability density function of harmonic voltage; E Z,l is the origin moment of the harmonic voltage; l is the order of the origin moment; The constraints of the maximum entropy model are shown in formula (28); In formula (28), H is the entropy function; st is the constraint condition; f Z (z) is the probability density function of harmonic voltage; ∞ is an infinite number; -∞ is an infinite negative number; E Z,l is the origin moment of the harmonic voltage; l is the order of the origin moment; The Lagrange multiplier method is used for the entropy function H to construct a correction function, as shown in formula (29); In formula (29), L is the correction equation; H is the entropy function; f Z (z) is the probability density function of harmonic voltage; E Z,l is the origin moment of the harmonic voltage; +∞ is an infinite positive number; -∞ is an infinite negative number; l is the order of the origin moment; λ0,λ1,...,λ4 are constants to be determined; By solving the nonlinear equation, the Lagrange coefficients λ0, λ1, ..., λ4 are obtained, and the maximum entropy probability density function is finally obtained as shown in formula (30): In formula (30), f Z (z) is the probability density function of harmonic voltage; Z is the variable harmonic voltage to be determined; e is the base of the natural logarithm function; l is the order of the origin moment; λ0, λ1, ..., λ4 are the coefficients to be determined by the Lagrange multiplier method.

2. The probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory according to claim 1 is characterized by: In step 2, wind turbines, photovoltaics, and loads are connected as input random variables to generate the estimated point x generated by the system. i,j , as shown in formula (9): x i,j =μ i +ξ i,j s i i=1,2,...,n,j=1,2,...,m (9); In formula (9), x i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; n is the total number of connected wind turbines, photovoltaics and loads as input random variables; m is the total number of estimated points; i is the i-th random variable among the n connected wind turbines, photovoltaics and loads as input random variables; j is the j-th estimated point among the m estimated points; μ i and σ i are the expectation and standard deviation of each wind turbine, photovoltaic and load connected as input random variables; Weight w i,j It is expressed as formula (10): In formula (10), n is the total number of input random variables generated by uncertain wind turbines, photovoltaics and loads; m is the total number of estimation points; i is the i-th random variable among the n random variables generated by uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; w i,j The weight coefficient for the random variables generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken; is the sum of the weight coefficients of the random variables generated by all uncertain wind turbines, photovoltaics and loads at all estimation points; ξ i,j and w i,j The calculation of satisfies formula (11): In formula (11), w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; n is the number of wind turbines, photovoltaic and load connected as input random variables; ξ i,j For the estimated point x i,j The position coefficient of i,k is the normalized k-th order central moment, i.e., the k-th order central moment and standard deviation σ of the input random variables of wind turbines, photovoltaics and loads i The ratio of to the kth power; Normalized k-th order central moment λ i,k As shown in formula (12): In formula (12), λ i,k is the normalized k-th order central moment; M k (X i ) is the input random variable X, which is connected to wind turbines, photovoltaics and loads i The k-th central moment of i The input random variable X is the connection of wind turbines, photovoltaics and loads i The expectation and standard deviation of ; k is the order of the central moment; Connect wind turbines, photovoltaics and loads as input random variables X i The k-th central moment M k (X i ), as shown in formula (13); In formula (13), M k (X i ) is the input random variable X, which is connected to wind turbines, photovoltaics and loads i The k-th central moment of i The input random variable X is the connection of wind turbines, photovoltaics and loads i Expectations; x i The wind turbine, photovoltaic and load are connected as one of the random variables in the input random variables; +∞ is an infinite positive number; -∞ is an infinite negative number; f(x i ) is the probability density function of the random variable; k is the order of the central moment; The calculation process of point estimation is shown in formula (14): In formula (14), X nm is the mth estimated point of the nth random variable; Assume that there are n random variables in the system, among which there are m probability sets. If each random variable can construct m estimation points, the total number of estimation points is m×n. Then the first 2m-1 order moments of the random variable to be determined in the function can be obtained, thereby determining its probability distribution. When each connected wind turbine, photovoltaic and load is used as the input random variable X i , take the jth estimated point x i,j When , the remaining n-1 input random variables all take their expected values, then the digital characteristics of the harmonic voltage function to be determined are obtained by formula (15): In formula (15), X(i,j) is the input random variable of wind turbine, photovoltaic and load, and takes its jth estimated point x i,j When , the remaining n-1 random input variables are taken from the estimated point vectors at each mathematical expectation; μ n is the expected value of the input random variables of wind turbine, photovoltaic and load; i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; Z is the harmonic voltage of the random variable to be determined; w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; h(i) is the function value corresponding to the mathematical expectation of the wind turbine, photovoltaic and load as input random variables, which reflects the nonlinear relationship between the input variable and the output variable; l is the origin moment order of the variable to be determined, the harmonic voltage; E(Z l ) is the l-order origin moment of the random variable harmonic voltage to be determined; According to the first two order origin moments, the mathematical expectation μ of the variable harmonic voltage Z can be obtained: Z and standard deviation σ Z As shown in formula (16): In formula (16), Z is the random variable harmonic voltage to be determined; E(Z 2 ) is the second-order origin moment of the random variable harmonic voltage to be determined; μ Z is the expected value of the random variable harmonic voltage to be determined; σ Z is the variance of the random variable harmonic voltage to be determined.

3. The probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory according to claim 2 is characterized by: In the step 2, the Taylor series expansion of multivariate function is combined with formula (11) to convert the output variable harmonic voltage Z into the input random variable X when the wind turbine, photovoltaic and load are used. i The expected value of is expanded to obtain the position coefficient ξ i,j And each estimated point x i,j The corresponding weight w i,j ; Each connected wind turbine, photovoltaic and load is used as the input random variable X i Position coefficient ξ i,j The calculation is shown in formula (17); In formula (17), ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; k is the order of the central moment; λ i,k is the kth-order central moment of the normalized i-th random variable; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; Weight w i,j The calculation of is shown in formula (18); In formula (18), ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; k is the order of the central moment; λ i,k is the kth-order central moment of the normalized i-th random variable; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; When the estimated point m=3, according to formula (11), the position coefficient ξ is obtained i,j And each estimated point x i,j The corresponding weight w i,j , as shown in formula (19) and formula (20): In formula (19), λ i,4 is the normalized 4th-order central moment, also the input random variable X i The skewness coefficient can be obtained based on the first four origin moments of the random variable; ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; λ i,3 is the normalized third-order central moment; ξ i,3 is the location coefficient of the third estimation point; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; In formula (20), w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; w i,3 is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the third estimation point is taken; i,4 is the normalized 4th-order central moment, also the input random variable X i The skewness coefficient can be obtained based on the first four origin moments of the random variable; ξ i,j For each wind turbine, photovoltaic and load connected as the input random variable X i Position coefficient; λ i,3 is the normalized third-order central moment; ξ i,3 is the location coefficient of the third estimation point; n is the number of wind turbines, photovoltaics and loads connected as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; When m=3, there are 3n estimation points in total. According to formula (9), each random variable X i There is an estimated point in the mathematical expectation μ i So there are n points corresponding to the same estimated vector; these n identical estimated point vectors only need to calculate the function value once, and the sum of their weight coefficients can be obtained, as shown in formula (21): In formula (21), w0 is the estimated point in the mathematical expectation μ i The weight at the point; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; w i,3 is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the third estimation point is taken; i,4 is the normalized fourth-order central moment; λ i,3 is the normalized third-order central moment; Therefore, the statistical characteristics of the output variable corresponding to the three-point estimation method can be expressed as formula (22): In formula (22), X(i,j) is the input random variable of wind turbine, photovoltaic and load, and its jth estimated point x is taken. i,j When , the remaining n-1 random input variables are all taken from the estimated point vectors at each mathematical expectation time; i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimated point among the m estimated points; μ n is the expected value of the input random variables of wind turbine, photovoltaic and load; i,j is the estimated point generated by connecting wind turbines, photovoltaics and loads as input random variables to the system; Z is the harmonic voltage of the random variable to be determined; w i,j is the weight coefficient when the wind turbine, photovoltaic and load are connected as input random variables and the jth estimation point is taken; w0 is the weight coefficient when the random variable generated by the i-th uncertain wind turbine, photovoltaic and load takes the estimation point as the expected value; h(i) is the function value corresponding to the mathematical expectation of the connected wind turbine, photovoltaic and load as input random variables, which reflects the nonlinear relationship between the input variable and the output variable; l is the origin moment order of the variable to be determined, the harmonic voltage; E(Z l ) is the l-order origin moment of the random variable harmonic voltage to be determined.

4. The probabilistic harmonic power flow calculation method based on improved three-point estimation and maximum entropy theory according to claim 3 is characterized by: In the step 2, In formula (23), i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimation point among the m estimation points; ξ′ i,j is the new position coefficient corresponding to each node; w′ i,j is the new weight corresponding to each node; n is the number of wind turbines, photovoltaics and loads connected as input random variables; Substitute ξ′ in formula (23) i,1 ,ξ′ i,2 Substitute it back into formula (9) to obtain a new set of estimated points x′ i,j , as shown in formula (24): In formula (24), x′ i,j is the newly added estimation point; μ i is the expected value of the i-th connected wind turbine, photovoltaic and load as input random variables; i is the i-th random variable among the random variables generated by n uncertain wind turbines, photovoltaics and loads; j is the j-th estimation point among m estimation points; n is the number of connected wind turbines, photovoltaics and loads as input random variables; Then, formula (22) estimates the origin moment of the output variable harmonic voltage to be determined, which can be updated to formula (25); In formula (25), X(i,j) is the random input variable generated by uncertain wind turbines, photovoltaics and loads when its jth estimated point x i,j When , the remaining n-1 random input variables are all taken from the estimated point vectors at each mathematical expectation time; i is the i-th random variable among the n uncertain random variables generated by wind turbines, photovoltaics and loads; j is the j-th estimated point among the m estimated points; μ n is the expectation of the random variables generated by uncertain wind turbines, photovoltaics and loads; i,j is the estimated point generated by the random input variable generated by each uncertain wind turbine, photovoltaic and load connected to the system; Z is the harmonic voltage of the random variable to be determined; E'(Z l ) is the l-order origin moment of the newly added random variable harmonic voltage; w i,j is the weight coefficient of the random variable generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken; w0 is the weight coefficient of the random input variable generated by the uncertain wind turbine, photovoltaic and load when the estimation point is the expected value; h(i) is the function value corresponding to the mathematical expectation of the random input variable generated by the uncertain wind turbine, photovoltaic and load, which reflects the nonlinear relationship between the input variable and the output variable; X μ is the estimated point vector at the expected value; l is the origin moment order of the variable harmonic voltage to be determined; w i ' ,j The newly added weight coefficient is the random variable generated by the i-th uncertain wind turbine, photovoltaic and load when the j-th estimation point is taken.

Images