Probability energy flow calculation method for electricity-gas-heat area comprehensive energy system based on improved random response surface method

Abstract

The invention provides an electricity-gas-heat region comprehensive energy system probability energy flow calculation method based on an improved random response surface method, and belongs to the technical field of comprehensive energy system research, and the method comprises the following steps: building steady-state energy flow models of an electric power system, a thermodynamic system and a natural gas system; considering the influence of the correlation between photovoltaic output and wind power output in the same region on a system state variable, and performing standardization processing on an input variable by adopting Nataf transformation; and the electricity-gas-heat RIES probability energy flow is solved by adopting a Garleerkin projection method to optimize a random response surface method. According to the method, the probability energy flow of the regional integrated energy system is calculated, and the method has effectiveness, applicability and solution accuracy and rapidness.

Description

Technical Field

[0001] The invention relates to a probabilistic energy flow calculation method for an electric-gas-heat regional integrated energy system based on an improved random response surface method, and belongs to the research field of integrated energy systems. Background Art

[0002] Under the "dual carbon" goal, regional integrated energy systems containing renewable energy have the advantages of multi-energy complementarity, green and clean, and pollution-free, and have received widespread attention and attention. Since RIES covers energy systems in various energy forms such as power supply, heat supply, and gas supply, and there are uncertain factors such as intermittent output of renewable energy and load fluctuations, higher requirements are placed on RIES modeling and energy flow calculation.

[0003] For power systems, the methods for calculating probabilistic power flows are mainly divided into three categories: simulation method, approximate method, and analytical method. Among the simulation methods, Monte Carlo simulation (MCS) has a simple principle and is easy to implement. After calculating a large number of samples, the calculation results are accurate and used to verify the accuracy of power flow calculation. The representative point estimate method (PEM) among the approximate methods uses the probability distribution of the input random variables to approximate the digital characteristics of the output random variables to be calculated. The analytical method represented by the semi-invariant method obtains the semi-invariant of the output variable through calculation on the basis of the semi-invariants of each order of the input variables, thereby obtaining the probability distribution of the output variable. In addition to the above methods, the random response surface method has also gradually attracted attention in recent years. The random response surface method establishes a functional relationship model between the input random variables and the output variables, and uses fewer sample points to approximate the probability distribution characteristics of the system. Compared with the traditional MCS, SRSM greatly reduces the amount of calculation and can better maintain the calculation accuracy.

[0004] RIES has load fluctuations and uncertainties in renewable energy output, and their impact on the probabilistic energy flow of the system cannot be ignored. Therefore, when calculating the probabilistic energy flow of RIES, it is necessary to analyze the uncertainty factors. Chen Sheng et al. extended the concept of probabilistic power flow to the electric-gas hybrid IES and proposed a probabilistic power flow calculation method based on the Monte Carlo simulation method; Sun Juan et al. established a probabilistic energy flow model for the electric-thermal interconnected IES and proposed a Monte Carlo simulation method based on Latin hypercube sampling (LHS); Liao Xingxing et al. proposed a multi-point linear method based on the electric-thermal interconnected IES model to address the nonlinear problems existing in the system using the semi-invariant method; Zhang Rufeng et al. constructed a wind power and load uncertainty model based on the electric-gas IES steady-state energy flow model and proposed a calculation method for the semi-invariant of the output state variable; HRMassrur et al. proposed a probabilistic power flow calculation method based on the 2m+1 point estimation strategy based on the electric-gas-thermal IES model.

[0005] The randomly fluctuating wind and solar output in RIES is a system random variable, which not only has the characteristics of uncertainty, but also has correlation between photovoltaic and wind power generation in the actual system, and this correlation will affect the energy flow characteristics of the system. Li Chenyang et al. considered the impact of the correlation between the load of the power system and the natural gas system on the electric-gas IES; Tang Jianqing et al. considered the impact of the correlation between the load of the power system and the thermal system on the electric-thermal interconnected IES; Yang Jiahao considered the impact of the correlation between cold, hot, electric, gas loads and photovoltaics on the regional IES state variables; Qiu Yibin et al. analyzed the impact of the wind speed correlation of wind farms in different regions on the electric-gas-thermal IES; X. Zhu et al. proposed a probabilistic power flow analysis model for power systems that considers the correlation between various wind speeds, and introduced the association vector machine into the correlation modeling of wind speed through historical learning samples; Y. Xie et al. proposed an interval point estimation method based on the affine algorithm that considers the correlation of injected power based on the Nataf inverse transform and the p-box model, which takes into account the uncertain correlation coefficient. Considering that there is a certain degree of correlation between the wind and solar power outputs of the same system, if the correlation between wind power and photovoltaic power is ignored when studying the probabilistic energy flow of the electricity-gas-heat RIES, it will affect the stable operation of the electricity-gas-heat RIES.

[0006] In view of this, the invention firstly establishes a steady-state energy flow model of the power system, the thermal system and the natural gas system, and realizes the coupling of the power system, the thermal system and the natural gas system by using gas units, cogeneration units and power-to-gas equipment; secondly, the random variables such as photovoltaic output, wind power output and load are analyzed, and an electric-gas-heat RIES probabilistic energy flow model containing photovoltaic and wind power is established; thirdly, considering the influence of the correlation between photovoltaic output and wind power output in the same region on the system state variables, it is proposed to use the Garlerkin projection method to optimize the random response surface method (G-SRSM) to solve the electric-gas-heat RIES probabilistic energy flow; finally, in the example, by comparing G-SRSM with other methods to solve the electric-gas-heat RIES energy flow calculation model containing photovoltaic and wind power, the influence of the change of the correlation between photovoltaic output and wind power output on the system operation is analyzed, the effectiveness and applicability of the proposed model are verified, and at the same time, it is proved that G-SRSM has high accuracy and rapidity in solving the electric-gas-heat RIES energy flow calculation model containing photovoltaic and wind power. Summary of the invention

[0007] The present invention provides a probabilistic energy flow calculation method for an electric-gas-heat regional integrated energy system based on an improved random response surface method, which is used for probabilistic energy flow calculation of an integrated energy system, improves calculation efficiency and accuracy to ensure operation planning and optimized scheduling of RIES.

[0008] The technical solution of the present invention is: a probabilistic energy flow calculation method for an electric-gas-heat regional integrated energy system based on an improved random response surface method, and the specific steps of the method are as follows:

[0009] Step 1. Establish steady-state energy flow models for power system, thermal system and natural gas system;

[0010] Step 2: Considering the impact of the correlation between photovoltaic power output and wind power output in the same region on the system state variables, the input variables are standardized using Nataf transformation;

[0011] Step 3: Use the Garlerkin projection method to optimize the random response surface method to solve the electric-gas-thermal RIES probabilistic energy flow.

[0012] As a further solution of the present invention, the specific steps of Step 1 are:

[0013] Step 1.1, construct a steady-state energy flow model for the thermal system;

[0014] Step1.1.1, build hydraulic model;

[0015] The hydraulic process refers to the process in which hot water flows from the heat source side to the heat user side in the hot water pipe network. This process satisfies the flow balance equation and the loop pressure balance equation, as shown in equations (1)-(3).

