A comprehensive energy system planning method considering the credibility of demand response

Abstract

The present invention discloses a comprehensive energy system planning method considering the credibility of demand response. First, the comprehensive energy load data, wind power data, generator data, natural gas source data, power grid data, natural gas network parameter data and planning cost parameter data on the energy hub side are obtained; then a demand response model is established, and a comprehensive energy system planning model considering the energy hub is established. The lower layer of the comprehensive energy system is the transmission line and the gas transmission line, and the upper layer is the energy hub; finally, the mixed integer linear programming is solved to obtain the corresponding planning results. The present invention establishes a model to describe the performance of customer demand response based on the Z number, which solves the problem of incomplete uncertainty modeling in traditional methods. The present invention can simultaneously consider the impact of uncertain factors and information credibility on customer demand response performance, so as to make a more accurate and comprehensive estimate of the actual situation.

Description

Technical Field

[0001] The invention belongs to the field of integrated energy system planning and operation, and proposes an integrated energy system planning method taking into account the credibility of demand response. Background Art

[0002] With the rapid economic development, human society is facing an increasingly serious energy depletion crisis and environmental pollution problems. On the one hand, the non-renewable nature of primary energy requires the accelerated development of new energy and the improvement of energy efficiency; on the other hand, the traditional energy system that only manages single energy sources such as electricity, gas, and heat in isolation greatly limits the comprehensive and coordinated utilization of energy and is not conducive to the improvement of energy efficiency. In this context, the concept of energy hub (EH) came into being.

[0003] The concept of energy hub comes from the research project "Vision of Future Energy Network" conducted by ETH Zurich in 2007. It is defined as an input-output port model used to describe the exchange and coupling relationship between energy, load and network in the integrated energy system. The coupling matrix describing the input energy and output load ports can briefly represent the various coupling relationships such as conversion, storage and transmission between various forms of energy such as electricity, heat and gas, and plays an important role in the planning and operation research of the integrated energy system.

[0004] At the same time, with the rapid development of smart grids, demand response (DR) technology has gradually become a research hotspot in the power industry. Users participating in demand response can adjust their consumption patterns during peak hours to alleviate power shortages based on the electricity price signals or incentive mechanisms released by the power market. Typical incentive-based demand-side responses include the PJM capacity market in the United States and the ISO England power market. The PJM market is dominated by load reduction suppliers, who obtain load reduction amounts by signing contracts with users. Unlike traditional power generation resources, the capacity benefits of DR are actually mined from changes in customer energy consumption patterns. Therefore, the extent to which DR enhances system reliability depends not only on the technical characteristics of the end-user load appliances, but also on the behavioral preferences of individuals on the demand side.

[0005] Through the analysis of PJM and ISO England customer demand response performance data, it can be seen that customers do not always fulfill their contract commitments 100%. The present invention takes into account the Z number to describe the user's demand response performance, considers the uncertainty of the user's demand response performance and the reliability of information, and combines the comprehensive energy system planning method considering the energy hub. On this basis, the present invention proposes a comprehensive energy system planning method considering the credibility of demand response. Summary of the invention

[0006] The purpose of this invention is to analyze the uncertainty of customer demand response performance and propose a modeling method for integrated energy system planning that considers the credibility of demand response. First, based on the Z number, this method establishes the final probability density model of customer demand response performance; then, according to the probability density function, the Monte Carlo method is used for sampling, the performance is output in n states and the average is obtained, and finally, combined with the integrated energy system planning model considering the energy hub, a comprehensive energy system planning method that considers the credibility of demand response is proposed.

[0007] The present invention is based on the following steps:

[0008] According to the uncertain performance of customer demand response, the present invention proposes a comprehensive energy system planning method considering the credibility of demand response, and the steps are as follows:

[0009] Step (1), obtaining comprehensive energy load data, wind power data, generator data, natural gas source data, power grid data, natural gas grid parameter data and planning cost parameter data on the energy hub side.

[0010] Step (2): Establish a demand response model.

[0011] Step (3) establishes a comprehensive energy system planning model considering energy hubs. The lower layer of the comprehensive energy system is the power transmission line and gas transmission line, and the upper layer is the energy hub.

[0012] Step (4), solve the mixed integer linear programming and obtain the corresponding planning results.

[0013] Furthermore, the specific method of step (1) is as follows:

[0014] Obtain the cooling, heating and electricity load data of typical days in summer, transition season and winter on the energy hub side, the wind power output data on three typical days, the output upper limit, ramp rate and unit output cost data of the generator, the output upper limit and unit output cost data of the natural gas source, the power upper limit, reactance and investment cost data of the power grid, the flow upper limit, pipeline coefficient and investment cost data of the natural gas network, DR incentives, DR penalties and unit cost data of wind abandonment.

[0015] Furthermore, the specific steps of step (2) are as follows:

[0016] Step (2.1), use Z number to establish customer demand response performance uncertainty model. Introduce variable u b The customer's demand response performance is expressed as the ratio of the customer's actual load reduction to the amount contracted by the customer. The mathematical expression is:

[0017]

[0018] Among them, u b represents the demand response performance of the customer at node b, represents the actual load reduction of the customer on node b at time t, represents the load reduction contracted by the customer at node b at time t. b Lower customers compared to u b Customers with higher load shedding will be more likely to comply with their contractual commitments. b =0, which means the customer does not perform any load reduction at all.

[0019] u b The value of is uncertain. To take this uncertainty into account, a method based on the Z number is used. According to the Z number method, for any uncertain variable, its impact can be expressed as an ordered pair of fuzzy numbers. To indicate that Is the restriction part, indicating that u b The assumptions of restriction, The reliability part indicates Therefore, the mathematical formula of the customer's demand response performance Z number can be expressed as:

[0020]

[0021] in,

[0022]

[0023]

[0024] Among them, x1-x4 and x′1-x′4 are the trapezoidal membership functions defined and triangle membership function In order to convert the model into a mixed integer linear programming for easy solution, the fuzzy expectation method is used to convert the Z number in the model into a conventional fuzzy quantity, and then the centroid method is used to convert the obtained fuzzy variables into corresponding probability variables. The implementation steps are as follows:

[0025] (a) Convert the reliability part into a clear number, namely the reliability weight:

[0026]

[0027] Where β is the reliability weight, Representation function The integral value within its domain is Representation function The value of the integral within its domain.

[0028] (b) Add the obtained reliability weight to the restricted part, and the weighted Z number can be expressed as:

[0029]

[0030] in, is the weighted Z number, It is weighted Membership function.

[0031] (c) Convert the irregular fuzzy number (weighted fuzzy number) into a regular fuzzy number:

[0032]

[0033] in, is the regular fuzzy number of the weighted fuzzy number.

[0034] (d) Determined by the centroid method The corresponding equivalent probability density function is:

[0035]

[0036] where f(x) is the corresponding probability density function.

[0037] Step (2.2), simplify the model. Use the Monte Carlo method to sample the corresponding probability density function f(x) n times, and finally take the average to obtain u b For the trapezoidal probability density function f(x) of the above demand response performance, the parameters are x1, x2, x3, and x4. Discretize it first and insert the number x between x1 and x2. 1,2 , insert the number x between x2 and x3 2,3 , insert the number x between x3 and x4 3,4 . Divide x1-x4 into six segments evenly. The interval [x1,x 1,2 ]、[x 1,2 ,x2], [x2,x 2,3 ]、[x 2,3 ,x3]、 [x3,x 3,4 ]、[x 3,4 ,x4] The corresponding probabilities P are P1-P6 respectively. Monte Carlo sampling is performed. If the generated random number r is between [0,P1], then the corresponding u b for If the random number r is between [P1, P1+P2], then u b for If the random number r is between [P1+P2,P1+P2+P3], then u b for If the random number r is between [P1+P2+P3,P1+P2+P3+P4], then u b for If the random number r is between [P1+P2+P3+P4,P1+P2+P3+P4+P5], then u b for If the random number r is between [P1+P2+P3+P4+P5,P1+P2+P3+P4+P5+P6], then u b for As shown below:

[0038]

[0039] Perform the Monte Carlo sampling process as shown above n times, take the average of the results, and obtain the customer demand response performance u b .

[0040] Furthermore, the specific steps of step (3) are as follows:

[0041] Step (3.1): Establish an energy hub model: In the energy hub model, there are three types of component models:

[0042] Energy conversion components include combined heat and power (CHP), compression chillers (CERG) and electric heat pumps (HP), energy storage components include heat storage (HS) and cold storage (CS) devices, and renewable energy components include wind power (WT).

[0043] (a): Energy conversion element:

[0044] The input-output relationship of the energy hub is expressed as:

[0045]

