Comprehensive energy system toughness evaluation method based on Kolmogorov equation

Abstract

The invention relates to the technical field of comprehensive energy systems, in particular to a comprehensive energy system toughness evaluation method based on a Kolmogorov equation. The method comprises the following steps of: constructing a multi-energy system Iki stochastic differential equation, decomposing a system change process into a system slow change process and a system fast change process, deducing characteristics by utilizing historical data of extreme disasters, calculating expectation and variance of the extreme disasters, obtaining probability distribution of a system operation process, constructing a Kolmogorov forward equation, and calculating a Kolmogorov entropy value, so as to obtain a Kolmogorov entropy value; the method comprises the following steps of: calculating a drift coefficient and a diffusion coefficient, inputting data of an imminent predicted extreme disaster into a Kolmogorov forward equation, calculating probability distribution of emergency load shedding of the system by combining an optimal scheduling strategy of the multi-energy system, and outputting system toughness uncertainty evaluation according to the probability distribution. The method supports subsequent toughness improvement means, and is of great significance to guarantee safe and stable operation of social livelihood and comprehensive energy systems.

Description

Technical Field

[0001] The present invention relates to the technical field of integrated energy systems, and in particular to an integrated energy system resilience assessment method based on the Kolmogorov equation. Background Art

[0002] Frequent extreme disasters lead to huge losses due to power grid failures. The integrated energy system has attracted widespread attention because it can improve the system's ability to resist extreme disasters through the regulation and conversion of various resources. However, due to the involvement of multiple energy sources, the increase in the penetration rate of new energy, and the occurrence of multiple extreme disasters, it is difficult to evaluate the operating status of the integrated energy system in the face of extreme disasters, which greatly increases the difficulty of building subsequent prevention and control measures and system recovery strategies under failures.

[0003] In the past, resilience assessments of integrated energy systems were mostly based on a fixed extreme disaster as the background, and a certain fault scenario as the input of the calculation case to obtain a deterministic resilience assessment of the system. The shortcoming of this method is that the temporal and spatial distribution of extreme disasters is probabilistic, and the corresponding damage to the multi-energy system is also uncertain. Calculating system resilience based on a certain fault scenario alone and formulating response strategies based on this is not comprehensive and reference-based. Therefore, it is necessary to consider the impact of different extreme disasters on the integrated energy system and make a more reasonable multi-energy system resilience assessment to provide a reference for decision makers. Summary of the invention

[0004] The purpose of the present invention is to overcome the defects of the existing technology and propose a comprehensive energy system resilience assessment method based on the Kolmogorov equation. The Kolmogorov forward equation is used to learn the system state characteristics from historical data, calculate the drift coefficient and the diffusion coefficient, construct a stochastic differential equation, and calculate the probability distribution of the resilience index, thereby obtaining an uncertainty assessment of the system resilience.

[0005] To achieve the above object, the present invention adopts the following specific technical solutions:

[0006] The comprehensive energy system resilience assessment method based on the Kolmogorov equation provided by the present invention comprises the following steps:

[0007] S1. Collect extreme disaster parameters in the system area, use Monte Carlo simulation to simulate extreme disaster parameters, and obtain the probability distribution of the characteristics of the extreme disaster that will occur;

[0008] S2. Build an integrated energy system of electricity, gas, heat and storage, using cogeneration units, gas boilers, gas turbines and power-to-gas energy hubs to couple the energy systems;

[0009] S3, input the distribution and characteristics of extreme disasters into the integrated energy system as parameters, calculate the impact of extreme disasters on the integrated energy system, and obtain the probability distribution of system failure scenarios;

[0010] S4. Construct the Ito stochastic differential equation of the integrated energy system, decompose the change process of the integrated energy system into slow-changing process and fast-changing process to solve the probability distribution, construct the Kolmogorov forward equation to solve it, calculate the Kolmogorov entropy value, and obtain the drift coefficient and diffusion coefficient;

[0011] S5, inputting the predicted extreme disaster data into the Kolmogorov forward equation to solve the probability distribution of emergency load shedding of the integrated energy system;

[0012] S6. Construct system resilience indicators, draw resilience curves, and output uncertainty assessment of comprehensive energy system resilience.

[0013] Furthermore, in step S1, the system area extreme disaster parameters include wind speed, rainfall rate, ice accumulation, and atmospheric turbidity parameters, and the probability distribution of the characteristics of the extreme disaster that will occur is obtained through Monte Carlo simulation.

[0014] Further, in step S2, constructing the electricity-gas-heat-storage integrated energy system includes constructing a power grid system, a gas grid system, a heat grid system, and an energy transmission equation;

[0015] The operation constraints of the power grid system include AC power flow constraints, unit ramp constraints, and energy storage device constraints, as follows:

[0016] AC power flow constraints:

[0017]

[0018]

[0019] Where P i and Q i are the active and reactive power injection of node i, respectively. Injection is positive and output is negative. V i and V j is the voltage amplitude at nodes i and j; G ij and B ij They are the network impedance matrix Y bus The real and imaginary parts of represent the elements of the admittance matrix; θ ij is the difference in voltage phase angle between nodes i and j; n is the total number of nodes in the network; g ij and b ij The conductance and susceptance of the circuit ij, b sh,ij is the line shunt susceptance, P ij and Q ijare the active and reactive power flows of line ij respectively; V i min 、V i max are the lower and upper limits of the voltage amplitude; The upper limit of the active and reactive transmission capacity of the line;

[0020] Unit climbing constraints:

[0021] P i,t -P i,t-1 ≤ΔP i,up ;

[0022] P i,t-1 -P i,t ≤ΔP i,down ;

[0023] In the formula, ΔP i,up , ΔP i,down Indicates the upper and lower limits of the ramp rate of unit i per unit time;

[0024] Constraints on electrical energy storage devices:

[0025]

[0026] u n,c,t +u n,d,t ≤1;

[0027]

[0028] Where P ESS,n,t is the power input to the nth energy storage device; η n,c is the charging efficiency of the nth energy storage device; u n,c,t is the charging status of the energy storage device n, 1 for yes and 0 for no; u n,d,t is the discharge state of the energy storage device n, 1 for yes and 0 for no; are the charging and discharging power of energy storage device n; SOC n,t The amount of energy stored in the energy storage device; are the lower and upper limits of the charging power of the energy storage device n, respectively; are the lower and upper limits of the discharge power of the energy storage device n respectively; are the lower and upper limits of the storage capacity of the energy storage device n respectively.

[0029] Furthermore, in step S2, the operation of the gas network system includes the natural gas pipeline network equation and the gas compressor constraints, which are as follows:

[0030] Natural gas pipeline network equation:

[0031]