[0016] A h m=m h (1)

[0017] Bh f =0(2)

[0018] h f =Km|m|(3)

[0019] Where: A h is the node-branch association matrix of the heating network; m is the pipeline flow of the heating network; m h is the flow column vector of the injection node; B is the heating network loop-branch association matrix; h f is the pressure head loss caused by the friction of water flowing through the pipe; K is the resistance coefficient of the heating network pipe.

[0020] Step1.1.2, build a thermal model;

[0021] The heat equation, pipe temperature drop equation and node mixed temperature equation of the thermal model are shown in equations (4)-(6).

[0022] Φ=C p m h (T s -T o )(4)

[0023]

[0024] ∑(m in T in )=(∑(m out )T out (6)

[0025] Where: Φ is the thermal power of each node in the thermal system; C p is the specific heat capacity of water; T s is the heating temperature; T o is the output temperature; T start and T end are the temperature at the beginning and end of the pipeline respectively; T a is the ambient temperature; λ is the heat transfer coefficient per unit length of the pipeline; L is the length of the heating network pipeline; m in and m out is the hot water flow rate of the injection node and the outflow node; T in and T out are the hot water temperatures at the injection node and outflow node, respectively.

[0026] Step 1.2, construct a steady-state energy flow model for the natural gas system;

[0027] When calculating the electricity-gas-heat RIES probabilistic energy flow, it is necessary to establish a natural gas system pipeline model and analyze the flow conditions of the natural gas system. The pipeline flow equation is shown in Equation (7)-Equation (8).

[0028]

[0029]

[0030] Where: G ij is the flow rate of pipeline ij under steady-state conditions; k ij is the pipeline parameter; s ij Indicates the flow direction of natural gas in the pipeline; p i and p j are the pressures at node i and node j respectively.

[0031] The natural gas network needs to satisfy the flow continuity equation, that is,

[0032] A g G=G g (9)

[0033] Where: A g is the node-branch association matrix of the natural gas network; G is the pipeline flow of the natural gas network; G g is the column vector of the flow rate injected into the natural gas network node.

[0034] Step 1.3, construct a steady-state energy flow model of the power system;

[0035] The steady-state energy flow model of the power system adopts the AC power flow model. The state variables of the system are mainly the node voltage amplitude and phase angle. The power flow equation expressed in polar coordinates is:

[0036]

[0037] Where: P i and Q i are the active power and reactive power injected into node i respectively; U i and U j are the voltage amplitudes of nodes i and j respectively; G ij and B ij The conductance and susceptance corresponding to node i and node j of the power system node admittance matrix are respectively; θ ij is the phase angle difference between node i and node j; n1 is the number of system nodes.

[0038] Step 1.4, build coupling element model;

[0039] Step1.4.1, build CHP unit model;

[0040] The CHP unit is an important coupling part of the power-gas-heat RIES. The unit consumes natural gas to generate electricity and supplies the waste heat in the production process to the heat network to achieve the goal of maximizing energy utilization. The back-pressure CHP unit meets (11).

[0041] C CHP =Φ CHP P CHP (11)

[0042] Where: C CHP is the thermal output of the CHP unit; Φ CHP The heat-to-electricity ratio of the unit; P CHP It is the electrical output of CHP unit.

[0043] Step 1.4.2, build the gas unit;

[0044] The gas unit is an important device for coupling the power system with the natural gas system. The unit is not only a load node of the natural gas system, but also a power node of the power system. The power generation of the gas unit supplied to the power system is proportional to the natural gas flow rate flowing into the unit from the natural gas system. The gas-to-electricity conversion relationship is shown in formula (12).

[0045] P gas =η G2P G gas (12)

[0046] Where: P gas is the power generation of the gas unit; η G2P The efficiency of gas-fired units in converting natural gas and electricity, G gas It is the gas consumption of gas unit.

[0047] Step1.4.3, build P2G device model;

[0048] When the natural gas load in RIES is at a low point or the output of renewable energy such as wind and solar power in the power system is high, the P2G device converts excess electricity into natural gas, achieving energy complementarity and improving energy utilization efficiency. The conversion process is shown in formula (13).

[0049]

[0050] Where: Q P2G is the amount of natural gas generated; P P2G The power consumed by the P2G device; H CH4 is the high calorific value of natural gas, which is 39MJ / m 3 ; η P2GConversion efficiency of P2G equipment.

[0051] Step 1.5, construct the electricity-gas-heat RIES energy flow model;

[0052] The RIES model is established by combining the power system, natural gas system, thermal system and coupling elements, realizing the multi-directional coupling between the subsystems in RIES. In addition, wind power and photovoltaic power generation have the characteristics of uncertainty and correlation, which increases the complexity of energy flow. In order to avoid the generation of local optimal solutions, the unified solution method is used to model and solve RIES as a whole, considering the mutual coupling between the subsystems, which can accurately describe and process the complex coupling relationship between the subsystems, avoid suboptimal solutions, and meet the accuracy requirements of energy flow calculation. Based on the unified solution method, a unified model for steady-state energy flow calculation of electric-gas-thermal RIES is invented and established, and the Newton-Raphson method is used to solve it. Its modified equations are shown in equations (14) and (15).

[0053]

[0054] Where: P and Q are the active and reactive power of the power system nodes respectively; ΔW is the input variable correction; ΔX is the state variable correction; J is the Jacobian matrix, which is composed of the power sub-matrix J ee 、Electrothermal arrayJ eh 、Electrical sub-array J eg 、Thermionic Array he 、Thermal arrayJ hh 、Hot gas array J hg 、Gas Electron Array J ge 、Qire array gh 、Gas Array J gg composition; V is voltage; θ is phase angle.

[0055] As a further solution of the present invention, the specific steps of Step 2 are:

[0056] In the electric-gas-heat RIES, as the penetration rate of wind power and photovoltaic power in the system continues to increase, the correlation between the output of the two is getting stronger and stronger, which has an important impact on the calculation of probabilistic energy flow. The output of renewable energy such as wind and light is closely related to the ambient temperature, wind speed, and weather, and has the characteristics of randomness, volatility, and seasonality. The wind and light output in the same area is correlated, and the degree of this correlation will change with the geographical location and regional energy consumption characteristics. At the same time, the wind-light correlation shows different significant change trends with seasonal changes. Ignoring this correlation will make it difficult to accurately describe the accurate scene of the output of new energy in RIES, affecting the calculation accuracy. Therefore, in RIES, the wind turbines and photovoltaic units connected to the power system at the same time have a connection in time and space distribution, and the output wind power and photovoltaic power will have a high correlation. The correlation is characterized by the correlation coefficient. Considering the influence of the correlation coefficient between the two, it is impossible to meet the requirement that the input variable is a standard input variable. Nataf transformation is often used to standardize it, mapping the input random variable to the standard normal distribution η.

[0057] Step 2.1. Input variables are transformed into non-independent standard normal distribution

[0058] For the input random variable matrix X, the correlation coefficient matrix of X can be calculated from the wind-solar correlation. Let the correlation coefficient matrix of the m-order random variable matrix X be C ρ , use the Spearman rank correlation coefficient calculation method to calculate the correlation coefficient ρ between random variables ij .

[0059]

[0060] Where: i , σ j X i , X j The standard deviation of cov(X i ,X j ) is X i , X j The covariance of .

[0061] The correlation coefficient matrix of the standard normal distribution vector η is C′ ρ , ρ′ ij is the correlation coefficient between random variables. Based on the principle of equal probability transformation, through the Nataf transformation theory, X that satisfies the non-normal distribution can be obtained through transformation.

[0062]

[0063] Where: μ i , μ j are random variables x i and xj The mean value of F i 1 、F j 1 are random variables x i and x j The inverse function of the cumulative distribution function; η i 、n j are random variables x i and x j The corresponding standard normal distribution; φ(·) is the standard normal distribution density function.

[0064] Step 2.2, obtain independent standard normal distribution;

[0065] Considering C′ ρ is a symmetric matrix, and C' is decomposed using the Cholesky decomposition method. ρ By decomposing, we can obtain the required lower triangular matrix L, and then we can obtain the required independent standard normal distribution ξ, as shown in formula (18).

[0066]

[0067] As a further solution of the present invention, the specific steps of Step 3 are:

[0068] Step3.1, selection of standard input variables;

[0069] Compared with the stochastic response surface method, the improved stochastic response surface method (SRSM) based on the Galerkin projection method has obvious advantages in solving the probabilistic energy flow of the integrated energy system. The following is the specific process of improving the stochastic response surface method:

[0070] Y=F(X)(19)

[0071] As shown in formula (19), in the SRSM-PLF calculation problem, the state variable to be calculated is the output response Y of model F, that is, the output random variable Y is the flow to be calculated, which is the node voltage, node pressure, pipeline flow and other parameters. The input random variable X is the uncertainty factor in the system, that is, the wind speed of the wind farm and the light intensity of the photovoltaic station studied by the invention. The calculation process is as follows.

[0072] Wind power output is mainly determined by wind speed. The uncertainty of wind speed is generally described by Weibull distribution, and the probability density function f(v) is expressed as formula (20).

[0073]

[0074] Where: k wis the shape coefficient; c is the scale coefficient, indicating the average wind speed; v is the wind speed.

[0075] According to the relationship between wind speed and wind turbine output power, the wind power output P(v) can be approximately expressed as formula (21).

[0076]

[0077] k2=-k1v ci (twenty three)

[0078] Where: P N is the rated power of the fan; v is the actual wind speed; v ci 、v N 、v co They are cut-in wind speed, rated wind speed and cut-out wind speed respectively; k1 and k2 are wind turbine power characteristic parameters.

[0079] Photovoltaic output is mainly determined by light intensity. According to historical data, light intensity satisfies Beta distribution in a certain time scale. Therefore, it can be inferred that photovoltaic output obeys Beta distribution in a certain time scale, and the probability density function is:

[0080]

[0081] Where: α and β are the shape parameters of the Beta distribution; P a P is the actual active power output of photovoltaic power; max It is the maximum active output of photovoltaic.

[0082] Step 3.2, express the output response Y with the Hermite chaotic polynomial about ξ and determine the number of unknown coefficients, as shown in equation (25) and equation (26).

[0083]

[0084] Where: a=[a0,a i ,a ij ,...] are the coefficients to be determined for the polynomial expansion; ξ is the n-dimensional independent standard normal distribution; n is the number of random variables input to the system; H m (·) is the m-th order Hermite orthogonal polynomial; H(ξ) is the Hermite orthogonal polynomial.

[0085] The number of coefficients to be determined in formula (25) is

[0086] N a =(m+n)! / m! n! (26)

[0087] Where N a is the number of coefficients to be determined in the expansion.

[0088] Step 3.3, determine the unknown coefficients of the polynomial expansion terms. The accuracy of SRSM is affected by random sampling points, and with the increase of sample points and the increase of dimension, the coefficient matrix is ​​prone to ill-conditioning, and the accuracy will deteriorate when the least squares method is used for solution. Therefore, the invention uses the Garlerkin projection method for calculation. First, a basic approximate model is established based on the least squares method, which can roughly describe the input-output relationship of the system; secondly, the Galerkin projection method is used to simultaneously project both sides of the above equation onto each orthogonal polynomial, that is, to project the complex nonlinear problem of the initial model onto a low-dimensional, easy-to-handle space, so as to more accurately solve non-convex problems and realize the linearization of the model, and equation (27) can be obtained.

[0089]

[0090] Where: <> is the inner product operation; H i , H j is the orthogonal polynomial of the i-th and j-th terms; y represents the output response of the system; a i are the unknown coefficients of the polynomial expansion. According to the definition of the inner product, we can use the orthogonality of orthogonal polynomials to get equation (28).

[0091]

[0092] Where: ij Indicates the orthogonality between orthogonal polynomials. When i=j, δ ij The value is 1, i≠j when δ ij The value is 0.

[0093] Rearranging the above formula yields:

[0094] a i =E[yH i (ξ)] / E[H i (ξ)H i (ξ)],(i=0,1,…,n)(29)

[0095] Where: E[·] represents the expected value operation.

[0096] Then the PCE coefficients can be obtained in turn through the above formula, where the denominator E[H i (ξ)H i (ξ)] is the expectation of the orthogonal polynomial, since the orthogonal polynomial is about each standard random variable ξ=[ξ1,…,ξ n ], the form of which is known, so the value of the denominator can be calculated quickly and easily. i (ξ)H i (ξ)] is about the multidimensional mixed orthogonal polynomial H i(ξ), so the calculation of the denominator can be regarded as the integral of each orthogonal polynomial E{[hi k(ξ k ) 2 ]}, that is