[0046] Among them, I b,e (t) and I b,g (t) is the electric power and natural gas flow transmitted from the lower grid to the energy hub connected to node b at time t, P b,wt (t) is the output power of the wind power plant at the energy hub connected to node b at time t, and is the efficiency of natural gas conversion into electricity and heat, η HP and η CE are the efficiency of electric heating by electric heat pump and electric cooling by refrigerant, P b,hs (t) and P b,cs (t) are the charging power of heat storage and cold storage in the energy hub connected to node b at time t. b,1 (t), v b,2 (t), and (1-v b,1 (t)-vb,2 v′ is the dispatch factor for the electric power transmitted from the lower grid to the energy hub connected to node b at time t, which is allocated to the electric load, compression refrigerator and electric heat pump; b,1 (t), v′ b,2 (t), and (1-v′ b,1 (t)-v′ b,2 (t)) is the dispatch factor for distributing the electric power output of the cogeneration in the energy hub connected to node b at time t to the electric load, compression refrigerator and electric heat pump; v″ b,1 (t), v b,2 (t)″, and (1-v″ b,1 (t)-v″ b,2 (t)) is the dispatch factor for allocating the wind power output power in the energy hub connected to node b at time t to the electrical load, compression refrigerator and electric heat pump.

[0047] (b): Energy storage element:

[0048] For the energy storage element, the model is as follows:

[0049]

[0050] Among them, S b,j (t) represents the state of charge (SOC) of the energy hub system component j connected to node b at time t; and Represent the charging and discharging efficiency of system component j respectively; Δt represents the time span of a single time period; and They represent the input and output power of the system element j connected to the energy hub of node b at time t in the energy form k (k = c is the energy form of cold, k = h is the energy form of heat, and k = e is the energy form of electricity). For the heat storage element, if P b,hs (t)<0, then otherwise, The same applies to cold storage elements. j Table system component j storage capacity. sta_cha b (t) and sta_dis b (t) are all 0-1 variables, representing the charging and discharging states of the energy storage element in the energy hub connected to node b at time t.

[0051] (c): Renewable energy components:

[0052] Distributed renewable energy components are zero-input single-output components. The output of wind turbines is shown in the formula:

[0053]

[0054] Among them, P r is the rated output of the wind turbine, v is the actual wind speed, and v i is the cut-in wind speed, v0 is the cut-out wind speed, v r is the rated wind speed.

[0055] Step (3.2): Establish the energy network model of the integrated energy system:

[0056] (a): The DC power flow model used in the power grid is as follows:

[0057]

[0058]

[0059] Among them: α represents the form of energy transmitted from the lower grid to the energy hub, α = e represents electric power, and α = g represents natural gas flow. represents the power on line L at time t; θ L (t) and θ L ′(t) represents the phase of the starting point and the end point of the transmission line L at time t; x L Represents the reactance of line L; LE α and LC α Represent the collection of built and to-be-built lines respectively; is a 0-1 variable indicating whether the line L to be built is constructed or not, M e is a constant of 1000. Equation (14) uses the Big M method to establish the power flow formula for the newly built line.

[0060] (b): Gas network model:

[0061] The flow rate of natural gas pipelines is related to many factors such as the diameter, temperature, and pressure of the pipeline, and is nonlinear. There are many formulas to describe the flow equation of natural gas pipelines, including the American Gas Association formula, the Colebrook-White formula, the Panhandle A formula, the Weymouth formula, etc. The Weymouth steady-state model is used to describe the flow rate of natural gas pipelines, that is, the flow rate is only related to the pressure at both ends of the pipeline, and flows from the node with high pressure to the node with low pressure. For the natural gas pipeline L, the first and last nodes are n and n′ respectively. The specific expression of the natural gas flow rate and the node pressure flowing through the pipeline is:

[0062]

[0063] in, Indicates the natural gas flow in the pipeline; c n-n′ is the pipeline transmission parameter, which is related to the pipeline length, diameter and temperature; pn and p n′ denote the gas pressure of nodes n and n′ in the natural gas grid respectively; sgn(p n ,p n′ ) is a direction variable, representing the flow direction of natural gas in the pipeline. n ,p n′ )=1, natural gas flows from node n to node n′, sgn(p n ,p n′ )=-1, natural gas flows from node n′ to node n.

[0064] In order to reduce the difficulty of solving the natural gas subsystem model and the complexity of its coupling with the power subsystem, the nonlinear natural gas pipeline flow equation is linearized by piecewise linearization, and then the nonlinear model is transformed into a mixed integer linear programming model. The nonlinear term in the natural gas subsystem model, that is, the square of the node pressure in equation (15), is converted into a mixed integer linear programming model by introducing a new variable p′ and setting p′=p 2 By removing the square term, equation (15) becomes equation (16):