[0032] In the formula, G ij,t For natural gas pipelines ij Natural gas flow rate in m 3 / h, it is stipulated that: if the flow direction of natural gas in pipeline ij is: i→j, then G ij,t >0; is the gas output of the gas source; is the gas production of PtG; is the gas output of the gas storage device, where gas storage is negative and gas discharge is positive; is the demand for gas load; Pipeline ij Maximum and minimum values ​​of airflow; are the maximum and minimum values ​​of gas production of the gas source, respectively; is the gas consumption of GT; p i,t is the air pressure at the node; C ij For pipeline ij Correlation coefficient of physical properties; are the maximum and minimum values ​​of the air pressure at node i, respectively; is the set of nodes connected to node i; is the set of nodes connected to the gas source; is the set of nodes connected to the gas load; is the set of nodes connected to PtG; is a set of nodes connected to the gas storage device; is the set of nodes connected to GT; S G for gas pipeline assembly; A collection of natural gas pipelines with gas compressors. is the gas load loss;

[0033] Gas compressor constraints:

[0034]

[0035] In the formula, Λ n is the compression factor of the gas compressor n, and it is stipulated that: the gas flow in the pipeline containing the gas compressor is i→j; is the electrical power consumed by the gas compressor n, kW; 0-1 variable is the switch state of the gas compressor n, 1 is on, 0 is off; σ n is the power conversion factor of gas compressor n, kW / m 3 .

[0036] Furthermore, in step S2, the operation of the heat network system includes heat medium balance constraints, temperature mixing, pipeline heat loss and transmission delay, and heat exchange, as follows:

[0037] Heat medium balance constraints:

[0038]

[0039] In the formula, h k,t is the mass flow rate of the heat medium in the heat supply / return pipe k; are the sets of heat supply / return pipes starting and ending at node i; N H It is the node set of the heat supply and heat recovery network;

[0040] Temperature Mixing:

[0041] The heat supply / return pipe network contains three types of nodes: ordinary nodes, branch nodes and intersection nodes. The mixed relationship of the temperature of each node is as follows:

[0042] Ordinary nodes:

[0043]

[0044] Branch Node:

[0045]

[0046] Junction Node:

[0047]

[0048] In the formula, are the temperatures of the heat medium flowing into and out of the heat supply / return pipe k, ℃; S H It is a collection of supply / return heat pipes;

[0049] Temperature range of heat medium in supply / return pipe:

[0050]

[0051] In the formula, They are the upper and lower limits of the heat medium temperature in the supply / return heat pipes, respectively. The limit of the temperature in the return heat pipe is lower than that in the supply heat pipe.

[0052] Pipeline heat loss and transmission delay:

[0053] The heat power carried by the heat medium in the heat supply / return pipe is proportional to the mass flow rate and temperature of the pipe.

[0054]

[0055] In the formula, are the inflow and outflow heat powers in the supply / recovery pipe k, kW; are the inflow and outflow temperatures in the supply / return heat pipe, respectively, ℃; c is the specific heat capacity, kJ / (kg·℃);

[0056] The heat power loss and delay of the supply / return heat pipe is:

[0057]

[0058] In the formula, μ loss is the heat loss coefficient per unit length; l k is the length of the supply / return heat pipe, m; τ k,t is the transmission delay time of pipeline k at time t;

[0059] Heat exchange:

[0060] The heat power generated by the heat source is sent to the heat supply / recovery network through heat exchange.

[0061]

[0062] In the formula, is the thermal power flowing into the heating pipe k1; is the thermal power flowing out of the heat recovery pipe k2; η ex,1 The heat exchange efficiency between the heat source and the heat supply / return transmission pipeline; is the thermal power generated by the heat source m; They are respectively a collection of heat supply and heat return pipes connected to the heat source;

[0063] Heat users obtain heat from the supply / recovery network through heat exchange.

[0064]

[0065] Furthermore, in step S2, the energy transfer equation includes the electric-thermal characteristic equation of the cogeneration unit, the gas turbine model, and the power-to-gas model, which are as follows:

[0066] The electric-thermal characteristic equation of the cogeneration unit:

[0067]

[0068] In the formula, are the electric power and thermal power generated by the cogeneration unit connected to the grid node i and the heat network node j, kW; are the upper and lower limits of electric power and thermal power, kW respectively; are the slopes of the maximum and minimum steam inlet operating condition curves respectively; φ i,j are the slope and intercept of the minimum condensing condition curve, kW;

[0069] Gas Turbine Model:

[0070] Gas turbines generate electricity by consuming natural gas, providing power, gas loads,

[0071]

[0072] In the formula, The gas turbine connected to the grid node k and the natural gas node n generates electrical power, kW; for the efficiency of the gas turbine; is the gas consumed by the gas turbine, m 3 / h; are the upper and lower limits of the generated electric power respectively; is the conversion factor, which is used to unify the dimensions; Π is the higher calorific value of natural gas, kBtu / m 3 ; Φ is the energy conversion coefficient, kBtu / (kW·h);

[0073] Power-to-gas model:

[0074] Power-to-gas generates natural gas by consuming electricity, which is the power load, gas source,

[0075]

[0076] In the formula, is the gas consumed by the gas turbine, m 3 / h; is the electrical power generated by the gas turbine, kW; η n ptg is the efficiency of the gas turbine; G n ptg,max , G n ptg,min They are the upper and lower limits of gas generation respectively.

[0077] Furthermore, in step S3, the impact of extreme disasters on the integrated energy system is calculated, including the impact of extreme disasters on the source, the impact of extreme disasters on the grid, and the impact of extreme disasters on the load, as follows:

[0078] Impact of extreme disasters on sources:

[0079] The impact of ice disasters on wind turbines is as follows:

[0080]

[0081] P ice,t =P(1-g1M turbine );

[0082] Where M turbine is the ice accumulation of the fan, k turbine is the ice accumulation coefficient, S is the precipitation rate, V is the wind speed, P is the wind turbine output before correction, P ice is the corrected fan output, g1 is the ice accumulation factor;

[0083] The impact of sandstorms on the power generation of photovoltaic units is as follows:

[0084]

[0085] Where P pv Output for photovoltaic power; Y pv is the rated capacity of the photovoltaic array under standard test conditions; f pv is the photovoltaic derating factor, which is related to the surface conditions of the photovoltaic array; G T is the solar radiation intensity currently received by the photovoltaic array, kW / m 2 ; G T,STC is the incident light intensity under standard test conditions, constant, 1kW / m 2 ; T c is the current surface temperature of the photovoltaic array, °C; T c,STC is the temperature of the photovoltaic array under standard test conditions, °C; α P is the temperature correction factor;

[0086] Impact of extreme disasters on the Internet:

[0087] The impact of the ice disaster on the line is as follows:

[0088]

[0089] In the formula, a w1 and b w1 are two threshold values ​​of ice wind load, L WI is the load under the combined action of wind and ice, P f is the line disconnection failure rate;

[0090] Impact of extreme disasters on the Netherlands:

[0091] The impacts of the ice disaster on the Netherlands are as follows:

[0092]

[0093] In the formula, is the heat network transmission loss, is the heat loss of the building; T amb is the ambient temperature; K ij ,cn ij is the heat network transmission pipeline parameter, R i is the thermal resistance of the building, is the air transmission parameter.

[0094] Furthermore, in step S4, the Ito stochastic differential equation is as follows:

[0095]

[0096] Where X(t) is an n×1 dimensional random process, dX(t) is the random differential in the sense of Ito, f(X(t), t) is called the drift coefficient, which is an n-dimensional vector function of X(t), and is a slow-varying process that does not contain white noise; G(X(t), t) is called the diffusion matrix, which is an n×m matrix function; B(t) is an m-dimensional unit Wiener process, and its covariance matrix is ​​V; W(t) = dB(t) is a Gaussian white noise process with a spectral density of K, and is a fast-varying process;

[0097] Transform to integral form:

[0098]

[0099] The integrated energy system is affected by line failure rate, wind and solar output, and load changes. The system state X(t) is represented by X(t, x1, x2..., x n );

[0100] The process of constructing the Kolmogorov forward equation solution is as follows:

[0101] The random process X(t) is a Markov process. For any t1, t2, …, t n (t1<t2<…<t n ), whose conditional probability density is for all X(t1), X(t2), ..., X(t n )satisfy:

[0102] P[X(t n )≤x n |X(t n-1 )≤x n-1 ,…,X(t1)≤x1]=P[X(t n )≤x n |X(t n-1 )≤x n-1 ];

[0103] The probability density is of the form:

[0104] p(x n , t n |x n-1 , t n-1 ; ... ; x1, t1) = p(x n , t n |x n-1 , t n-1 );

[0105] p(x n , t n )=p(x n , t n |x n-1 , t n-1 )p(xn-1 , t n-1 |x n-2 , t n-2 )...p(x1,t1);

[0106] For the three moments t1<t<t2, we can get:

[0107] p(X(t2),t2|X(t),t;X(t1),t1)=p(X(t2),t2|X(t),t)p(X(t),t|

[0108] Let X(t1), X(t), X(t2) be simply x1, y, x2, and integrate Y to get:

[0109] p(x2,t2|x1,t1)=∫p(x2,t2|y,t)p(y,t|x1,t1);

[0110] Converting the integral equation to an equivalent differential equation gives the Kolmogorov forward equation:

[0111]

[0112] Among them, a j , b ij is the derivative moment, and the solution is as follows:

[0113]

[0114] Using the collected extreme disaster parameters, Monte Carlo simulation is used to solve a j , b ij ;

[0115]

[0116] δ is the Dirac function, indicating that at t = t0, the system is in the initial state X0 with probability 1;

[0117] Boundary conditions are divided into non-infinite and infinite boundaries. The infinite boundary can be expressed as:

[0118]

[0119] The normalization conditions are as follows:

[0120]

[0121] Since the Markov process is a stationary process,

[0122]

[0123] p(x) is the stationary probability density, a i and b ij Independent of time t; input the initial state of the power system operation, construct a stochastic differential equation for uncertainties including new energy and load, and substitute it into the Kolmogorov forward equation to calculate the probability distribution of the real-time state of the integrated energy system.

[0124] Furthermore, in step S6, the system resilience index is as follows:

[0125] Load reduction as an indicator of toughness:

[0126] R=Q w ;

[0127] In the formula, Q W Reduce the load on the system;

[0128] The load shedding amount multiplied by the weighted sum of the scenario probabilities is used as the resilience indicator:

[0129] R=∑P w Q w ;

[0130] Where P w is the scenario probability;

[0131] The integral of the system performance under disasters compared to the system performance without disasters is used as a resilience indicator:

[0132]

[0133] Where P(t) is the load of the system when an extreme disaster occurs, and P0(t) is the load of the system when no extreme disaster occurs;

[0134] Stage performance as an indicator of toughness:

[0135]

[0136] In the formula, R represents system resilience, t r Indicates the time when the system starts to recover; t r +T0 indicates the time when the system is fully restored.

[0137] The present invention can achieve the following technical effects:

[0138] The comprehensive energy system resilience assessment method based on the Kolmogorov equation provided by the present invention utilizes the Kolmogorov forward equation to learn system state characteristics from historical data, calculates drift coefficients and diffusion coefficients, constructs stochastic differential equations, and calculates the probability distribution of resilience indicators, thereby obtaining an uncertainty assessment of system resilience. BRIEF DESCRIPTION OF THE DRAWINGS

[0139] Figure 1 is a flow chart of a comprehensive energy system resilience assessment method based on the Kolmogorov equation provided according to an embodiment of the present invention;

[0140] Figure 2 is a structural diagram of a comprehensive energy system provided according to an embodiment of the present invention;

[0141] Figure 3 is a load diagram of each energy system provided according to an embodiment of the present invention;

[0142] Figure 4 is a probability distribution curve of system resilience provided according to an embodiment of the present invention. DETAILED DESCRIPTION

[0143] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. In the following description, the same modules are represented by the same reference numerals. In the case of the same reference numerals, their names and functions are also the same. Therefore, the detailed description thereof will not be repeated.

[0144] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and do not constitute a limitation of the present invention.

[0145] The embodiment of the present invention provides a comprehensive energy system resilience assessment method based on the Kolmogorov equation, the process of which is as follows: Figure 1 As shown, including:

[0146] 1. Collect historical data on extreme disasters in the system area and use Monte Carlo to simulate extreme disaster parameters.

[0147] 11. Collect historical data on regional extreme disasters.

[0148] Quantify extreme disaster parameters from both temporal and spatial characteristics, and collect key parameters that have a greater impact on multi-energy systems, such as wind speed, rainfall rate, ice accumulation, atmospheric turbidity and other parameters under extreme disasters.

[0149] 12. Monte Carlo simulation of the spatiotemporal distribution and parameter characteristics of extreme disasters.

[0150] Using the data collected above, we can calculate the parameters such as the types and intensities of disasters that may occur in the region through statistical calculations, and use Monte Carlo simulation to obtain the probability distribution of the characteristics of the impending extreme disasters.

[0151] 2. Build an integrated energy system of electricity, gas, heat and storage, and use energy hubs such as cogeneration units, gas boilers, gas turbines and power-to-gas to couple different energy systems together.

[0152] 21. Construct the power grid system, gas grid system and heat grid system separately.

[0153] The power grid system can adopt IEEE standard calculation examples, and it is necessary to additionally obtain the actual load demand, installed capacity of conventional units, wind and solar power output, and spatial distribution of topology, and must meet constraints such as power flow and power balance during operation; the gas network system selects the appropriate topology structure based on the actual area, and also needs to meet constraints such as pipeline flow during operation; the heat network system uses gas boilers to generate high-temperature steam for heating, and needs to meet constraints such as water supply and return water heat during operation.

[0154] The following constraints must be met when the power grid system is running:

[0155] I. AC power flow constraints:

[0156]

[0157] Wherein, formula (1) is the node power balance equation, P i and Q i are the active and reactive power injection of node i, respectively. Injection is positive and output is negative. V i and V j is the voltage amplitude at nodes i and j; G ij and B ij They are the network impedance matrix Y bus The real and imaginary parts of represent the elements of the admittance matrix; θ ij is the difference in voltage phase angle between nodes i and j; n is the total number of nodes in the network; equation (2) is the line power flow constraint, g ij and b ij The conductance and susceptance of the circuit ij, b sh,ij is the line shunt susceptance, P ij and Q ij are the active and reactive power flows of line ij respectively; Equation (3) is the voltage amplitude constraint, V i min 、V i max are the lower and upper limits of the voltage amplitude; Formula (4) is the power flow direction and capacity constraints, It is the upper limit of the active and reactive power transmission capacity of the line.

[0158] II. Unit climbing constraints:

[0159] P i,t -P i,t-1 ≤ΔP i,up

[0160] P i,t-1 -P i,t ≤ΔPi,down ; (5)

[0161] In the formula, ΔP i,up , ΔP i,down Indicates the upper and lower limits of the ramp rate of unit i per unit time.

[0162] III. Energy Storage Device Constraints.

[0163] The energy storage device only contains electrical energy storage devices and needs to meet the following constraints:

[0164]

[0165] u n,c,t +u n,d,t ≤1; (7)

[0166]

[0167] Wherein, formula (6) is the energy storage charge and discharge balance equation, P ESS,n,t The power input to the nth energy storage device (relative to the energy storage device); η n,c is the charging efficiency of the nth energy storage device; Formula (7) is the device charging and discharging state constraint, indicating that charging and discharging cannot be performed at the same time, u n,c,t is the charging status of the energy storage device n, 1 for yes and 0 for no; u n,d,t is the discharge state of the energy storage device n, 1 for yes, 0 for no. Formula (8) is the energy storage device power constraint, are the charging and discharging power of energy storage device n (the power consumed by the grid); SOC n,t is the storage capacity of the energy storage device; equations (9) and (10) are the charge and discharge rate constraints of the energy storage device, are the lower and upper limits of the charging power of the energy storage device n, respectively; are the lower and upper limits of the discharge power of the energy storage device n respectively; Formula (11) is the power constraint, are the lower and upper limits of the storage capacity of the energy storage device n respectively.

[0168] During the operation of the gas network, the following constraints need to be met:

[0169] I. Natural gas pipeline network equation:

[0170]

[0171]

[0172] Among them, formula (12) is the flow balance constraint, formula (13) is the pipeline flow constraint, formula (14) is the gas source gas production constraint, (15) and (16) are gas load loss constraints, formula (17) is the pipeline physical property constraint, and formula (18) is the node gas pressure constraint.

[0173] G ij,t For natural gas pipelines ij Natural gas flow rate in, m 3 / h, regulations: If the pipeline ij The flow direction of natural gas is: i→j, then G ij,t >0; is the gas output of the gas source; is the gas production of PtG; is the gas output of the gas storage device, where gas storage is negative and gas discharge is positive; is the demand for gas load; Pipeline ij Maximum and minimum values ​​of airflow; are the maximum and minimum values ​​of gas production of the gas source, respectively; is the gas consumption of GT; p i,t is the air pressure at the node; C ij For pipeline ij Correlation coefficient of physical properties; are the maximum and minimum values ​​of the air pressure at node i, respectively; is the set of nodes connected to node i; is the set of nodes connected to the gas source; is the set of nodes connected to the gas load; is the set of nodes connected to PtG; is a set of nodes connected to the gas storage device; is the set of nodes connected to GT; S G for gas pipeline assembly; A collection of natural gas pipelines with gas compressors. is the gas load loss.

[0174] II. Gas Compressor Constraints.

[0175] The gas compressor in the natural gas pipeline consumes electricity to resist the pressure drop on the pipeline to maintain the flow of natural gas. Its model is:

[0176]

[0177] Where, formula (19) is the gas pressure constraint of natural gas compression, Λ nis the compression coefficient of the gas compressor n, and it is stipulated that the gas flow in the pipeline containing the gas compressor is i→j; Formula (20) is the energy conversion constraint of the compressor, is the electrical power consumed by the gas compressor n, kW; 0-1 variable is the switch state of the gas compressor n, 1 is on, 0 is off; σ n is the power conversion factor of gas compressor n, kW / m 3 .

[0178] The heating network system needs to meet the following constraints during operation:

[0179] I. Heat medium balance constraints.

[0180] The flow of heat medium in the supply and return pipes needs to meet the following constraints:

[0181]

[0182] Wherein, equation (21) is the heat medium flow constraint. k,t is the mass flow rate of the heat medium in the heat supply / return pipe k; are the sets of heat supply / return pipes starting and ending at node i; N H It is a collection of nodes of the supply and heat recovery network.

[0183] II. Temperature Mixing.

[0184] The heat supply / return pipe network contains three types of nodes: ordinary nodes, branch nodes and intersection nodes. The mixed relationship of the temperature of each node is as follows:

[0185] Ordinary nodes:

[0186]

[0187] Branch Node:

[0188]

[0189] Junction Node:

[0190]

[0191] In the formula, are the temperatures of the heat medium flowing into and out of the heat supply / return pipe k, ℃; S H It is a collection of supply / return heat pipes;

[0192] At the same time, the temperature of the heat medium in the supply / return pipe needs to be limited to a reasonable range, that is:

[0193]

[0194] In the formula, They are the upper and lower limits of the heat medium temperature in the supply / return pipes respectively; usually, the temperature limit in the return pipe is lower than that in the supply pipe.

[0195] III. Pipeline heat loss and transmission delay.

[0196] The heat power carried by the heat medium in the heat supply / return pipe is proportional to the mass flow rate and temperature of the pipe, as follows:

[0197]

[0198] In the formula, are the inflow and outflow heat powers in the supply / recovery pipe k, kW; are the inflow and outflow temperatures in the supply / return heat pipe, respectively, ℃; c is the specific heat capacity, kJ / (kg·℃).

[0199] The heat power loss and delay of the supply / return heat pipe is:

[0200]

[0201] In the formula, μ loss is the heat loss coefficient per unit length; l k is the length of the supply / return heat pipe, m; τ k,t is the transmission delay time of pipeline k at time t.

[0202] IV. Heat exchange.

[0203] The heat power generated by the heat source is sent into the heat supply / recovery network through heat exchange, which is:

[0204]

[0205] In the formula, is the thermal power flowing into the heating pipe k1; is the thermal power flowing out of the heat recovery pipe k2; η ex,1 The heat exchange efficiency between the heat source and the heat supply / return transmission pipeline; is the thermal power generated by the heat source m; They are respectively a collection of supply and return heat pipes connected to the heat source.

[0206] Heat users (buildings) obtain heat from the supply / recovery network through heat exchange, including:

[0207]

[0208] Among them, equation (31) is a heat exchange equation, that is, the heat supply is equal to the heat demand.

[0209] 22. Construct the energy transmission equation of the energy hub between multi-energy systems. In a multi-energy system, the conversion of different energies needs to be carried out through the energy hub. According to the unified energy path theory, the energy conversion effect can be achieved by establishing the transmission equation.

[0210] I. Electric-thermal characteristic equation of cogeneration unit:

[0211]

[0212] In the formula, are the electric power and thermal power generated by the cogeneration unit connected to the grid node i and the heat network node j, kW; are the upper and lower limits of electric power and thermal power, kW respectively; are the slopes of the maximum and minimum steam inlet operating condition curves respectively; φ i,j are the slope and intercept of the minimum condensing operating condition curve, kW respectively.

[0213] II. Gas turbine.

[0214] Gas turbines generate electricity by consuming natural gas, which is a power source and gas load. Its model is:

[0215]

[0216] In the formula, The gas turbine connected to the grid node k and the natural gas node n generates electrical power, kW; for the efficiency of the gas turbine; is the gas consumed by the gas turbine, m 3 / h; are the upper and lower limits of the generated electric power respectively; is the conversion factor, which is used to unify the dimensions; Π is the higher calorific value of natural gas, kBtu / m 3 ; Φ is the energy conversion coefficient, kBtu / (kW·h).

[0217] III. Power-to-gas.

[0218] Power-to-gas generates natural gas by consuming electricity, with electricity as the load and gas as the source. The model is:

[0219]

[0220] In the formula, is the gas consumed by the gas turbine, m 3 / h; is the electrical power generated by the gas turbine, kW; η n ptg is the efficiency of the gas turbine; G nptg,max , G n ptg,min They are the upper and lower limits of gas generation respectively.

[0221] 3. The distribution and characteristics of extreme disasters are input into the multi-energy system as parameters, and the probability distribution of system failure scenarios is obtained.

[0222] 31. Comprehensively consider the impact of extreme disasters on the system, mainly including the following aspects:

[0223] I. Impact of extreme disasters on sources.

[0224] Extreme disasters mainly affect the power generation of renewable energy. Taking ice and snow disasters and sandstorms as examples, the impact of ice disasters on wind turbines is as follows:

[0225]

[0226] P ice,t =P(1-g1M turbine ); (41)

[0227] Where M turbine is the ice accumulation of the fan, k turbine is the ice accumulation coefficient, S is the precipitation rate, V is the wind speed, P is the wind turbine output before correction, P ice is the corrected fan output, and g1 is the ice accumulation factor.

[0228] The sandstorm mainly affects the power generation of photovoltaic units:

[0229]

[0230] The above formula is the photovoltaic output formula, where P pv Output for photovoltaic power; Y pv is the rated capacity of the photovoltaic array under standard test conditions; f pv is the photovoltaic derating factor, which is mainly related to the surface conditions of the photovoltaic array; G T is the solar radiation intensity currently received by the photovoltaic array, kW / m 2 ; G T,STC is the incident light intensity under standard test conditions, constant, 1kW / m 2 ; T c is the current surface temperature of the photovoltaic array, °C; T c,STC is the temperature of the photovoltaic array under standard test conditions, °C; α Pis the temperature correction coefficient, which is related to the material of the photovoltaic array. When a sandstorm comes, the increase in atmospheric turbidity will affect the intensity of solar radiation and surface temperature received by the current photovoltaic array. Among them, aerosols have a greater impact on photovoltaic output. When sandstorm weather occurs, the dust particles in the atmosphere will greatly increase, thereby reducing the output of the photovoltaic unit.

[0231] II. The impact of extreme disasters on the Internet.

[0232] When extreme disasters occur, the pressure on the lines exceeds the upper limit, and there may be a risk of disconnection. Taking ice disasters as an example, the impact of ice disasters on lines is as follows:

[0233]

[0234] In the formula, a w1 and b w1 are two threshold values ​​of ice wind load, L WI is the load under the combined action of wind and ice, P f is the line disconnection failure rate.

[0235] III. Impact of extreme disasters on the load.

[0236] Take ice disaster as an example. When ice disaster occurs, the temperature drops sharply, the transmission loss of heat network increases greatly, and the heat loss of buildings increases. The specific performance is shown in the following formula:

[0237]

[0238] In the formula, is the heat network transmission loss, is the heat loss of the building; T amb is the ambient temperature; K ij ,cn ij is the heat network transmission pipeline parameter, R i is the thermal resistance of the building, is the air transmission parameter.

[0239] 4. Construct the Ito stochastic differential equation for the multi-energy system, decompose the system change process into the system slow-changing process and the system fast-changing process, perform mean clustering to divide the scenario probability distribution of the slow-changing process, and deduce its characteristics through historical data for the fast-changing process. Assume that the fast-changing process is Gaussian white noise, calculate its expectation and variance, and obtain the probability distribution of the fast-changing process. Construct the Kolmogorov forward equation (FPK), substitute the system operation data sample obtained by the above calculation into the FPK equation, calculate the Kolmogorov entropy value, and obtain the drift coefficient and diffusion coefficient.

[0240] 41. Construct Ito stochastic differential equation.

[0241] The multi-energy system under the influence of extreme disasters is described as the following stochastic differential equation:

[0242]

[0243] Where X(t) is an n×1 dimensional random process, dX(t) is the random differential in the sense of Ito, f(X(t), t) is called the drift coefficient, which is an n-dimensional vector function of X(t), and is a slow-varying process that does not contain white noise. G(X(t), t) is called the diffusion matrix, which is an n×m matrix function. B(t) is an m-dimensional unit Wiener process, and its covariance matrix is ​​V; W(t) = dB(t) is a Gaussian white noise process with a spectral density of K, which is a fast-varying process.

[0244] Write the above formula in integral form:

[0245]

[0246] The above is the Ito random state evolution equation. Since the multi-energy system is affected by multiple variables including line failure rate, wind and solar output, and load changes, the system state X(t) should be represented as X(t, x1, x2..., x n ) and substitute it into the above formula to solve it.

[0247] 42. Solve the probability distribution of system states.

[0248] Due to the uncertainty of the system operation state caused by extreme disasters, it is necessary to further analyze the Ito stochastic differential equation. The above differential equations have been divided into slow-changing processes and fast-changing processes, that is, mathematically expressed as, the solution of the Ito stochastic differential equation is a continuous-time Markov process, which is a diffusion process. According to the theory of stochastic dynamics, the probability density of the solution of the system state satisfies the Kolmogorov forward equation (FPK), and the probability density function of the solution can be obtained by solving the Kolmogorov forward equation. The deduction is as follows:

[0249] A random process X(t) is a Markov process. For any t1, t2, ..., t n (t1<t2<…<t n ), whose conditional probability density is for all X(t1), X(t2), ..., X(t n )satisfy:

[0250] P[X(t n )≤x n |X(t n-1 )≤x n-1 ,…,X(t1)≤x1]=P[X(t n )≤x n |X(t n-1 )≤x n-1]; (48)

[0251] The above formula is called the "no aftereffect" of the Markov process, and its probability density is in the form of:

[0252] p(x n , t n |x n-1 , t n-1 ; ... ; x1, t1) = p(x n , t n |x n-1 , t n-1 ); (49)

[0253] So we can get:

[0254] p(x n , t n )=p(x n , t n |x n-1 , t n-1 )p(x n-1 , t n-1 |x n-2 , t n-2 )...p(x1,t1);(50)

[0255] Solve the Kolmogorov forward equation for a Markov process.

[0256] For the three moments t1<t<t2, we can get:

[0257] p(X(t2),t2|X(t),t;X(t1),t1)=p(X(t2),t2|X(t),t)p(X(t),t|

[0258] Let X(t1), X(t), X(t2) be abbreviated as x1, y, x2. Integrating the above formula with respect to Y yields:

[0259] p(x2,t2|x1,t1)=∫p(x2,t2|y,t)p(y,t|x1,t1); (52)

[0260] Converting the integral equation to an equivalent differential equation gives the Kolmogorov forward equation:

[0261]

[0262] Among them, a j , b ij is the derivative moment, and the solution is as follows:

[0263]

[0264] Using the collected extreme disaster parameters, Monte Carlo simulation is used to solve a j 、b ij Solving this equation also requires initial conditions, boundary conditions, and normalization conditions. When the initial state is fixed, the initial conditions are:

[0265]

[0266] δ is the Dirac function, which means that at t=t0, the system is in the initial state X0 with probability 1.

[0267] Boundary conditions are divided into non-infinite and infinite boundaries. The infinite boundary can be expressed as:

[0268]

[0269] The normalization conditions are as follows:

[0270]

[0271] In order to solve the probability density distribution when the system reaches a stable state, since the Markov process is a stable process, The equation now becomes:

[0272]

[0273] The above equation is called the stationary Kolmogorov forward equation or the simplified Kolmogorov forward equation, p(x) is the stationary probability density, a i and b ij It is independent of time t. The initial state of the power system operation is used as input, and a stochastic differential equation is constructed to describe the uncertainties including new energy and load, which is substituted into the FPK equation to calculate the probability distribution of the real-time state of the system.

[0274] 5. Input the data of the upcoming predicted extreme disasters into the FPK equation to solve the probability distribution of emergency load shedding of the system.

[0275] The probability distribution of the system's operating state is obtained from the above steps. The system parameters and extreme disaster parameters are input. The system starts to operate in a certain initial state. The dispatching strategy of the multi-energy system is input. When the dispatching strategy cannot meet the system load demand, emergency load shedding is performed. According to the probability distribution of the system's operating state in a certain initial state, the probability distribution of emergency load shedding can be obtained.

[0276] 6. Construct system resilience indicators, draw resilience curves, and output system resilience uncertainty assessment.

[0277] 61. Propose toughness index. The toughness index can be the maximum load shedding capacity of the system, the system recovery rate, and the area of ​​the system toughness curve.

[0278] The main toughness index R can be expressed by the following formula:

[0279] 1) Directly use the load reduction as the toughness index:

[0280] R=Q w ; (60)

[0281] Among them, Q W Reduce the load on the system.

[0282] 2) The weighted sum of load reduction multiplied by scenario probability is used as the resilience index:

[0283] R=∑P w Q w ; (61)

[0284] Among them, P w is the scene probability.

[0285] 3) The integral of the system performance under disasters and the system performance without disasters is used as the resilience index:

[0286]

[0287] Among them, P(t) is the load of the system when an extreme disaster occurs, and P0(t) is the load of the system when no extreme disaster occurs.

[0288] 4) Performance at a certain stage as an indicator of resilience:

[0289]

[0290] Where R represents system resilience, t r Indicates the time when the system starts to recover; t r +T0 indicates the time when the system is fully restored.

[0291] 62. Obtain uncertainty assessment of system resilience. Substitute the probability distribution of the system's operating state under extreme disasters, use the probability distribution of system load shedding obtained in step 5, bring it into the resilience index framework, and obtain uncertainty assessment of the system resilience index.

[0292] This example uses the "14+7+6" electricity-gas-heat integrated energy system example to simulate and verify the effectiveness of the proposed method. Figure 2 shown.

[0293] The load of each energy system is as follows: Figure 3Taking ice and snow disasters as an example, the probability of line disconnection of each line is obtained according to the present invention as follows:

[0294] P f =[0.09,0.35,0.03,0.64,0.19,0.47,0.23,0.37,0.39,0.14,0.69,0.55,0.66,0.63,0.43,0.65,0.07,0.14,0.04,0.23];