[0097]

[0098] Where: k represents the kth input random variable; hi k(ξ k ) is a one-dimensional orthogonal polynomial.

[0099] Among them, for each one-dimensional orthogonal polynomial integral E{[hi k(ξ k ) 2 ]} has a closed analytical solution, which is very easy to calculate, so the calculation is mainly from the numerator, E[yH i The calculation of (ξ)] is equivalent to the calculation of yH i (ξ). Due to the existence of the output response function y, it is almost impossible to obtain a closed solution for the molecule, and the output response function needs to be called a certain number of times. E[yH i (ξ)] can be solved by numerical methods or sampling methods. Considering that the sampling method has poor efficiency when the dimension is high, and the input variable is a two-dimensional random variable, the full factor numerical integration method is used to solve the numerator. The essence of the full factor integration method is the Gauss quadrature formula, which is characterized by using different weight functions for different random variables, and the nodes and weights of each dimension can be obtained through the Gaussian distribution table.

[0100] The Gauss quadrature formula is shown in formula (31).

[0101] E i (ξ)]=∫yH i (ξ)W(ξ)dξ (31)

[0102] Where: W(ξ) is the weight function of the standard normal distribution.

[0103]

[0104] Combining equation (27) with equation (31) and equation (32), the numerator of the unknown coefficient can be solved.

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

[0106] (1) Compared with MCS, SRSM, PEM, and UD-QMC, the G-SRSM algorithm proposed in the present invention can quickly and accurately obtain the RIES probability energy flow results.