[0065]

[0066] The square term on the left side of equation (16) Perform piecewise linearization, the specific formula is equation (17)-(21), and we get:

[0067]

[0068]

[0069]

[0070]

[0071]

[0072] Among them, N represents the total number of linear segments. The larger the value of N, the more linear segments there are and the higher the accuracy. l represents the current number of segments. Represents the slope of a linear piecewise function; represents the intercept of the linear piecewise function on the y-axis, It is a 0-1 variable, indicating whether the point in question falls within the linear segmentation interval. 1 means it falls within the interval, and 0 means it does not fall within the interval. Represents the length of the linear piecewise function.

[0073] Furthermore, the specific method of step (4) is as follows:

[0074] The planning model constraints of the system are divided into node energy conservation constraints, power / gas source constraints, line power constraints, three types of component model constraints, and power grid and gas grid model constraints, as follows:

[0075] (a): Node energy conservation constraint:

[0076] An energy hub can be viewed as a load of multiple energy forms on a lower grid node.

[0077]

[0078] in, represents the output of the generator or natural gas source at node b at time t, L1 b and L2 b They represent the set of lines starting and ending at node b, respectively. represents the output power of the renewable energy source connected to node b at time t.

[0079] (b): Power supply / gas source constraints:

[0080] Formula (23) is the upper and lower limit constraints of generator / natural gas source output; Formula (24) is the generator / natural gas source output ramp constraint; Formula (25) is the centralized renewable energy output constraint:

[0081]

[0082]

[0083]

[0084] in, and They represent the lower and upper limits of the output power of the generator / natural gas source located at node b respectively; Indicates the maximum ramping value of the generator / natural gas source output at node b; Represents the upper limit of renewable energy output of node b at time t.

[0085] (c): Line power constraints:

[0086]

[0087]

[0088] in, is the upper limit of the transmission power of line L.

[0089] (d): Energy hub component model constraints:

[0090] For energy conversion elements, the input and output power satisfy equation (28)

[0091]

[0092] in, and They represent the lower and upper limits of the input power of energy form k of element j respectively.

[0093] The input and output constraints of the energy storage element are as follows:

[0094]

[0095]

[0096] in, and They represent the lower and upper limits of the output power of energy form k of element j respectively.

[0097] The state of charge (SOC) of the energy storage element satisfies the constraints of equations (31)-(32):

[0098]

[0099] S j (0) = S j (T) (32)

[0100] Among them, S j and They represent the lower and upper limits of the SOC of element j respectively; T is the total number of time periods in a day. In order to ensure that the energy storage elements do not charge and discharge at the same time, formula (33) is introduced:

[0101] sta_cha b (t)+sta_dis b (t)≤1 (33)

[0102] (e): Constraints of power grid and gas grid models:

[0103] The DC power flow constraints of the power grid are given by equations (13) and (14); the flow equation of the gas grid and the piecewise linearization model constraints of the flow equation are given by equations (16)-(21).

[0104] (f): objective function:

[0105] For the above-mentioned integrated energy system planning model taking into account demand response and taking into account the Z number, the economic minimum is set as the objective function. The objective function cost is divided into four parts. The first part is the equal annual value F1 of the cost of expanding the power grid and gas grid, as shown in formula (34):

[0106]

[0107] The second part is the annual operating cost of the entire system, which is the output cost of the generator and natural gas source, as shown in formula (35):

[0108]

[0109] The third part is the annual wind curtailment cost of the system, and the specific formula is shown in formula (36):

[0110]

[0111] The fourth part is the demand response cost, and the specific formula is shown in formula (37):

[0112]

[0113] The objective function requires the sum of the four costs to be minimized, that is,

[0114] Obj=min(F1+F2+F3+F4) (38)

[0115] Among them, Obj is the objective function, K1 represents the set of energy forms of energy input to the energy hub; represents the construction cost of the line L to be built, R represents the residual value rate of fixed assets, r represents the discount rate, τ l represents the operating life of the line L to be built; D represents the total number of days in a year, N S represents the total number of scenes, ω(s) represents the probability of scene s, represents the unit operating cost of the generator or natural gas source located at node b, They represent the unit wind curtailment cost of renewable energy at node b, the unit reward cost of customer demand response, and the unit default penalty cost, respectively. b RE and They represent the wind curtailment of renewable energy at node b and the committed load reduction signed by the customer.

[0116] The Gurobi toolbox is used to solve the above mixed integer linear programming and obtain the corresponding planning results.

[0117] The beneficial effects of the present invention are as follows:

[0118] The integrated energy system planning method considering the credibility of demand response proposed by the present invention has the following advantages over the traditional integrated energy system planning method considering demand response: the present invention establishes a model to describe the performance of customer demand response based on the Z number, which solves the problem of incomplete uncertainty modeling in the traditional method. In actual application scenarios, the performance of power user demand response may depend on a variety of uncertain factors. In traditional methods, in order to effectively describe this uncertainty, a thorough demand-side survey is generally required, and then modeling is performed using probability or possibility methods based on the obtained data. Therefore, traditional methods require that the information for uncertainty modeling must be completely accessible and credible. However, in practical applications, it is often difficult to obtain completely reliable information about each user. Therefore, the present invention solves the problem that traditional statistical models can only describe uncertainty characteristics that depend on historical data and cannot consider the impact of future influencing factors by establishing a model to describe the performance of customer demand response based on the Z number. This method can simultaneously consider the impact of uncertainty factors and information credibility on customer demand response performance, thereby making a more accurate and comprehensive estimate of the actual situation. BRIEF DESCRIPTION OF THE DRAWINGS

[0119] Figure 1 is the trapezoidal probability density function of the demand response performance.

[0120] Figure 2 The figure is a flowchart of a solution according to an embodiment of the present invention.

[0121] Figure 3 It is comprehensive energy load and wind power data.

[0122] Figure 4 The weighted membership function of scenarios case 3 and case 4 designed to solve the demand response model. DETAILED DESCRIPTION

[0123] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0124] like Figure 2 As shown, the present invention provides a comprehensive energy system planning method considering the credibility of demand response, and the specific implementation steps are as follows:

[0125] Step 1: Obtain comprehensive energy load data, wind power data and other system input data on the energy hub side such as Figure 3 As shown. Three scenarios are selected: summer, winter and transition season. The cold, hot and electric load data adopt the load demand in the Mediterranean region. The output of renewable energy is also divided into the above three scenarios. The output curve comes from the National Renewable Energy Laboratory of the United States. The electric load of the grid node is 50% of the electric load on the energy hub side. The parameters of the three generators are shown in Table 1:

[0126] Table 1: Generator parameters

[0127]

[0128] The natural gas source parameters are shown in Table 2:

[0129] Table 2: Natural gas source parameters

[0130]

[0131] The parameters of the built power transmission lines and gas pipelines are shown in Tables 3 and 4. The parameters of the lines to be built are the same:

[0132] Table 3: Grid parameters

[0133]

[0134] Table 4: Natural gas grid parameters

[0135]

[0136] Other planning cost parameters are shown in Table 5:

[0137] Table 5: Planning cost parameters

[0138] type Cost (yuan / kw / h) DR Incentive 2 DR Penalty 1.2 Wind curtailment 2.5

[0139] Step 2: Solve the demand response model;

[0140] The data on user demand response performance from 2008 to 2022 published on the PJM official website and ISO New England official website are selected. In order to highlight the effect of comparison, four scenarios are set: case 1: DR is not considered; case 2: DR is considered, but the uncertainty of customer demand response performance is not considered, that is, it is assumed that customers will fulfill their load reduction according to the contract; case 3: The uncertainty of customer demand response performance is described by the Z number. Part of it is represented by a trapezoidal function, Part of it is expressed by trigonometric functions; case 4: using Z number to describe the uncertainty of customer demand response performance, Part of it is represented by a trapezoidal function, Part is 1, that is, For case 3 and case 4, the anti-fuzzy number method and centroid method can be used to convert them into probability density functions, and then Monte Carlo sampling is used for 1000 times to take the average value, and matlab is used for simulation. The weighted membership functions of case 3 and case 4 are as follows Figure 4As shown in the figure, the mean value of the probability density sampling in case 3 is 0.75, and that in case 4 is 0.7. This is because the reliability of the probability density function of the customer demand response performance is not 100%, and the customer demand response performance will be worse in practice. Therefore, using the traditional probability method to describe the customer demand response performance will lead to over-optimism, and the planning results will have the risk of line overload.

[0141] Figure 1 is the trapezoidal probability density function of the demand response performance.

[0142] Step 3: Solve the integrated energy system planning model considering energy hubs;

[0143] Using the demand response results under several scenarios obtained in step 2, the integrated energy system planning model considering the energy hub is solved. The yalmip toolbox is used to model and call the Gurobi solver for solving. The simulation results are shown in Table 6 below:

[0144] Table 6: Planning results

[0145]