[0295] Taking the total amount of load lost by the multi-energy system within a certain period of time as the resilience index, the system resilience index is calculated with and without new energy installed, and the probability distribution curve of the resilience of the multi-energy system is obtained, as shown in Figure 2. Figure 4 shown.

[0296] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "example", "specific example", or "some examples" etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification and the features of the different embodiments or examples, without contradiction.

[0297] Although the embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and are not to be construed as limitations of the present invention. A person skilled in the art may change, modify, replace and vary the above embodiments within the scope of the present invention.

[0298] The above specific implementations of the present invention do not constitute a limitation on the protection scope of the present invention. Any other corresponding changes and modifications made based on the technical concept of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A comprehensive energy system resilience assessment method based on the Kolmogorov equation, characterized in that: The steps include: S1. Collect extreme disaster parameters in the system area, use Monte Carlo simulation to simulate extreme disaster parameters, and obtain the probability distribution of the characteristics of the extreme disaster that will occur; S2. Build an integrated energy system of electricity, gas, heat and storage, using cogeneration units, gas boilers, gas turbines and power-to-gas energy hubs to couple the energy systems; S3, input the distribution and characteristics of extreme disasters into the integrated energy system as parameters, calculate the impact of extreme disasters on the integrated energy system, and obtain the probability distribution of system failure scenarios; S4. Construct the Ito stochastic differential equation of the integrated energy system, decompose the change process of the integrated energy system into slow-changing process and fast-changing process to solve the probability distribution, construct the Kolmogorov forward equation to solve it, calculate the Kolmogorov entropy value, and obtain the drift coefficient and diffusion coefficient; S5, inputting the predicted extreme disaster data into the Kolmogorov forward equation to solve the probability distribution of emergency load shedding of the integrated energy system; S6. Construct system resilience indicators, draw resilience curves, and output uncertainty assessment of comprehensive energy system resilience.

2. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 1 is characterized in that: In step S1, the system area extreme disaster parameters include wind speed, rainfall rate, ice accumulation, and atmospheric turbidity parameters, and the probability distribution of the characteristics of the extreme disaster to be occurred is obtained through Monte Carlo simulation.

3. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 1 is characterized in that: In step S2, the construction of the electricity-gas-heat-storage integrated energy system includes constructing a power grid system, a gas grid system, a heat grid system and an energy transmission equation; The power grid system operation constraints include AC power flow constraints, unit ramp constraints, and energy storage device constraints, which are as follows: The AC power flow constraints are: Where P i and Q i are the active and reactive power injection of node i, respectively. Injection is positive and output is negative. V i and V j is the voltage amplitude at nodes i and j; G ij and B ij They are the network impedance matrix Y bus The real and imaginary parts of represent the elements of the admittance matrix; θ ij is the difference in voltage phase angle between nodes i and j; n is the total number of nodes in the network; g ij and b ij The conductance and susceptance of the circuit ij, b sh,ij is the line shunt susceptance, P ij and Q ij are the active and reactive power flows of line ij respectively; are the lower and upper limits of the voltage amplitude; The upper limit of the active and reactive transmission capacity of the line; The unit climbing constraints: P i,t -P i,t-1 ≤ΔP i,up ; P i,t-1 -P i,t ≤ΔP i,down ; In the formula, ΔP i,up , ΔP i,down Indicates the upper and lower limits of the ramp rate of unit i per unit time; The electrical energy storage device constraints: in n,c,t +in n,d,t ≤1; Where P ESS,n,t is the power input to the nth energy storage device; η n,c is the charging efficiency of the nth energy storage device; u n,c,t is the charging status of the energy storage device n, 1 for yes and 0 for no; u n,d,t is the discharge state of the energy storage device n, 1 for yes and 0 for no; are the charging and discharging power of energy storage device n; SOC n,t The amount of energy stored in the energy storage device; are the lower and upper limits of the charging power of the energy storage device n, respectively; are the lower and upper limits of the discharge power of the energy storage device n respectively; are the lower and upper limits of the storage capacity of the energy storage device n respectively.

4. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 3 is characterized in that: In step S2, the operation of the gas network system includes the natural gas pipeline network equation and the gas compressor constraints, which are as follows: The natural gas pipeline network equation: In the formula, G ij,t is the natural gas flow rate in the natural gas pipeline ij, m 3 / h, it is stipulated that if the flow direction of natural gas in pipeline ij is: i→j, then G ij,t >0; is the gas output of the gas source; is the gas production of PtG; is the gas output of the gas storage device, where gas storage is negative and gas discharge is positive; is the demand for gas load; are the maximum and minimum values ​​of the airflow in pipeline ij, respectively; are the maximum and minimum values ​​of gas production of the gas source, respectively; is the gas consumption of GT; p i,t is the air pressure at the node; C ij is the correlation coefficient with the physical characteristics of pipeline ij; are the maximum and minimum values ​​of the air pressure at node i, respectively; is the set of nodes connected to node i; is the set of nodes connected to the gas source; is the set of nodes connected to the gas load; is the set of nodes connected to PtG; is a set of nodes connected to the gas storage device; is the set of nodes connected to GT; S G for natural gas pipeline assembly; A collection of natural gas pipelines with gas compressors. is the gas load loss; The gas compressor constraints: In the formula, Λ n is the compression factor of the gas compressor n, and it is stipulated that: the gas flow in the pipeline containing the gas compressor is i→j; is the electrical power consumed by the gas compressor n, kW; 0-1 variable is the switch state of the gas compressor n, 1 is on, 0 is off; σ n is the power conversion factor of gas compressor n, kW / m 3 .

5. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 3 is characterized in that: In step S2, the operation of the heat network system includes heat medium balance constraints, temperature mixing, pipeline heat loss and transmission delay, and heat exchange, as follows: The heat medium balance constraint: In the formula, h k,t is the mass flow rate of the heat medium in the heat supply / return pipe k; are the sets of heat supply / return pipes starting and ending at node i; N H It is the node set of the heat supply and heat recovery network; The temperature mix: The heat supply / return pipe network contains three types of nodes: ordinary nodes, branch nodes and intersection nodes. The mixed relationship of the temperature of each node is as follows: Ordinary nodes: Branch Node: Junction Node: In the formula, are the temperatures of the heat medium flowing into and out of the heat supply / return pipe k, ℃; S H It is a collection of supply / return heat pipes; Temperature range of heat medium in supply / return pipe: In the formula, They are the upper and lower limits of the heat medium temperature in the supply / return pipes, respectively. The limit of the temperature in the return pipe is lower than that in the supply pipe. The pipeline heat loss and transmission delay: The heat power carried by the heat medium in the heat supply / return pipe is proportional to the mass flow rate and temperature of the pipe. In the formula, are the inflow and outflow heat powers in the supply / recovery pipe k, kW; are the inflow and outflow temperatures in the supply / return heat pipe, respectively, ℃; c is the specific heat capacity, kJ / (kg·℃); The heat power loss and delay of the supply / return heat pipe is: In the formula, μ loss is the heat loss coefficient per unit length; l k is the length of the supply / return heat pipe, m; τ k,t is the transmission delay time of pipeline k at time t; The heat exchange: The heat power generated by the heat source is sent to the heat supply / recovery network through heat exchange. In the formula, is the thermal power flowing into the heating pipe k1; is the thermal power flowing out of the heat recovery pipe k2; η ex,1 The heat exchange efficiency between the heat source and the heat supply / return transmission pipeline; is the thermal power generated by the heat source m; They are respectively a collection of heat supply and heat return pipes connected to the heat source; Heat users obtain heat from the supply / recovery network through heat exchange.

6. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 3 is characterized in that: In step S2, the energy transfer equation includes the electric-thermal characteristic equation of the cogeneration unit, the gas turbine model, and the power-to-gas model, which are as follows: The electric-thermal characteristic equation of the cogeneration unit: In the formula, are the electric power and thermal power generated by the cogeneration unit connected to the grid node i and the heat network node j, kW; are the upper and lower limits of electric power and thermal power, kW respectively; are the slopes of the maximum and minimum steam inlet operating condition curves respectively; φ i,j are the slope and intercept of the minimum condensing condition curve, kW; The gas turbine model: Gas turbines generate electricity by consuming natural gas, providing power, gas loads, In the formula, The gas turbine connected to the grid node k and the natural gas node n generates electrical power, kW; for the efficiency of the gas turbine; is the gas consumed by the gas turbine, m 3 / h; are the upper and lower limits of the generated electric power respectively; θ = Π / Φ is the conversion factor, which serves to unify the dimensions; Π is the higher calorific value of natural gas, kBtu / m 3 ; Φ is the energy conversion coefficient, kBtu / (kW·h); The power-to-gas model: Power-to-gas generates natural gas by consuming electricity, which is the power load, gas source, In the formula, is the gas consumed by the gas turbine, m 3 / h; is the electrical power generated by the gas turbine, kW; η n ptg is the efficiency of the gas turbine; G n ptg,max , G n ptg,min They are the upper and lower limits of gas generation respectively.

7. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 1 is characterized in that: In step S3, the impact of extreme disasters on the integrated energy system is calculated, including the impact of extreme disasters on the source, the impact of extreme disasters on the grid, and the impact of extreme disasters on the load, as follows: The impact of extreme disasters on sources: The impact of ice disasters on wind turbines is as follows: P ice,t =P(1-g1M turbine ); Where M turbine is the ice accumulation of the fan, k turbine is the ice accumulation coefficient, S is the precipitation rate, V is the wind speed, P is the wind turbine output before correction, P ice is the corrected fan output, g1 is the ice accumulation factor; The impact of sandstorms on the power generation of photovoltaic units is as follows: Where P pv For photovoltaic output; Y pv is the rated capacity of the photovoltaic array under standard test conditions; f pv is the photovoltaic derating factor, which is related to the surface conditions of the photovoltaic array; G T is the solar radiation intensity currently received by the photovoltaic array, kW / m 2 ; G T,STC is the incident light intensity under standard test conditions, constant, 1kW / m 2 ; T c is the current surface temperature of the photovoltaic array, °C; T c,STC is the temperature of the photovoltaic array under standard test conditions, °C; α P is the temperature correction factor; The impact of extreme disasters on the network: The impact of the ice disaster on the line is as follows: In the formula, a w1 and b w1 are two threshold values ​​of ice wind load, L WI is the load under the combined action of wind and ice, P f is the line disconnection failure rate; The impact of extreme disasters on the load: The impact of the ice disaster on the Netherlands is as follows: In the formula, is the heat network transmission loss, is the heat loss of the building; T amb is the ambient temperature; k ij ,cn ij is the heat network transmission pipeline parameter, R i is the thermal resistance of the building, is the air transmission parameter.

8. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 1 is characterized in that: In step S4, the Ito stochastic differential equation is as follows: Where X(t) is an n×1 dimensional random process, dX(t) is the random differential in the sense of Ito, f(X(t), t) is called the drift coefficient, which is an n-dimensional vector function of X(t), and is a slow-varying process that does not contain white noise; G(X(t), t) is called the diffusion matrix, which is an n×m matrix function; B(t) is an m-dimensional unit Wiener process, and its covariance matrix is ​​V; W(t) = dB(t) is a Gaussian white noise process with a spectral density of K, and is a fast-varying process; Transform to integral form: The integrated energy system is affected by line failure rate, wind and solar output, and load changes. The system state X(t) is represented by X(t, x1, x2…, x n ); The process of constructing the Kolmogorov forward equation solution is as follows: The random process X(t) is a Markov process. For any t1, t2, …, t n (t1<t2<…<t n ), whose conditional probability density is for all X(t1), X(t2), ..., X(t n )satisfy: P[X(t n )≤x n |X(t n-1 )≤x n-1 ,…,X(t1)≤x1]=P[X(t n )≤x n |X(t n-1 )≤x n-1 ]; The probability density is of the form: p(x n ,t n |x n-1 ,t n-1 ;…;x1,t1)=p(x n ,t n |x n-1 ,t n-1 ); p(x n ,t n )=p(x n ,t n |x n-1 ,t n-1 )p(x n-1 ,t n-1 |x n-2 ,t n-2 )…p(x1,t1); For the three moments t1<t<t2, we can get: p(X(t2),t2|X(t),t;X(t1),t1)=p(X(t2),t2|X(t),t)p(X(t),t| Let X(t1), X(t), X(t2) be simply x1, y, x2, and integrate Y to get: p(x2,t2|x1,t1)=∫p(x2,t2|y,t)p(y,t|x1,t1); Converting the integral equation to an equivalent differential equation gives the Kolmogorov forward equation: Among them, a j , b ij is the derivative moment, and the solution is as follows: Using the collected extreme disaster parameters, Monte Carlo simulation is used to solve a j 、b ij ; δ is the Dirac function, indicating that at t = t0, the system is in the initial state X0 with probability 1; Boundary conditions are divided into non-infinite and infinite boundaries. The infinite boundary can be expressed as: The normalization conditions are as follows: Since the Markov process is a stationary process, p(x) is the stationary probability density, a i and b ij Independent of time t; input the initial state of the power system operation, construct a stochastic differential equation for uncertainties including new energy and load, and substitute it into the Kolmogorov forward equation to calculate the probability distribution of the real-time state of the integrated energy system.

9. The comprehensive energy system resilience assessment method based on the Kolmogorov equation according to claim 1 is characterized in that: In step S6, the system resilience index is as follows: Load reduction as an indicator of toughness: R=Q w ; In the formula, Q W Reduce the load on the system; The load shedding amount multiplied by the weighted sum of the scenario probabilities is used as the resilience indicator: R=∑P w Q w ; Where P w is the scenario probability; The integral of the system performance under disasters compared to the system performance without disasters is used as a resilience indicator: In the formula, P(t) is the load of the system when an extreme disaster occurs, and P0(t) is the load of the system when no extreme disaster occurs; Stage performance as an indicator of toughness: In the formula, R represents system resilience, t r Indicates the time when the system starts to recover; t r +T0 indicates the time when the system is fully restored.

Images