[0107] (2) When gas units, CHP units and P2G equipment are connected to the power system, they participate in balancing the voltage fluctuations caused by the uncertainty of wind power output and photovoltaic output in the power system. The standard deviation of the voltage amplitude changes accordingly, affecting the operating status of the natural gas system and the thermal system, resulting in an increase in the fluctuation range of node gas pressure and pipeline flow.

[0108] (3) In order to avoid changes in the standard deviation caused by different wind and solar power fluctuation absorption capabilities due to differences in system structure, two scenarios are set up, and the correlation coefficients of photovoltaic and wind power are changed in the two scenarios respectively. The same rule is obtained, that is, as the positive and negative correlation between the two increases, the wind and solar power output increases or decreases, resulting in an increase in the standard deviation of the voltage amplitude of the power system node, the pipeline flow of the thermal system, and the gas pressure of the natural gas system node. BRIEF DESCRIPTION OF THE DRAWINGS

[0109] Figure 1 It is a flow chart of RIES probability energy flow calculation in the present invention;

[0110] Figure 2 Schematic diagram of the probability density of the voltage amplitude of the power system node 10 under different algorithms in the present invention;

[0111] Figure 3 Schematic diagram of gas pressure probability density of natural gas node 10 under different algorithms in the present invention;

[0112] Figure 4 Schematic diagram of the flow probability density of the thermal system pipeline 4 under different algorithms in the present invention;

[0113] Figure 5 It is a schematic diagram of the probability density of the voltage amplitude of the power system node 10 under different correlation coefficients in the present invention;

[0114] Figure 6 It is a schematic diagram of the cumulative probability of voltage amplitude of the power system node 10 under different correlation coefficients in the present invention;

[0115] Figure 7 It is a schematic diagram of the probability density of the flow rate of the thermal system pipeline 4 under different correlation coefficients in the present invention;

[0116] Figure 8 It is a schematic diagram of the cumulative probability of flow in the thermal system pipeline 4 under different correlation coefficients in the present invention;

[0117] Fig. 9 It is a schematic diagram of the probability density of gas pressure of natural gas node 10 under different correlation coefficients in the present invention;

[0118] Fig.10 It is a schematic diagram of the cumulative probability of gas pressure of the natural gas node 10 under different correlation coefficients in the present invention;

[0119] Fig.11It is a schematic diagram of the voltage amplitude distribution of power system nodes under different correlation coefficients in the present invention; DETAILED DESCRIPTION

[0120] Example 1: Figure 1-Figure 11 As shown in the figure, a probabilistic energy flow calculation method for the electric-gas-heat regional integrated energy system based on the improved random response surface method is proposed.

[0121] Step 1. Establish steady-state energy flow models for power system, thermal system and natural gas system;

[0122] Step 2: Considering the impact of the correlation between photovoltaic power output and wind power output in the same region on the system state variables, the input variables are standardized using Nataf transformation;

[0123] Step 3, use the Garlerkin projection method to optimize the random response surface method to solve the electric-gas-heat RIES probabilistic energy flow;

[0124] As a further solution of the present invention, the specific steps of Step 1 are:

[0125] Step 1.1, construct a steady-state energy flow model for the thermal system;

[0126] Step1.1.1, build hydraulic model;

[0127] The hydraulic process refers to the process in which hot water flows from the heat source side to the heat user side in the hot water pipe network. This process satisfies the flow balance equation and the loop pressure balance equation, as shown in equations (1)-(3).

[0128] A h m=m h (1)

[0129] Bh f =0(2)

[0130] h f =Km|m|(3)

[0131] Where: A h is the node-branch association matrix of the heating network; m is the pipeline flow of the heating network; m h is the flow column vector of the injection node; B is the heating network loop-branch association matrix; h f is the pressure head loss caused by the friction of water flowing through the pipe; K is the resistance coefficient of the heating network pipe.

[0132] Step1.1.2, build a thermal model;

[0133] The heat equation, pipe temperature drop equation and node mixed temperature equation of the thermal model are shown in equations (4)-(6).

[0134] Φ=C p m h (T s -T o )(4)

[0135]

[0136] ∑(m in T in )=(∑m out )T out (6)

[0137] Where: Φ is the thermal power of each node in the thermal system; C p is the specific heat capacity of water; T s is the heating temperature; T o is the output temperature; T start and T end are the temperature at the beginning and end of the pipeline respectively; T a is the ambient temperature; λ is the heat transfer coefficient per unit length of the pipeline; L is the length of the heating network pipeline; m in and m out is the hot water flow rate of the injection node and the outflow node; T in and T out are the hot water temperatures at the injection node and outflow node, respectively.

[0138] Step 1.2, construct a steady-state energy flow model for the natural gas system;

[0139] When calculating the electricity-gas-heat RIES probabilistic energy flow, it is necessary to establish a natural gas system pipeline model and analyze the flow conditions of the natural gas system. The pipeline flow equation is shown in Equation (7)-Equation (8).

[0140]

[0141] Where: G ij is the flow rate of pipeline ij under steady-state conditions; k ij is the pipeline parameter; s ij Indicates the flow direction of natural gas in the pipeline; p i and p j are the pressures at node i and node j respectively.

[0142] The natural gas network needs to satisfy the flow continuity equation, that is,

[0143] A g G=G g (9)

[0144] Where: A gis the node-branch association matrix of the natural gas network; G is the pipeline flow of the natural gas network; G g is the column vector of the flow rate injected into the natural gas network node.

[0145] Step 1.3, construct a steady-state energy flow model of the power system;

[0146] The steady-state energy flow model of the power system adopts the AC power flow model. The state variables of the system are mainly the node voltage amplitude and phase angle. The power flow equation expressed in polar coordinates is:

[0147]

[0148] Where: P i and Q i are the active power and reactive power injected into node i respectively; U i and U j are the voltage amplitudes of nodes i and j respectively; G ij and B ij The conductance and susceptance corresponding to node i and node j of the power system node admittance matrix are respectively; θ ij is the phase angle difference between node i and node j; n is the number of system nodes.

[0149] Step 1.4, build coupling element model;

[0150] Step1.4.1, build CHP unit model;

[0151] The CHP unit is an important coupling part of the power-gas-heat RIES. The unit consumes natural gas to generate electricity and supplies the waste heat in the production process to the heat network to achieve the purpose of maximizing energy utilization. The back-pressure CHP unit satisfies formula (11).

[0152] C CHP =Φ CHP P CHP (11)

[0153] Where: C CHP is the thermal output of the CHP unit; Φ CHP The heat-to-electricity ratio of the unit; P CHP It is the electrical output of CHP unit.

[0154] Step 1.4.2, build the gas unit;

[0155] The gas unit is an important device for coupling the power system with the natural gas system. The unit is not only a load node of the natural gas system, but also a power node of the power system. The power generation of the gas unit supplied to the power system is proportional to the natural gas flow rate flowing into the unit from the natural gas system. The gas-to-electricity conversion relationship is shown in formula (12).

[0156] P gas =η G2P G gas (12)

[0157] Where: P gas is the power generation of the gas unit; η G2P The efficiency of gas-fired units in converting natural gas and electricity, G gas It is the gas consumption of gas unit.