[0146] The simulation data shows that compared with case 1, case 2 does not consider the uncertainty of customer demand response. It is believed that customers will fulfill the load reduction amount signed in the contract as agreed. Reducing a small amount of load during peak load can reduce the construction investment of a line. The total cost of case 2 is the lowest. Compared with case 2, case 3 and case 4 consider the uncertainty of customer demand response performance. The difference is that case 3 considers the reliability of the probability density function of customer demand response performance, and the performance result will be slightly worse, and the cost will be higher than case 4. According to the cost analysis of planning results under different circumstances, it can be concluded that the method proposed in the present invention can reduce economic costs and can be applied in the field of integrated energy system planning. Therefore, the above data analysis shows the effectiveness of this method.

Claims

1. A comprehensive energy system planning method considering the credibility of demand response, characterized in that: Here are the steps: Step (1), obtaining comprehensive energy load data, wind power data, generator data, natural gas source data, power grid data, natural gas grid parameter data and planning cost parameter data on the energy hub side; Step (2), establishing a demand response model; Step (3), establishing a comprehensive energy system planning model considering energy hubs; the lower layer of the comprehensive energy system is the power transmission line and the gas transmission line, and the upper layer is the energy hub; Step (3.1): Establish an energy hub model: In the energy hub model, there are three types of component models: Energy conversion components include combined heat and power (CHP), compression chillers (CERG) and electric heat pumps (HP), energy storage components include heat storage (HS) and cold storage (CS) devices, and renewable energy components include wind power (WT); (a): Energy conversion element: The input-output relationship of the energy hub is expressed as: Among them, I b,e (t) and I b,g (t) is the electric power and natural gas flow transmitted from the lower grid to the energy hub connected to node b at time t, P b,wt (t) is the output power of the wind power plant at the energy hub connected to node b at time t, and is the efficiency of natural gas conversion into electricity and heat, η HP and η CE are the efficiency of electric heating by electric heat pump and electric cooling by refrigerant, P b,hs (t) and P b,cs (t) are the charging power of heat storage and cold storage in the energy hub connected to node b at time t; v b,1 (t), v b,2 (t), and (1-v b,1 (t)-v b,2 v′ is the dispatch factor for the electric power transmitted from the lower grid to the energy hub connected to node b at time t to distribute to the electric load, compression refrigerator and electric heat pump; b,1 (t), v′ b,2 (t), and (1-v′ b,1 (t)-v′ b,2 v″ is the dispatch factor for distributing the electric power output of the cogeneration in the energy hub connected to node b at time t to the electric load, compression refrigerator and electric heat pump; b,1 (t), v b,2 (t)″, and (1-v″ b,1 (t)-v″ b,2 (t)) is the dispatch factor for the distribution of wind power output electric power to the electric load, compression refrigerator and electric heat pump in the energy hub connected to node b at time t; (b): Energy storage element: For the energy storage element, the model is as follows: Among them, S b,j (t) represents the charge state of the energy hub system component j connected to node b at time t; and Represent the charging and discharging efficiency of system component j respectively; Δt represents the time span of a single time period; and They represent the input and output power of the system element j connected to the energy hub of node b at time t in the form of energy k; k = c is the energy form of cold, k = h is the energy form of heat, and k = e is the energy form of electricity; for the heat storage element, if P b,hs (t)<0, then otherwise, The same applies to cold storage elements; E j Table system component j storage capacity; sta_cha b (t) and sta_dis b (t) are all 0-1 variables, representing the charging and discharging states of the energy storage element in the energy hub connected to node b at time t; (c): Renewable energy components: Distributed renewable energy components are zero-input single-output components. The output of wind turbines is shown in the formula: Among them, P r is the rated output of the wind turbine, v is the actual wind speed, and v i is the cut-in wind speed, v0 is the cut-out wind speed, v r is the rated wind speed; Step (3.2): Establish the energy network model of the integrated energy system: (a): The DC power flow model used in the power grid is as follows: Where: α represents the energy form transmitted from the lower grid to the energy hub, α = e represents the electric power, and α = g represents the natural gas flow; represents the power on line L at time t; θ L (t) and θ′ L (t) respectively represent the phases of the starting point and the end point of the transmission line L at time t; x L Represents the reactance of line L; LE α and LC α Represent the collection of built and to-be-built lines respectively; is a 0-1 variable indicating whether the line L to be built is constructed or not, M e is a constant of 1000; Equation (14) uses the Big M method to establish the power flow formula for the new line; (b): Gas network model: The Weymouth steady-state model is used to describe the flow of natural gas pipelines, that is, the flow is only related to the pressure at both ends of the pipeline, and flows from the node with high pressure to the node with low pressure; for the natural gas pipeline L, the first and last nodes are n and n′ respectively, and the specific expression of the natural gas flow through the pipeline and the node pressure is: in, Indicates the natural gas flow in the pipeline; c n-n′ is the pipeline transmission parameter, which is related to the pipeline length, diameter and temperature; p n and p n′ denote the gas pressure of nodes n and n′ in the natural gas grid respectively; sgn(p n ,p n′ ) is a direction variable, representing the flow direction of natural gas in the pipeline. n ,p n′ )=1, natural gas flows from node n to node n′, sgn(p n ,p n′ ) = -1, natural gas flows from node n′ to node n; In order to reduce the difficulty of solving the natural gas subsystem model and the complexity of its coupling with the power subsystem, the nonlinear natural gas pipeline flow equation is linearized by piecewise linearization, and then the nonlinear model is transformed into a mixed integer linear programming model; the nonlinear term in the natural gas subsystem model, that is, the square of the node pressure in equation (15), is converted into a mixed integer linear programming model by introducing a new variable p′ and setting p′=p 2 By removing the square term, equation (15) becomes equation (16): The square term on the left side of equation (16) Perform piecewise linearization, the specific formula is equation (17)-(21), and we get: Among them, N represents the total number of linear segments. The larger the value of N, the more linear segments there are and the higher the accuracy. l represents the current number of segments. Represents the slope of a linear piecewise function; represents the intercept of the linear piecewise function on the y-axis, It is a 0-1 variable, indicating whether the point in question falls within the linear segmentation interval. 1 means it falls within the interval, and 0 means it does not fall within the interval. represents the length of the linear piecewise function; Step (4), solve the mixed integer linear programming and obtain the corresponding planning results.

2. A comprehensive energy system planning method considering demand response credibility according to claim 1, characterized in that: The specific method of step (1) is as follows: Obtain the cooling, heating and electricity load data of typical days in summer, transition season and winter on the energy hub side, the wind power output data on three typical days, the output upper limit, ramp rate and unit output cost data of the generator, the output upper limit and unit output cost data of the natural gas source, the power upper limit, reactance and investment cost data of the power grid, the flow upper limit, pipeline coefficient and investment cost data of the natural gas network, DR incentives, DR penalties and unit cost data of wind abandonment.

3. A comprehensive energy system planning method considering demand response credibility according to claim 2, characterized in that: The specific steps of step (2) are as follows: Step (2.1), use Z number to establish customer demand response performance uncertainty model; introduce variable u b The customer's demand response performance is expressed as the ratio of the customer's actual load reduction to the amount contracted by the customer. The mathematical expression is: Among them, u b represents the demand response performance of the customer at node b, represents the actual load reduction of the customer on node b at time t, represents the load reduction contracted by the customer at node b at time t; if u b =0, which means that the customer does not perform load reduction at all; u b The value of is uncertain. In order to take this uncertainty into account, a method based on the Z number is used. According to the Z number method, for any uncertain variable, its impact can be expressed by an ordered pair of fuzzy numbers. To indicate that Is the restriction part, indicating that u b The assumptions of restriction, The reliability part indicates The data credibility is a measure of the customer's demand response performance; therefore, the mathematical formula of the customer's demand response performance Z number can be expressed as: in, Among them, x1-x4 and x′1-x′4 are the trapezoidal membership functions defined and triangle membership function characteristic parameters; In order to convert the model into a mixed integer linear programming for easy solution, the fuzzy expectation method is used to convert the Z number in the model into a conventional fuzzy quantity, and then the centroid method is used to convert the obtained fuzzy variables into corresponding probability variables. The implementation steps are as follows: (a) Convert the reliability part into a clear number, namely the reliability weight: Where β is the reliability weight, Representation function The integral value within its domain is Representation function The value of the integral within its domain; (b) Add the obtained reliability weight to the restricted part, and the weighted Z number can be expressed as: in, is the weighted Z number, It is weighted Membership function; (c) Convert irregular fuzzy numbers into regular fuzzy numbers: in, is the regular fuzzy number of the weighted fuzzy number; (d) Determined by the centroid method The corresponding equivalent probability density function is: Where f(x) is the corresponding probability density function; Step (2.2), simplify the model; use the Monte Carlo method to sample the corresponding probability density function f(x) n times, and finally take the average to obtain u b ; For the trapezoidal probability density function f(x) of the above demand response performance, the parameters are x1, x2, x3, x4 respectively; first discretize it and insert the number x between x1 and x2 1,2 , insert the number x between x2 and x3 2,3 , insert the number x between x3 and x4 3,4 ; Divide x1-x4 into six segments evenly; the interval [x1,x 1,2 ]、[x 1,2 ,x2], [x2,x 2,3 ]、[x 2,3 ,x3]、[x3,x 3,4 ]、[x 3,4 ,x4] The corresponding probabilities P are P1-P6 respectively. Monte Carlo sampling is performed. If the generated random number r is between [0,P1], then the corresponding u b for If the random number r is between [P1, P1+P2], then u b for If the random number r is between [P1+P2,P1+P2+P3], then u b for If the random number r is between [P1+P2+P3,P1+P2+P3+P4], then u b for If the random number r is between [P1+P2+P3+P4,P1+P2+P3+P4+P5], then u b for If the random number r is between [P1+P2+P3+P4+P5,P1+P2+P3+P4+P5+P6], then u b for As shown below: Perform the Monte Carlo sampling process as shown above n times, take the average of the results, and obtain the customer demand response performance u b .

4. A comprehensive energy system planning method considering demand response credibility according to claim 3, characterized in that: The specific method of step (4) is as follows: The planning model constraints of the system are divided into node energy conservation constraints, power / gas source constraints, line power constraints, three types of component model constraints, and power grid and gas grid model constraints, as follows: (a): Node energy conservation constraint: The energy hub can be considered as a load of multiple energy forms on the lower grid node; in, represents the output of the generator or natural gas source at node b at time t, L1 b and L2 b They represent the set of lines starting and ending at node b, respectively. represents the output power of the renewable energy source connected to node b at time t; (b): Power supply / gas source constraints: Formula (23) is the upper and lower limit constraints of generator / natural gas source output; Formula (24) is the generator / natural gas source output ramp constraint; Formula (25) is the centralized renewable energy output constraint: in, and They represent the lower and upper limits of the output power of the generator / natural gas source located at node b respectively; Indicates the maximum ramping value of the generator / natural gas source output at node b; represents the upper limit of renewable energy output of node b at time t; (c): Line power constraints: in, is the upper limit of the transmission power of line L; (d): Energy hub component model constraints: For energy conversion elements, the input and output power satisfy equation (28) in, and They represent the lower and upper limits of the input power of energy form k of element j respectively; The input and output constraints of the energy storage element are as follows: in, and They represent the lower and upper limits of the output power of energy form k of element j, respectively; The state of charge SOC of the energy storage element satisfies the constraints of equations (31)-(32): S j (0)=S j (T) (32) in, S j and They represent the lower and upper limits of the SOC of element j respectively; T is the total number of time periods in a day; in order to ensure that the energy storage elements do not charge and discharge at the same time, formula (33) is introduced: sta_cha b (t)+sta_dis b (t)≤1 (33) (e): Constraints of power grid and gas grid models: The DC power flow constraints of the power grid are given by equations (13) and (14); the flow equation of the gas grid and the piecewise linearization model constraints of the flow equation are given by equations (16)-(21); (f): objective function: For the integrated energy system planning model considering demand response and taking into account the Z number, the economic minimum is set as the objective function. The objective function cost is divided into four parts. The first part is the equal annual value F1 of the expansion cost of the power grid and gas grid, as shown in formula (34): The second part is the annual operating cost of the entire system, which is the output cost of the generator and natural gas source, as shown in formula (35): The third part is the annual wind curtailment cost of the system, and the specific formula is shown in formula (36): The fourth part is the demand response cost, and the specific formula is shown in formula (37): The objective function requires the sum of the four costs to be minimized, that is, Obj=min(F1+F2+F3+F4) (38) Among them, Obj is the objective function, K1 represents the set of energy forms of energy input to the energy hub; represents the construction cost of the line L to be built, R represents the residual value rate of fixed assets, r represents the discount rate, τ l represents the operating life of the line L to be built; D represents the total number of days in a year, N S represents the total number of scenes, ω(s) represents the probability of scene s, represents the unit operating cost of the generator or natural gas source located at node b, c b cut , c b boc They represent the unit wind curtailment cost of renewable energy at node b, the unit reward cost of customer demand response, and the unit default penalty cost, respectively. and They represent the wind curtailment of renewable energy on node b and the committed load reduction signed by the customer; The Gurobi toolbox is used to solve the above mixed integer linear programming and obtain the corresponding planning results.

Images