[0158] Step1.4.3, build P2G device model;

[0159] When the natural gas load in RIES is at a low point or the output of renewable energy such as wind and solar power in the power system is high, the P2G device converts excess electricity into natural gas, achieving energy complementarity and improving energy utilization efficiency. The conversion process is shown in formula (13).

[0160]

[0161] Where: Q P2G is the amount of natural gas generated; P P2G The power consumed by the P2G device; H CH4 is the high calorific value of natural gas, which is 39MJ / m 3 ; η P2G Conversion efficiency of P2G equipment.

[0162] Step 1.5, construct the electricity-gas-heat RIES energy flow model;

[0163] The RIES model is established by combining the power system, natural gas system, thermal system and coupling elements, realizing the multi-directional coupling between the subsystems in RIES. In addition, wind power and photovoltaic power generation have the characteristics of uncertainty and correlation, which increases the complexity of energy flow. In order to avoid the generation of local optimal solutions, a unified solution method is used to model and solve RIES as a whole, considering the mutual coupling between the subsystems, and can accurately describe and process the complex coupling relationship between the subsystems, avoid suboptimal solutions, and meet the accuracy requirements of energy flow calculation. Based on the unified solution method, the present invention establishes a unified model for steady-state energy flow calculation of electric-gas-thermal RIES, and uses the Newton-Raphson method to solve it. Its modified equations are shown in equations (14) and (15).

[0164]

[0165] Where: P and Q are the active and reactive power of the power system nodes respectively; ΔW is the input variable correction; ΔX is the state variable correction; J is the Jacobian matrix, which is composed of the power sub-matrix J ee 、Electrothermal arrayJ eh 、Electrical sub-array Jeg 、Thermionic Array he 、Thermal arrayJ hh 、Hot gas array J hg 、Gas Electron Array J ge 、Qire array gh 、Gas Array J gg composition; V is voltage; θ is phase angle.

[0166] As a further solution of the present invention, the specific steps of Step 2 are:

[0167] In the electric-gas-heat RIES, as the penetration rate of wind power and photovoltaic power in the system continues to increase, the correlation between the output of the two is getting stronger and stronger, which has an important impact on the calculation of probabilistic energy flow. The output of renewable energy such as wind and light is closely related to the ambient temperature, wind speed, and weather, and has the characteristics of randomness, volatility, and seasonality. The wind and light output in the same area is correlated, and the degree of this correlation will change with the geographical location and regional energy consumption characteristics. At the same time, the wind-light correlation shows different significant change trends with seasonal changes. Ignoring this correlation will make it difficult to accurately describe the accurate scene of the output of new energy in RIES, affecting the calculation accuracy. Therefore, in RIES, the wind turbines and photovoltaic units connected to the power system at the same time have a connection in time and space distribution, and the output wind power and photovoltaic power will have a high correlation. The correlation is characterized by the correlation coefficient. Considering the influence of the correlation coefficient between the two, it is impossible to meet the requirement that the input variable is a standard input variable. Nataf transformation is often used to standardize it, mapping the input random variable to the standard normal distribution η.

[0168] Step 2.1. Input variables are transformed into non-independent standard normal distribution

[0169] For the input random variable matrix X, the correlation coefficient matrix of X can be calculated from the wind-solar correlation. Let the correlation coefficient matrix of the m-order random variable matrix X be C ρ , use the Spearman rank correlation coefficient calculation method to calculate the correlation coefficient ρ between random variables ij .

[0170]

[0171] Where: i , σ j X i , X j The standard deviation of cov(X i ,X j ) is X i , X j The covariance of .

[0172] The correlation coefficient matrix of the standard normal distribution vector η is C′ ρ , ρ′ ij is the correlation coefficient between random variables. Based on the principle of equal probability transformation, through the Nataf transformation theory, X that satisfies the non-normal distribution can be obtained through transformation.

[0173]

[0174] Where: μ i , μ j are random variables x i and x j The mean of are random variables x i and x j The inverse function of the cumulative distribution function; η i 、n j are random variables x i and x j The corresponding standard normal distribution; φ(·) is the standard normal distribution density function.

[0175] Step 2.2, get independent standard normal distribution

[0176] Considering C′ ρ is a symmetric matrix, and C' is decomposed using the Cholesky decomposition method. ρ By decomposing, we can obtain the required lower triangular matrix L, and then we can obtain the required independent standard normal distribution ξ, as shown in formula (18).

[0177]

[0178] As a further solution of the present invention, the specific steps of Step 3 are:

[0179] Step 3.1. Selection of standard input variables

[0180] Compared with the stochastic response surface method, the improved stochastic response surface method (SRSM) based on the Galerkin projection method has obvious advantages in solving the probabilistic energy flow of the integrated energy system. The following is the specific process of improving the stochastic response surface method:

[0181] Y=F(X)(19)

[0182] As shown in formula (19), in the SRSM-PLF calculation problem, the state variable to be calculated is the output response Y of model F, that is, the output random variable Y is the flow to be calculated, which is the node voltage, node pressure, pipeline flow and other parameters. The input random variable X is the uncertainty factor in the system, that is, the wind speed of the wind farm and the light intensity of the photovoltaic station studied by the invention. The calculation process is as follows.

[0183] Wind power output is mainly determined by wind speed. The uncertainty of wind speed is generally described by Weibull distribution, and the probability density function f(v) is expressed as formula (20).

[0184]

[0185] Where: k w is the shape coefficient; c is the scale coefficient, indicating the average wind speed; v is the wind speed.

[0186] According to the relationship between wind speed and wind turbine output power, the wind power output P(v) can be approximately expressed as formula (21).

[0187]

[0188] Where: P N is the rated power of the fan; v is the actual wind speed; v ci 、v N 、v co They are cut-in wind speed, rated wind speed and cut-out wind speed respectively; k1 and k2 are wind turbine power characteristic parameters.

[0189] Photovoltaic output is mainly determined by light intensity. According to historical data, light intensity satisfies Beta distribution in a certain time scale. Therefore, it can be inferred that photovoltaic output obeys Beta distribution in a certain time scale, and the probability density function is:

[0190]

[0191] Where: α and β are the shape parameters of the Beta distribution; P a P is the actual active power output of photovoltaic power; max It is the maximum active output of photovoltaic.

[0192] Step 3.2, express the output response Y with the Hermite chaotic polynomial about ξ and determine the number of unknown coefficients, as shown in equation (25) and equation (26).

[0193]

[0194] Where: a=[a0,a i ,a ij,...] are the coefficients to be determined for the polynomial expansion; ξ is the n-dimensional independent standard normal distribution; n is the number of random variables input to the system; H m (·) is the m-th order Hermite orthogonal polynomial; H(ξ) is the Hermite orthogonal polynomial.

[0195] The number of coefficients to be determined in formula (25) is

[0196] N a =(m+n)! / m! n! (26)

[0197] Where N a is the number of coefficients to be determined in the expansion.

[0198] Step 3.3, determine the unknown coefficients of the polynomial expansion terms. The accuracy of SRSM is affected by random sampling points, and with the increase of sample points and the increase of dimension, the coefficient matrix is ​​prone to ill-conditioning, and the accuracy will deteriorate when the least squares method is used for solution. Therefore, the invention uses the Garlerkin projection method for calculation. First, a basic approximate model is established based on the least squares method, which can roughly describe the input-output relationship of the system; secondly, the Galerkin projection method is used to simultaneously project both sides of the above equation onto each orthogonal polynomial, that is, to project the complex nonlinear problem of the initial model onto a low-dimensional, easy-to-handle space, so as to more accurately solve non-convex problems and realize the linearization of the model, and equation (27) can be obtained.

[0199]

[0200] Where: <> is the inner product operation; H i , H j is the orthogonal polynomial of the i-th and j-th terms; y represents the output response of the system; a i are the unknown coefficients of the polynomial expansion. According to the definition of the inner product, we can use the orthogonality of orthogonal polynomials to get equation (28).

[0201]

[0202] Where: ij Indicates the orthogonality between orthogonal polynomials. When i=j, δ ij The value is 1, i≠j when δ ij The value is 0.

[0203] Rearranging the above formula yields:

[0204] a i =E[yH i (ξ)] / E[H i (ξ)H i (ξ)],(i=0,1,…,n)(29)

[0205] Where: E[·] represents the expected value operation.

[0206] Then the PCE coefficients can be obtained in turn through the above formula, where the denominator E[H i (ξ)H i (ξ)] is the expectation of the orthogonal polynomial, since the orthogonal polynomial is about each standard random variable ξ=[ξ1,…,ξ n ], the form of which is known, so the value of the denominator can be calculated quickly and easily. i (ξ)H i (ξ)] is about the multidimensional mixed orthogonal polynomial H i (ξ), so the calculation of the denominator can be regarded as the integral of each orthogonal polynomial E{[hi k(ξ k ) 2 ]}, that is

[0207]

[0208] Where: k represents the kth input random variable; hi k(ξ k ) is a one-dimensional orthogonal polynomial.

[0209] Among them, for each one-dimensional orthogonal polynomial integral E{[hi k(ξ k ) 2 ]} has a closed analytical solution, which is very easy to calculate, so the calculation is mainly from the numerator, E[yH i The calculation of (ξ)] is equivalent to the calculation of yH i (ξ). Due to the existence of the output response function y, it is almost impossible to obtain a closed solution for the molecule, and the output response function needs to be called a certain number of times. E[yH i (ξ)] can be solved by numerical methods or sampling methods. Considering that the sampling method has poor efficiency when the dimension is high, and the input variable is a two-dimensional random variable, the full factor numerical integration method is used to solve the numerator. The essence of the full factor integration method is the Gauss quadrature formula, which is characterized by using different weight functions for different random variables, and the nodes and weights of each dimension can be obtained through the Gaussian distribution table.

[0210] The Gauss quadrature formula is shown in formula (31).

[0211] E i (ξ)]=∫yH i (ξ)W(ξ)dξ (31)

[0212] Where: W(ξ) is the weight function of the standard normal distribution.

[0213]

[0214] Combining equation (27) with equation (31) and equation (32), the numerator of the unknown coefficient can be solved.

[0215] As a further solution of the present invention, the present invention further comprises step Step 4, and the specific steps of Step 4 are:

[0216] Step 4.1, a simulation analysis is performed on the RIES consisting of IEEE-30-node power system / IEEE-69-node power system, 20-node thermal system and 24-node natural gas system. At the same time, two simulation scenarios are set to verify the reliability of the simulation.

[0217] Scenario 1: RIES consisting of IEEE-30-node system, 20-node thermal system, and 24-node natural gas system.

[0218] Scenario 2: RIES consisting of an IEEE-69-node system, a 20-node thermal system, and a 24-node natural gas system.

[0219] Among them, wind turbines and photovoltaic units are connected to nodes 19 and 30 of the power system respectively. In the power system, nodes 14 and 15 are connected to nodes 12 and 13 of the natural gas system respectively through a gas unit. In the power system, nodes 7, 8, 11, and 13 are connected to nodes 3, 4, 5, and 6 of the natural gas system respectively through a CHP unit, among which nodes 7 and 8 of the natural gas system are connected to the thermal system S1, and nodes 9 and 10 are connected to S2. In scenario 1, a P2G device is connected to natural gas system 2 and power system node 14, and in scenario 2, a P2G device is connected to natural gas system 2 and power system node 2. It is assumed that the node loads of the power system, natural gas system, and thermal system (non-gas unit loads) meet the normal distribution with the expected initial value and the standard deviation of 5% of the expected value; the sampling scale of MCS is set to 10,000 times, and the convergence criterion of variance 5% is met.

[0220] Step 4.2: Verification of algorithm effectiveness

[0221] In order to verify the effectiveness of the G-SRSM method, the probabilistic energy flow calculation is carried out in the IEEE-69 node system in combination with G-SRSM, and the calculation results are compared with those of MCS, SRSM, Quasi-Monte Carlo method based on uniform design (UD-QMC), and point estimation method. Figure 2-Figure 4 , as shown in Table 1 and Table 2.

[0222] Table 1 Mean and standard deviation of system parameters under different algorithms

[0223]

[0224]

[0225] Table 2 Computational performance and accuracy under different algorithms

[0226]

[0227] Figure 2 , Figure 3 and Figure 4 The probability density curves of the voltage amplitude of the power system node 10, the gas pressure of the natural gas system node 10, and the flow rate of the thermal system pipeline 4 are shown in Tables 1 and 2. The expectation and standard deviation of the voltage amplitude of the power system node 10, the gas pressure of the natural gas system node 10, and the flow rate of the thermal system pipeline 4 are shown in Tables 1 and 2. The calculation results of the probability energy flow based on G-SRSM are very different from those of the Monte Carlo method. The maximum errors of the expectation and standard deviation do not exceed 3%, and the average error does not exceed 2%. The calculation time is much shorter than the comparison method in the article, indicating that the G-SRSM method proposed in the present invention has reliable calculation accuracy and efficiency.

[0228] Step 5.3: Impact of wind-solar correlation on system probabilistic energy flow

[0229] Considering the wind farm and photovoltaic power station in the same regional integrated energy system, their geographical location and climatic conditions are the same, so the wind speed and light intensity are generally strongly correlated, and the output power of the wind farm and the power of the photovoltaic power station are also strongly correlated. In order to verify the impact of the change in wind-solar correlation on the operating parameters of the regional integrated energy system, the proposed electricity-gas-heat RIES model and G-SRSM calculation method are used to set different correlation coefficient scenarios (the scenarios are -0.8, -0.4, 0.4, and 0.8, respectively) to calculate the probabilistic energy flow. Figure 5 , 7 , 9 are the probability density curves of the voltage amplitude of the power system node 10, the flow rate of the thermal system pipeline 4, and the gas pressure of the natural gas system node 10, respectively. Figure 6 , Figure 8 , Fig.10 The figures in the figure are the cumulative probability curves of the voltage amplitude of node 10 in the power system, the flow rate of pipeline 4 in the thermal system, and the gas pressure of node 10 in the natural gas system. Tables 3, 4, and 5 are the means and standard deviations of the voltage at node 10, the flow rate of pipeline 4 in the thermal system, and the gas pressure at node 10 in the natural gas system under different correlation coefficients, respectively.

[0230] Table 3 Mean and standard deviation of node 10 voltage under different correlation coefficients

[0231]

[0232] Table 4 Mean and standard deviation of flow rate of pipeline 4 under different correlation coefficients

[0233]

[0234] Table 5 Mean and standard deviation of air pressure at node 10 under different correlation coefficients

[0235]

[0236]

[0237] Combine Table 3 with Figure 5 It can be seen that in the process of the correlation coefficient increasing from -0.8 to 0.8, when photovoltaic and wind power show a positive correlation, the greater the correlation, the wind power and photovoltaic output will increase or decrease at the same time, and the voltage amplitude fluctuation affecting the power system node 10 will increase significantly; when photovoltaic and wind power show a negative correlation, the uncertainty of wind and solar output will increase significantly, the active power fluctuation in the system will be large, the standard deviation will increase accordingly, and the fluctuation range of the voltage amplitude will increase significantly. At the same time, with the change of the correlation between wind power and photovoltaic output, the peak value of PDF changes significantly, and the fluctuation range of the voltage amplitude at node 10 also changes significantly.

[0238] Combine Table 3 with Figure 6 It can be seen that as the correlation coefficient increases from -0.8 to 0.8, with the change of the correlation coefficient between wind power and photovoltaic output, there is a significant shift in the CDF curve between scenario 1 and scenario 2, indicating that the voltage amplitude of node 10 in the power system fluctuates greatly, which affects the energy flow of the power system.

[0239] Combined with Table 4, Table 5 and Figure 7-Figure 10 It can be seen that when the correlation coefficient increases from -0.8 to 0.8, when photovoltaic and wind power show correlation changes, the power fluctuation in the power system increases significantly. In order to avoid the risk of voltage over-limit, the natural gas system and thermal system adjust the energy flow respectively, and smooth the part of the fluctuating power through the coupling links such as power-to-gas, CHP units, and gas turbines in the system. The change in correlation causes the voltage amplitude to fluctuate violently, and the fluctuation of energy flow is transmitted to the thermal system and natural gas system through the coupling link, resulting in an increase in the peak of the PDF curve, a significant shift in the CDF curve, and a large fluctuation in the overall RIES energy flow.

[0240] Since the change of wind-solar correlation will affect the power flow of the entire power system, the operation status of the system can be reflected by analyzing the change of voltage amplitude at each node of the system. The distribution of voltage amplitude at power system nodes under different correlation coefficients is shown in the figure below: Fig.11 As shown. Fig.11It can be seen that the 69-node system has a stronger ability to absorb the fluctuation of wind and solar power output caused by uncertainty than the 30-node system. Therefore, under the same wind and solar correlation, the voltage fluctuation in scenario 1 is larger than that in scenario 2. Under the same scenario and different correlation conditions, (1) the correlation and uncertainty of wind power and photovoltaic output will cause node voltage fluctuations. When the node voltage amplitude is high or low, especially when the voltage amplitude exceeds the upper and lower limits, it will affect the stable operation of the power system; (2) Under the premise of meeting the system stability requirements, due to the enhanced correlation between wind power and photovoltaic output, the voltage amplitude of each node fluctuates greatly as a whole.

[0241] In summary, in the RIES probabilistic energy flow calculation, the change in the correlation between wind power and photovoltaic output will cause fluctuations in state parameters such as power system node voltage, natural gas system node pressure, and thermal system pipeline flow, and the uncertainty of system operation will increase.

[0242] The specific implementation modes of the present invention are described in detail above in conjunction with the accompanying drawings, but the present invention is not limited to the above implementation modes, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A probabilistic energy flow calculation method for an electric-gas-heat regional integrated energy system based on an improved stochastic response surface method, characterized in that: The specific steps of the method are as follows: Step 1. Establish steady-state energy flow models for power system, thermal system and natural gas system; Step 2: Considering the impact of the correlation between photovoltaic power output and wind power output in the same region on the system state variables, the input variables are standardized using Nataf transformation; Step 3: Use the Garlerkin projection method to optimize the random response surface method to solve the electric-gas-thermal RIES probabilistic energy flow.

2. According to claim 1, a probabilistic energy flow calculation method for an electric-gas-heat regional integrated energy system based on an improved stochastic response surface method is characterized by: The specific steps of Step 1 are: Step 1.1, construct a steady-state energy flow model for the thermal system; Step1.1.1, build hydraulic model; The hydraulic process refers to the process in which hot water flows from the heat source side to the heat user side in the hot water pipe network. This process satisfies the flow balance equation and the loop pressure balance equation, as shown in equations (1)-(3); A h m=m h (1) Bh f =0 (2) h f =Km|m| (3) Where: A h is the node-branch association matrix of the heating network; m is the pipeline flow of the heating network; m h is the flow column vector of the injection node; B is the heating network loop-branch association matrix; h f is the pressure head loss caused by the friction of water flowing through the pipe; K is the resistance coefficient of the heating network pipe; Step1.1.2, build a thermal model; The heat equation, pipe temperature drop equation and node mixed temperature equation of the thermal model are shown in equations (4)-(6); Φ=C p m h (T s -T o ) (4) ∑(m in T in )=(∑m out )T out (6) Where: Φ is the thermal power of each node in the thermal system; C p is the specific heat capacity of water; T s is the heating temperature; T o is the output temperature; T start and T end are the temperature at the beginning and end of the pipeline respectively; T a is the ambient temperature; λ is the heat transfer coefficient per unit length of the pipeline; L is the length of the heating network pipeline; m in and m out is the hot water flow rate of the injection node and the outflow node; T in and T out are the hot water temperatures at the injection node and outflow node, respectively; Step 1.2, construct a steady-state energy flow model for the natural gas system; When calculating the electricity-gas-heat RIES probabilistic energy flow, it is necessary to establish a natural gas system pipeline model and analyze the flow conditions of the natural gas system; the pipeline flow equation is shown in equations (7)-(8); Where: G ij is the flow rate of pipeline ij under steady-state conditions; k ij is the pipeline parameter; s ij Indicates the flow direction of natural gas in the pipeline; p i and p j are the pressures at nodes i and j, respectively; The natural gas network needs to satisfy the flow continuity equation, that is, A g G=G g (9) Where: A g is the node-branch association matrix of the natural gas network; G is the pipeline flow of the natural gas network; G g is the column vector of flow rate injected into the natural gas network node; Step 1.3, construct a steady-state energy flow model of the power system; The steady-state energy flow model of the power system adopts the AC power flow model. The state variables of the system are the node voltage amplitude and phase angle. The power flow equation expressed in polar coordinates is: Where: P i and Q i are the active power and reactive power injected into node i respectively; U i and U j are the voltage amplitudes of nodes i and j respectively; G ij and B ij The conductance and susceptance corresponding to node i and node j of the power system node admittance matrix are respectively; θ ij is the phase angle difference between node i and node j; n1 is the number of system nodes; Step 1.4, build a coupling element model; Step1.4.1, build CHP unit model; The back pressure CHP unit satisfies equation (11); C CHP =Φ CHP P CHP (11) Where: C CHP is the thermal output of the CHP unit; Φ CHP The heat-to-electricity ratio of the unit; P CHP It is the electrical output of CHP unit; Step 1.4.2, build the gas unit; The power generation of the power supply system of the gas unit is proportional to the natural gas flow rate flowing into the unit from the natural gas system. The gas-to-electricity conversion relationship is shown in formula (12); P gas =the G2P G gas (12) Where: P gas is the power generation of the gas unit; η G2P The efficiency of gas-fired units in converting natural gas and electricity, G gas is the gas consumption of the gas unit; Step1.4.3, build P2G device model; The P2G device converts excess electricity into natural gas when the natural gas load is at a low point in RIES or when the renewable energy output in the power system is high. The conversion process is shown in formula (13); Where: Q P2G is the amount of natural gas generated; P P2G The power consumed by the P2G device; H CH4 is the high calorific value of natural gas, which is 39MJ / m 3 ; η P2G The conversion efficiency of P2G equipment; Step 1.5, construct the electricity-gas-heat RIES energy flow model; The RIES model is established by combining the power system, natural gas system, thermal system and coupling elements to realize the multi-directional coupling between the subsystems in RIES. The RIES is modeled and solved as a whole using a unified solution method. Based on the unified solution method, a unified model for the steady-state energy flow calculation of the electric-gas-thermal RIES is established and solved using the Newton-Raphson method. The modified equations are shown in Equations (14) and (15). Where: P and Q are the active and reactive power of the power system nodes respectively; ΔW is the input variable correction; ΔX is the state variable correction; J is the Jacobian matrix, which is composed of the power sub-matrix J ee 、Electrothermal arrayJ eh 、Electrical sub-array J eg 、Thermionic Array he 、Thermal array J hh 、Hot gas array J hg 、Gas Electron Array J ge 、Qire array gh 、Gas Array J gg composition; V is voltage; θ is phase angle.

3. The method for calculating probabilistic energy flow of an electric-gas-heat regional integrated energy system based on an improved stochastic response surface method according to claim 1 is characterized by: In the Step 2, the input variables are standardized by using Nataf transformation to map the input random variables into a standard normal distribution; The specific steps of Step 2 are: Step2.1, input variables are transformed into non-independent standard normal distribution; For the input random variable matrix X, the correlation coefficient matrix of X is calculated based on the wind-solar correlation. Let the correlation coefficient matrix of the m-order input random variable matrix X be C ρ , use the Spearman rank correlation coefficient calculation method to calculate the correlation coefficient ρ between random variables ij ; Where: i , σ j X i , X j The standard deviation of cov(X i ,X j ) is X i , X j The covariance of The correlation coefficient matrix of the standard normal distribution vector η is C′ ρ , ρ′ ij is the correlation coefficient between random variables; based on the principle of equal probability transformation, through the Nataf transformation theory, X that satisfies the non-normal distribution can be transformed to obtain: Where: μ i , μ j are random variables x i and x j The mean value of F i 1 、F j 1 are random variables x i and x j The inverse function of the cumulative distribution function; η i 、n j are random variables x i and x j The corresponding standard normal distribution; φ(·) is the standard normal distribution density function; Step 2.2, obtain independent standard normal distribution; Considering C′ ρ is a symmetric matrix, and C' is decomposed using the Cholesky decomposition method. ρ Decompose and obtain the required lower triangular matrix L, that is, the required independent standard normal distribution ξ can be obtained, as shown in formula (18); 4. The method for calculating probabilistic energy flow of an electric-gas-heat regional integrated energy system based on an improved stochastic response surface method according to claim 1 is characterized by: The specific steps of Step 3 are: Step3.1, selection of standard input variables; The specific process of improving random response surface method is as follows: Y=F(X)(19) As shown in formula (19), in the SRSM-PLF calculation problem, the state variable to be calculated is the output response Y of model F, that is, the output random variable Y is the flow quantity to be calculated, including node voltage, node air pressure, and pipeline flow parameters; the input random variable X is the uncertainty factor in the system, that is, the wind speed of the wind farm and the light intensity of the photovoltaic station studied in the invention. The calculation process is described as follows. Wind power output is determined by wind speed. The uncertainty of wind speed is described by Weibull distribution, and the probability density function f(v) is expressed as formula (20); Where: k w is the shape coefficient; c is the scale coefficient, indicating the average wind speed; v is the actual wind speed; According to the relationship between wind speed and wind turbine output power, the wind power output P(v) is approximately expressed as formula (21); k2=-k1v ci (23) Where: P N is the rated power of the fan; v is the actual wind speed; v ci 、v N 、v co They are cut-in wind speed, rated wind speed, and cut-out wind speed; k1 and k2 are wind turbine power characteristic parameters; The photovoltaic output is determined by the light intensity. The photovoltaic output obeys the Beta distribution within a certain time scale, and the probability density function is: Where: α and β are the shape parameters of the Beta distribution; P a P is the actual active power output of photovoltaic power; max The maximum active power output of photovoltaic power; Step 3.2, the output response Y is expressed by the Hermite chaotic polynomial about ξ and the number of coefficients to be determined is determined, as shown in equation (25) and equation (26); Where: a=[a0,a i ,a ij ,...] are the coefficients to be determined in the polynomial expansion; ξ is the n-dimensional independent standard normal distribution; n is the number of random variables input to the system; H m (·) is the m-order Hermite orthogonal polynomial; H(ξ) is the Hermite orthogonal polynomial; The number of undetermined coefficients in formula (25) is: N a =(m+n)! / m! n! (26) where N a is the number of coefficients to be determined in the expansion; Step 3.3, determine the coefficients of the polynomial expansion terms; use the Garlerkin projection method to calculate, firstly, establish a basic approximate model based on the least squares method, which can describe the input-output relationship of the system; secondly, use the Galerkin projection method to simultaneously project both sides of the above equation onto each orthogonal polynomial, that is, project the complex nonlinear problem of the initial model onto a low-dimensional, easy-to-handle space, realize the linearization of the model, and obtain equation (27); Where: <> is the inner product operation; H i , H j is the orthogonal polynomial of the i-th and j-th terms; y represents the output response of the system; a i are the unknown coefficients of the polynomial expansion; according to the definition of the inner product, using the orthogonality of the orthogonal polynomial, we get equation (28); Where: ij Indicates the orthogonality between orthogonal polynomials. When i=j, δ ij The value is 1, i≠j when δ ij The value is 0; Rearranging the above formula gives: a i =E[yH i (ξ)] / E[H i (ξ)H i (ξ)],(i=0,1,…,n) (29) Where: E[·] represents the expected value operation; The calculation of the denominator is regarded as the integral of each orthogonal polynomial E{[hik(ξ k ) 2 ]}, that is Where: k represents the kth input random variable; hik(ξ k ) is a one-dimensional orthogonal polynomial; The Gauss quadrature formula is shown in formula (31); E[yH i (ξ)]=∫yH i (ξ)W(ξ)dξ (31) Where: W(ξ) is the weight function of the standard normal distribution; Combining equation (27) with equation (31) and equation (32), we can solve the numerator of the unknown coefficient.

Images