Multi-energy coupling aggregator feasible region rapid identification method

Abstract

A multi-energy coupling aggregator feasible region rapid identification method relates to the technical field of comprehensive energy system operation optimization, and comprises the following steps: firstly, establishing an optimization operation model which aims at minimizing the total operation cost and covers electricity-gas-heat multi-energy coupling and equipment operation constraints; secondly, on the basis of the model, defining an interactive energy feasible region (IEFR) and a flexible climbing feasible region (FRFR) to be identified; then, a max-min robust optimization model used for depicting a feasible region boundary is constructed, and the max-min robust optimization model is converted into a mixed integer second-order cone programming MISOCP model capable of being efficiently solved through a strong dual theory and a KKT condition; and finally, carrying out iterative solution by adopting a polyhedral projection algorithm based on a dichotomy and a plane cutting mechanism, and rapidly obtaining an accurate boundary of the IEFR and the FRFR.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of integrated energy system operation optimization, and particularly relates to a multi-energy coupling aggregator feasible region rapid identification method. BACKGROUND

[0002] Under the "double carbon" target, the multi-energy aggregator with deep coupling of electricity-gas-heat has become a new subject of providing auxiliary services such as flexible ramping (FRP) to the power grid, due to its multi-time scale regulation capability and cross-energy complementary advantage. However, the internal multi-energy system structure is complex, and the traditional evaluation method has obvious limitations: the iterative solution method has the risk of data privacy leakage, and the calculation efficiency is low and the convergence speed is slow, which is difficult to meet the real-time scheduling requirements; and the traditional feasible region equivalent method has problems such as complex modeling and difficult calculation when dealing with high-dimensional multi-energy systems. Therefore, the existing technology cannot realize the rapid and accurate evaluation of the flexible resources of such aggregators, hindering their efficient participation in the power grid interaction and auxiliary service market. SUMMARY

[0003] The present application discloses a multi-energy coupling aggregator feasible region rapid identification method: under the electricity-gas coordination framework, an aggregator operation model is established to minimize the distance of the relaxation variable to the feasible region boundary, and a max-min robust model considering the flexible ramping constraint is constructed to locate the feasible region boundary; the nonlinear constraint of the natural gas network is processed by second-order cone relaxation, and then converted into a second-order cone programming containing complementary constraints by means of strong duality and KKT conditions, and then gradually approaching IEFR / FRFR by using the polyhedral projection algorithm of bisection-cutting plane, so as to avoid the dimension disaster of traditional enumeration in high-dimensional scenarios.

[0004] To achieve the above-mentioned purpose, a multi-energy coupling aggregator feasible region rapid identification method is disclosed, comprising the following steps:

[0005] Establishing a multi-energy coupling aggregator optimal operation model with the total operation cost minimization as the objective function;

[0006] Based on the optimal operation model, defining the interactive energy feasible region IEFR and the flexible ramping feasible region FRFR as the feasible region to be identified;

[0007] Constructing a max-min robust optimization model for identifying the boundary of the feasible region;

[0008] Converting the max-min model into a mixed integer second-order cone programming (MISOCP) model based on the mathematical transformation of the dual theory and KKT conditions;

[0009] The polyhedral projection algorithm based on dichotomy and cutting plane mechanism is used to solve the mixed integer second order cone programming MISOCP model, and the boundary of the interactive energy feasible region IEFR and the flexible ramping feasible region FRFR is iteratively calculated and output.

[0010] As a further scheme of the present application, a multi-energy coupling aggregator optimal operation model is established, with the total operation cost as the objective function, which is specifically expressed as:

[0011]

[0012] In the formula, the first term is the relevant cost of energy procurement, the second term is the operation and maintenance cost of four types of energy conversion devices, i.e., combined heat and power unit CHP, gas turbine GT, electricity-to-gas device P2G and electric heating boiler EB, the third term is the upward flexibility reserve cost provided by the four types of energy conversion devices, i.e., CHP, GT, P2G and EB, the fourth term is the downward flexibility reserve cost provided by the four types of energy conversion devices, i.e., CHP, GT, P2G and EB, the fifth term is the income of providing flexible ramping service, and the last term is the penalty cost of curtailment of wind and light, and are respectively the time-of-use electricity price and gas price of the upper grid and the natural gas network; and are respectively the electricity and natural gas purchased by the multi-energy coupling aggregator from the upper energy network, i.e., the energy interaction amount of the multi-energy coupling aggregator and the upper energy network; CHP , r GT , r P2G and r EB are respectively the unit output operation and maintenance cost coefficients of the four types of devices, i.e., CHP, GT, P2G and EB; flexi_RU is the upward flexibility reserve cost coefficient of the four types of devices, i.e., CHP, GT, P2G and EB; flexi_RD is the downward flexibility reserve cost coefficient of the four types of devices, i.e., CHP, GT, P2G and EB; and are respectively the output powers of the four types of devices, i.e., CHP, GT, P2G and EB; and represent the upward and downward flexibility reserve of the four types of devices, i.e., CHP, GT, P2G and EB; is the unit price of providing flexible ramping service for the upper grid; FRU t and FRDt To provide upward and downward flexible ramping services for the upper power grid for multi-coupling aggregators; cur To reduce the penalty cost of new energy; To reduce the amount of new energy.

[0013] As a further scheme of the present application, the constraint conditions of the optimization operation model specifically include:

[0014] The device operation constraints of four types of energy conversion devices, including combined heat and power units CHP, gas turbines GT, electricity-to-gas devices P2G, and electric heating boilers EB, including energy conversion constraints, output limit constraints, and device ramping constraints;

[0015] New energy output constraints, including distributed wind turbines WT and distributed photovoltaic PV;

[0016] Linearized power flow constraints of the distribution network;

[0017] Second-order cone relaxation constraints of the natural gas network;

[0018] Heat network balance constraints.

[0019] As a further scheme of the present application, the combined heat and power unit CHP generates electricity and produces heat by consuming natural gas, and the relationship between its electricity generation and natural gas consumption is as follows:

[0020]

[0021] The relationship between electricity generation and heat production is as follows:

[0022]

[0023] The output and ramping constraints are as follows:

[0024]

[0025] In the formula, and are the heat output and natural gas consumption of the combined heat and power unit CHP, respectively; and are the gas-electricity and electricity-heat conversion efficiencies of the combined heat and power unit CHP, respectively; and are the upper and lower limits of the output of the combined heat and power unit CHP; and are the upper and lower ramping limits of the combined heat and power unit CHP, respectively;

[0026] The gas turbine GT produces electricity by burning natural gas, and its energy conversion relationship is shown in the following formula:

[0027]

[0028] The operation constraints are as follows:

[0029]

[0030] wherein, is the gas consumption of the gas turbine GT; is the gas-electricity conversion efficiency of the gas turbine GT; and are the upper and lower output limits of the gas turbine GT, respectively; and are the upper and lower ramping limits of the gas turbine GT, respectively;

[0031] The electricity-to-gas device P2G converts electricity into natural gas, and the relationship between the electricity consumption and the natural gas output of the electricity-to-gas device P2G is as follows:

[0032]

[0033] The operation constraints are as follows:

[0034]

[0035] wherein, is the natural gas output of the electricity-to-gas device P2G; and are the upper and lower output limits of the electricity-to-gas device P2G, respectively; and are the upper and lower ramping limits of the electricity-to-gas device P2G, respectively;

[0036] The electric heating boiler EB meets the demand of the heat load by consuming electricity, and the relationship between the electricity consumption and the heat output of the electric heating boiler EB is as follows:

[0037]

[0038] The constraints are as follows:

[0039]

[0040]

[0041] wherein, is the heat output of the electric heating boiler EB; and are the minimum power and the maximum power of the electric heating boiler EB, respectively; and are the upper and lower ramping limits of the electric heating boiler EB, respectively;

[0042] The operation constraints of the distributed wind power WT and the distributed photovoltaic PV are as follows:

[0043]

[0044] In the formula: and are the maximum output power of distributed wind power WT and distributed photovoltaic PV respectively; and represent the reduction of distributed wind power WT and distributed photovoltaic PV.

[0045] As a further scheme of the application, the power distribution network linearization flow constraint is as follows:

[0046]

[0047] In the formula, and are the minimum and maximum power purchase amount of the multi-energy coupling aggregator; P ij,t is the active power flow on the line; P ij,max is the maximum active power flow of the line; is the electrical load, δ elec is a preset safety margin constant for sufficient operation buffer;

[0048] The energy balance formula is rewritten using the chance constraint:

[0049]

[0050] Where, pr{} represents the probability of the event in {} occurring, ζ is the confidence level, and γ is a random variable subject to a logistic distribution with a location parameter of and a scale parameter s;

[0051] The chance constraint is converted into a deterministic constraint using the deterministic equivalence method, as follows:

[0052]

[0053] K ζ =F -1 (1-ζ)

[0054] Where, F -1 is the inverse function of the cumulative distribution function of the logistic distribution.

[0055] As a further scheme of the application, the natural gas network second-order cone relaxation constraint: the power distribution network and the natural gas network are connected through multi-energy coupling devices, the influence of the compressor is not considered, the steady-state natural gas network model includes node flow balance constraints and steady-state transmission pipeline operation constraints, and the node flow balance constraint is as follows:

[0056]

[0057] The Weymouth equation is as shown in the following formula:

[0058]

[0059] The pipeline flow constraint is as shown in the following formula:

[0060] -G mn,max ≤G mn ≤G mn,max

[0061] The node pressure constraint is as shown in the following formula:

[0062] π m,min +δ gas ≤π m ≤π m,max -δ gas π m,min ≤π m ≤π m,max

[0063] The constraint of the multi-energy coupling aggregator purchasing natural gas is as shown in the following formula:

[0064]

[0065] In the formula: is the natural gas interaction amount of the multi-energy coupling aggregator and the upper natural gas network; and are the minimum and maximum amounts of natural gas purchased by the multi-energy coupling aggregator respectively; G mn is the tidal flow of the pipeline mn; is the natural gas load of the node m; π m is the gas pressure of the node m; m mn is the Weymouth characteristic parameter; G mn,max is the maximum tidal flow of the pipeline mn; π m,max and π m,min are the minimum and maximum gas pressures of the node m respectively, δ gas is a preset safety margin constant for leaving sufficient operation buffer;

[0066] The Weymouth equation is converted into a second-order cone form, and the problem is rewritten as a second-order cone programming problem:

[0067]

[0068] A new variable is introduced:

[0069]

[0070] The Weymouth equation is converted into the following formula:

[0071] (G mn ) 2 ≤τ mn

[0072] The transformation into the second-order cone constraint form is as follows:

[0073]

[0074] The heat network balance constraint is:

[0075]

[0076] In the formula, is the heat load in the multi-energy coupling aggregator, and ηtheat} is the heat network delivery efficiency coefficient, which is a constant close to 1.

[0077] As a further scheme of the present application, the multi-energy coupling aggregator constraint condition is expressed in a compact form as follows:

[0078]

[0079] In the formula, x is a variable related to energy interaction and FRP, that is, FRU t and FRD t ; y1 is a controllable variable of linear constraints, y2 is a controllable variable of nonlinear constraints, and the feasible region range of x is a convex hull formed by the projection of all constraints in the model on the x-axis. The mathematical expression of the feasible region is:

[0080]

[0081] For any x, there is a set of scheduling schemes y1 and y2 corresponding thereto. The definitions of the scheduling variables y1 and y2 for a fixed x' are as follows:

[0082]

[0083] As a further scheme of the present application, the max-min robust optimization model for identifying the boundary of the feasible region is specifically constructed as follows:

[0084] The following second-order cone programming problem is constructed by introducing a slack variable:

[0085]

[0086] s.t.Ax+By1-Iv≤b

[0087]

[0088] In the formula: v is a slack variable, and The following transformations are made:

[0089] (2y2) 2 +[(cy2+Iv)-1] 2 ≤[(cy2+Iv)+1] 2

[0090]

[0091] The compact form is:

[0092] ||Ey2+Iv+q i ||≤Fy2+Iv+d i

[0093] The objective function of the second-order cone programming problem is to minimize the adjustment that satisfies all constraints, f(x) = 0 represents that x is a point in the feasible region and Otherwise, if f(x) > 0, x is a point outside the feasible region and For a fixed x', Y FR is non-empty, the optimal value of the second-order cone programming problem is 0, and the model is rewritten as:

[0094]

[0095] Ax+By1-Iv≤b

[0096] ||Ey2+Iv+q i ||≤Fy2+Iv+d i

[0097] The concept of robust optimization is introduced, and the model is re-expressed as:

[0098]

[0099] s.t.Ax+By1-Iv≤b:π

[0100] ||Ey2+Iv+q i ||≤Fy2+Iv+d i :(ω i ,λ i ).

[0101] As a further scheme of the application, the method for converting the max-min model into a mixed integer second-order cone programming (MISOCP) model through mathematical transformation is as follows:

[0102] Based on the duality theory, the min problem inside the model is converted into a max problem for solving, and the dual problem of the min problem inside the model is:

[0103]

[0104]

[0105] ||ω i ||≤λ i

[0106] -1≤π≤0

[0107] There is a bilinear term in the objective function of the above formula, and the constraints of x and π are independent of each other. Based on robust optimization, the following form is obtained:

[0108]

[0109] For the max problem inside the model:

[0110] max(-Axπ)

[0111] s.t.Dx≤d:u

[0112] According to the strong duality theory, we get:

[0113] -Axπ=d·u

[0114] The KKT conditions for the max problem inside the model are:

[0115] A T π+D T u=0

[0116] 0≤u T ⊥(d-Dx)≤0

[0117] The constraint 0≤u T ⊥(d-Dx)≤0 is equivalent to the three constraints 0≤u T , (d-Dx)≤0 and u T (d-Dx)=0; obtain the mixed integer second order cone programming MISOCP model as follows:

[0118]

[0119] ||ω i ||≤λ i

[0120] -1≤π≤0

[0121] A T π+D T u=0

[0122] 0≤(d-Dx)≤M(δ-1)

[0123] 0≤u T ≤Mδ.

[0124] As a further scheme of the present application, the method for solving the MISOCP model by using a polyhedral projection algorithm based on dichotomy and cutting plane mechanism, iteratively calculating and outputting the boundary of the interactive energy feasible region IEFR and flexible ramping feasible region FRFR is as follows:

[0125] Step 1: for variable x, initialize a large enough region and select a point in the feasible region, denoted as x in ;

[0126] Step 2: by solving the mixed integer second-order cone programming MISOCP model, obtain the optimal solution and optimal value f1, and in the worst case, select an external point x out ; if f1(x out )>0, it means that x out is outside the feasible region, record x out and enter step 3; otherwise, f1(x out )=0, it means that x out is inside the feasible region, exit and record

[0127] Step 3: when x in is an internal point of the feasible region and x out is an external point of the feasible region, the connecting line between the internal point and the external point must intersect with the boundary of the feasible region, select dichotomy to find the feasible region boundary point x o , introduce auxiliary parameter γ, determine x o by γ=0.5*(γ1+γ2) and x in =(1-γ)·x out +γ·x o , solve whether x o is a feasible region boundary point, when x o is determined, the corresponding π o , λ oi and ω oi in the feasible region are also determined, record π o , λ oi and ω oi ;

[0128]

[0129] ||ω′ i ||≤λ′ i

[0130] -1≤h≤0

[0131] If f2(x o )≠0, it means that x o is not a boundary point of the feasible region, update γ and xo ; if f2(x o ) = 0, record the optimal solution of the above formula, i.e. π o , λ oi and ω oi , and go to step 4;

[0132] Step 4: after x o is determined, add a boundary constraint to the mixed integer second-order cone programming (MISOCP) model as a cutting plane using the known π o , λ oi and ω oi ; update the initial region return to step 2;

[0133] The iteration process of steps 2 to 4 continues until one of the following two conditions is met:

[0134] precision convergence: the polyhedral region defined by the cutting planes generated by two consecutive iterations has a volume change or a maximum displacement of boundary points less than a preset convergence tolerance;

[0135] iteration number convergence: the number of iterations reaches a preset maximum number of iterations.

[0136] The beneficial effects of the present application are:

[0137] 1. greatly improve the calculation efficiency and speed: by second-order cone relaxation (SOCP) of the complex nonlinear natural gas network constraint (Weymouth equation), the original problem is converted into a convex optimization form which is easier to solve.

[0138] Using strong duality theory and KKT conditions, the max-min bi-level robust optimization model which is difficult to solve directly is converted into a mixed integer second-order cone programming (MISOCP) model, which can be efficiently processed by modern commercial solvers.

[0139] The polyhedral projection algorithm based on bisection and cutting plane is adopted to gradually approach the boundary of the feasible region, avoiding the "dimension disaster" problem faced by traditional enumeration method in high-dimensional space. This method can quickly eliminate a large number of infeasible regions, has fast convergence speed, and is significantly superior to traditional iterative solution method, which can meet the requirements of real-time scheduling of power systems.

[0140] 2. accurate feasible region description: the max-min robust optimization model constructed in the present application can locate the boundary of the feasible region under the worst-case scenario, ensuring the conservativeness and reliability of the identified feasible region. The solution process of this method can finally output the boundary of the interactive energy feasible region (IEFR) and flexible ramping feasible region (FRFR) of the multi-energy coupling aggregator with high precision, providing clear and accurate flexibility resource evaluation basis for the dispatch center.

[0141] 3. Effective protection of data privacy:

[0142] By mathematical feasible region projection and boundary description, the dispatch center is interacted without disclosing the detailed device parameters, topology structure and real-time operation data inside the aggregator, which overcomes the privacy leakage risk caused by multiple data exchanges in traditional iterative solution.

[0143] 4. Enhance market participation and system flexibility:

[0144] The present application provides key technical support for multi-energy coupling aggregators participating in the electricity market, the reserve market and the flexible ramping (FRP) auxiliary service market. The dispatch center can quickly evaluate the potential of the aggregator to provide peak shaving, frequency modulation and other services based on the feasible region reported by the aggregator, thereby efficiently mobilizing distributed flexible resources and improving the operation reliability and new energy consumption capacity of the entire power system.

[0145] 5. Strong engineering applicability and robustness:

[0146] The present application fully considers the coupling characteristics of electricity-gas-heat multi-energy flow and device operation constraints, and the model construction is comprehensive. The robust optimization framework adopted can effectively cope with the uncertainty of renewable energy output and load demand, making the resulting feasible region result more robust and practical. BRIEF DESCRIPTION OF DRAWINGS

[0147] Figure 1 Flowchart of the polyhedral projection algorithm based on second-order cone programming.

[0148] Figure 2 Schematic diagram of finding boundary points by bisection method.

[0149] Figure 3 Schematic diagram of cutting plane action. DETAILED DESCRIPTION

[0150] The "multi-energy coupling device" described in the present application refers to the main unit inside the multi-energy coupling aggregator for realizing electricity-gas-heat energy conversion and regulation, which specifically includes the following six types:

[0151] Combined heat and power (CHP) unit, which outputs electric energy and heat energy by burning natural gas;

[0152] (1) Gas turbine (GT), which outputs electric energy by burning natural gas;

[0153] (2) Power-to-gas (P2G) device, which converts electric energy into natural gas;

[0154] (3) Electric boiler (EB), which converts electric energy into heat energy;

[0155] (4) Distributed wind power WT, output electric energy;

[0156] (5) Distributed photovoltaic PV, output electric energy.

[0157] As a preferred embodiment of the present application, a multi-energy coupling aggregator feasible region fast identification method has the following specific steps:

[0158] Step 1): A mathematical model with the minimum total operating cost as the objective function is established, which is expressed as:

[0159]

[0160] In the formula, the first term is the relevant cost of energy procurement, the second term is the operation and maintenance cost of the CHP, GT, P2G and EB four types of energy conversion equipment, the third term is the upward flexibility reserve cost provided by the CHP, GT, P2G and EB four types of energy conversion equipment, the fourth term is the downward flexibility reserve cost provided by the CHP, GT, P2G and EB four types of energy conversion equipment, the fifth term is the income of providing flexible ramp service, and the last term is the penalty cost of wind and light abandonment. and are the time-of-use electricity price and gas price of the upper grid and the natural gas network respectively; and are the electricity and natural gas purchased by the multi-energy coupling aggregator from the upper energy network, i.e., the energy interaction amount of the multi-energy coupling aggregator and the upper energy network; CHP , r GT , r P2G and r EB are the unit output operation and maintenance cost coefficients of the CHP, GT, P2G and EB four types of equipment; flexi_RU is the upward flexibility reserve cost coefficient of the CHP, GT, P2G and EB four types of equipment; flexi_RD is the downward flexibility reserve cost coefficient of the CHP, GT, P2G and EB four types of equipment; and are the output powers of the CHP, GT, P2G and EB four types of equipment; and represent the upward and downward flexibility reserve of the CHP, GT, P2G and EB four types of equipment; is the unit price of providing flexible ramp service for the upper grid; FRU t and FRD t are the upward and downward flexible ramp services provided by the multi-energy coupling aggregator for the upper grid; cur is the penalty cost of new energy reduction; is the new energy reduction amount.

[0161] The constraint conditions in step 2) include:

[0162] The operation constraints of the CHP, GT, P2G and EB energy conversion devices include energy conversion constraints, output limit constraints and device ramping constraints; the new energy output constraints include distributed wind power WT and distributed photovoltaic PV; the power distribution network linearization flow constraints; the natural gas network second-order cone relaxation constraints; and the heat network balance constraints.

[0163] The combined heat and power unit CHP can generate electricity and produce heat by consuming natural gas. The relationship between the electricity generation and the natural gas consumption is as follows:

[0164]

[0165] The relationship between the electricity generation and the heat production is as follows:

[0166]

[0167] The output and ramping constraints are as follows:

[0168]

[0169] In the formula, and are the heat output and the natural gas consumption of the combined heat and power unit CHP, respectively; and are the gas-electricity and electricity-heat conversion efficiencies of the combined heat and power unit CHP, respectively; and are the upper and lower limits of the output of the combined heat and power unit CHP; and are the upper and lower ramping limits of the combined heat and power unit CHP.

[0170] The gas turbine GT produces electricity by burning natural gas. The energy conversion relationship is as follows:

[0171]

[0172] The operation constraints are as follows:

[0173]

[0174] In the formula, is the gas consumption of the gas turbine GT; is the gas-electricity conversion efficiency of the gas turbine GT; and are the upper and lower limits of the output of the gas turbine GT; and are the upper and lower ramping limits of the gas turbine GT.

[0175] The P2G device P2G can convert electrical energy into natural gas, and the relationship between the P2G device P2G gas production and electricity consumption is shown in the following formula:

[0176]

[0177] The operation constraints are as follows:

[0178]

[0179] In the formula: is the natural gas output of the P2G device P2G; and are the upper and lower limits of the output of the P2G device P2G; and are the upper and lower limits of the ramp of the P2G device P2G.

[0180] The electric heating boiler EB meets the demand of heat load by consuming electrical energy. The relationship between the power consumption and heat output of the electric heating boiler EB is shown in the following formula:

[0181]

[0182] The constraints are as follows:

[0183]

[0184] In the formula: is the heat output of the electric heating boiler EB; and are the minimum power and maximum power of the electric heating boiler EB; and are the upper and lower limits of the ramp of the electric heating boiler EB.

[0185] The operation constraints of the distributed wind power WT and the distributed photovoltaic PV are as follows:

[0186]

[0187] In the formula: and are the maximum output power of the distributed wind power WT and the distributed photovoltaic PV, respectively; and represent the curtailment of the distributed wind power WT and the distributed photovoltaic PV.

[0188] Linearized distribution network model: The distribution network power flow model is simplified, the network line loss is ignored, and only the active power balance is considered, i.e. the total power on the generation side is balanced with the total power on the load side (including network power flow), as shown in the following formula:

[0189]

[0190] wherein, and Pmin and Pmax are the minimum and maximum power purchase amounts of the multi-energy coupling aggregator; P ij,t P is the active power flow on the line; P ij,max Pmax is the maximum active power flow of the line; P is the electrical load, δ elec δ is the preset safety margin constant for leaving sufficient operation buffer.

[0191] The energy balance equation is rewritten using the chance constraint:

[0192]

[0193] wherein, pr{} represents the probability of the event in {} occurring, ζ is the confidence level, and γ is a random variable subject to a logistic distribution with a location parameter of and a scale parameter s.

[0194] The chance constraint is converted into a deterministic constraint using the certainty equivalent method, as shown in the following equation:

[0195]

[0196] K ζ = F -1 (1-ζ)

[0197] wherein, F -1 is the inverse function of the cumulative distribution function of the logistic distribution.

[0198] Natural gas network model: The power distribution network and the natural gas network are connected through multi-energy coupling devices such as combined heat and power units CHP, gas turbines GT, and P2G devices. The devices connect the power distribution network and the natural gas network, and the impact of the compressor is not considered. The steady-state natural gas network model includes node flow balance constraints and steady-state pipeline operation constraints. The node flow balance constraint is shown in the following equation:

[0199]

[0200] The Weymouth equation is shown in the following equation:

[0201]

[0202] The pipeline flow constraint is shown in the following equation:

[0203] -G mn,max ≤ G mn ≤ G mn,max

[0204] The node pressure constraint is shown in the following equation:

[0205] πm,min + δ gas ≤ π m ≤ π m,max - δ gas π m,min ≤ π m ≤ π m,max

[0206] The constraints for the multi-energy coupled aggregator to purchase natural gas are shown in the following equation:

[0207]

[0208] where: is the natural gas interaction between the multi-energy coupled aggregator and the upper-level natural gas network; and are the minimum and maximum amounts of natural gas that the multi-energy coupled aggregator can purchase, respectively; mn is the tidal flow of the pipeline mn; is the natural gas load at node m; m is the gas pressure at node m; mn is the Weymouth characteristic parameter; mn,max is the maximum tidal flow of the pipeline mn; m,max and m,min are the minimum and maximum gas pressures at node m, respectively; gas is the preset safety margin constant, which is used to leave enough operating buffer.

[0209] The Weymouth equation is converted into a second-order cone form, and the problem is rewritten as a second-order cone programming problem:

[0210]

[0211] Some new variables are introduced:

[0212]

[0213] The Weymouth equation is converted into the following equation:

[0214] (G mn ) 2 ≤ τ mn

[0215] It is converted into a second-order cone constraint form:

[0216]

[0217] The heat network balance constraint: the multi-energy coupled aggregator heat load also depends on the electric heat boiler EB, which can convert electrical energy into thermal energy:

[0218]

[0219] In the formula, The heat load within the multi-energy coupled polymer is η{heat}, which is the heat network transport efficiency coefficient and is a constant close to 1.

[0220] Based on the optimized operation model, this embodiment defines the Interactive Energy Feasible Domain (IEFR) and the Flexible Climbing Feasible Domain (FRFR) as the feasible domains to be identified.

[0221] The IEF (Integrated Energy Free Rate) of a multi-energy coupling aggregator refers to the extent to which the aggregator purchases electricity and natural gas from the upstream power grid and natural gas grid. The FRFR (Fuel Free Rate) of a multi-energy coupling aggregator refers to the extent to which the aggregator can supply all its FRPs (Fuel-Retained Polymers).

[0222] The multi-energy coupled aggregator constraint is expressed in compact form as follows:

[0223]

[0224] In the formula, x is a variable related to energy interaction and FRP, i.e. FRU t and FRD t Let y1 be a controllable variable under linear constraints and y2 be a controllable variable under nonlinear constraints. Therefore, theoretically, the feasible region of x is a convex hull formed by the projections of all constraints in the model onto the x-axis. The mathematical expression for the feasible region is:

[0225]

[0226] For any x, there exists a corresponding set of scheduling schemes y1 and y2. The scheduling variables y1 and y2 for a fixed x′ are defined as follows:

[0227]

[0228] This embodiment constructs a max-min robust optimization model for identifying the feasible region boundary;

[0229] Finding the feasible region of the multi-energy coupled aggregator quotient essentially involves calculating the projections of x, y1, and y2 onto the variable x. Therefore, it is necessary to determine the feasible region for all x. FR x, Y FR Is Y non-empty? To verify Y... FR To assess the feasibility of fixing x′, we introduced slack variables and constructed the following second-order cone programming problem:

[0230]

[0231] stAx+By1-Iv≤b

[0232]

[0233] where v is a slack variable that can be forced to adjust to satisfy all operational constraints. Since is still nonlinear, the following transformation is made:

[0234] (2y2) 2 +[(cy2+Iv)-1] 2 ≤[(cy2+Iv)+1] 2

[0235]

[0236] Its compact form is:

[0237] ||Ey2+Iv+q i ||≤Fy2+Iv+d i

[0238] The objective function of the second-order cone programming problem is to minimize the adjustment to satisfy all constraints, so f(x) = 0 indicates that x is a point in the feasible region and Otherwise, if f(x) > 0, x is a point outside the feasible region and For a fixed x', Y FR is non-empty, and the optimal value of the second-order cone programming problem is 0, and the model is rewritten as:

[0239]

[0240] Ax+By1-Iv≤b

[0241] ||Ey2+Iv+q i ||≤Fy2+Iv+d i

[0242] Considering the possibility that Y FR is infeasible for x ∈ X FR in the worst case, the concept of robust optimization is introduced in this model, and the model is re-expressed as:

[0243]

[0244] s.t.Ax+By1-Iv≤b:π

[0245] ||Ey2+Iv+q i ||≤Fy2+Iv+d i :(ω i ,λ i )

[0246] The max-min model is converted into a mixed integer second-order conic programming (MISOCP) model by mathematical transformation in this embodiment.

[0247] Considering that commercial solvers cannot directly solve the max-min problem, the inner min problem is converted into a max problem based on duality theory. The dual problem of the inner min problem in the model is:

[0248]

[0249] ||ω i ||≤λ i

[0250] -1≤π≤0

[0251] There is a bilinear term in the objective function of the above formula, which leads to difficulty in calculation. The constraints of x and π are independent of each other, so based on robust optimization, the following form can be obtained:

[0252]

[0253] For the inner max problem:

[0254] max(-Axπ)

[0255] s.t.Dx≤d:u

[0256] According to strong duality theory, the following can be obtained:

[0257] -Axπ=d·u

[0258] The KKT conditions for the inner max problem are:

[0259] A T π+D T u=0

[0260] 0≤u T ⊥(d-Dx)≤0

[0261] The constraint 0≤u T ⊥(d-Dx)≤0 is equivalent to the three constraints 0≤u T , (d-Dx)≤0 and u T (d-Dx)=0. Therefore, a mixed integer second-order conic programming (MISOCP) model can be obtained as follows:

[0262]

[0263] ||ω i ||≤λi

[0264] -1≤π≤0

[0265] A T π+D T u=0

[0266] 0≤(d-Dx)≤M(δ-1)

[0267] 0≤u T ≤Mδ

[0268] This embodiment uses a polyhedral projection algorithm based on the bisection method and the cutting plane mechanism to solve the MISOCP model, iteratively calculates and outputs the boundaries of the IEFR and FRFR.

[0269] The basic idea of ​​the polyhedral projection algorithm based on the bisection method and the plane cutting mechanism is to find the boundary of the feasible region in the initial range, and then cut off the infeasible region until there are no infeasible solutions in the space.

[0270] Step 1 (Initialization and Input Data): For variable x, initialize a sufficiently large region. And select a point within the feasible region, denoted as x. in .

[0271] Step 2 (Determine if f is positive or negative): By solving the mixed-integer second-order cone programming (MISOCP) model, obtain its optimal solution and optimal value f1. In the worst case, select an external point x. out If f1(x) out If x > 0, it means x out Outside the feasible region, record x out And proceed to step 3; otherwise, f1(x) out ) = 0 means x out Within the feasible region, exit and log.

[0272] Step 3 (Finding Boundary Points): When x in Let x be an interior point of the feasible region. out When the points are outside the feasible region, the lines connecting them must intersect the boundary of the feasible region. Therefore, the bisection method is used to find the boundary point x of the feasible region. o An auxiliary parameter γ is introduced, which is obtained by using γ = 0.5*(γ1+γ2) and x o =(1-γ)·x in +γ·x out Determine x o Determine x by solving the problem. o Is it a boundary point of the feasible region? When x o When determined, the corresponding π within the feasible regiono , λ oi and ω oi may also be determined, record π o , λ oi and ω oi .

[0273]

[0274] ||ω′ i ||≤λ′ i

[0275] -1≤h≤0

[0276] If f2(x o )≠0, it means x o is not a boundary point of the feasible region, update γ and x o ; if f2(x o )=0, record the optimal solution of the above equation, i.e. π o , λ oi and ω oi , and go to step 4.

[0277] Step 4 (generate a cutting plane): after x o is determined, use the known π o , λ oi and ω oi to add a boundary constraint to the MISOCP model as a cutting plane update the initial region go to step 2.

[0278] The above iterative process (step 2 to step 4) will continue until one of the following two conditions is met:

[0279] precision convergence: the volume change or the maximum displacement of the boundary points of the polyhedral region defined by the cutting planes generated by two consecutive iterations is less than the preset convergence tolerance.

[0280] iteration number convergence: the number of iterations reaches the preset maximum number of iterations to ensure computational efficiency.

[0281] The above specific implementation can be adjusted in different ways by those skilled in the art without departing from the principles and purposes of the present application, the protection scope of the present application is subject to the claims and is not limited by the above specific implementation, each implementation within the scope is subject to the constraints of the present application.

Claims

1. A method for rapid identification of feasible regions of multi-energy coupled aggregation quotients, characterized in that, Includes the following steps: Establish a multi-energy coupled aggregator optimization operation model with the objective function of minimizing total operating cost; Based on the optimized operation model, the interactive energy feasible region (IEFR) and the flexible climbing feasible region (FRFR) are defined as feasible regions to be identified. Construct a max-min robust optimization model for identifying the boundaries of the feasible region; The max-min model is transformed into a mixed-integer second-order cone programming (MISOCP) model based on duality theory and KKT conditions. The mixed integer second-order cone programming (MISOCP) model is solved using a polyhedral projection algorithm based on the bisection method and the cutting plane mechanism. The boundaries of the interactive energy feasible region (IEFR) and the flexible climbing feasible region (FRFR) are calculated iteratively and output.

2. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 1, characterized in that, A multi-energy coupled aggregator quotient optimization operation model is established with the objective function of minimizing the total operating cost, specifically expressed as follows: In the formula, the first item is the cost related to energy procurement; the second item is the operation and maintenance cost of the four types of energy conversion equipment: CHP (Combined Heat and Power) units, GT (Gas Turbine), P2G (Electric Power to Gas), and EB (Electric Boiler); the third item is the upward flexibility reserve cost provided by these four types of energy conversion equipment; the fourth item is the downward flexibility reserve cost provided by these four types of energy conversion equipment; the fifth item is the revenue from providing flexible ramp-up services; and the last item is the penalty cost for wind and solar curtailment. and These are the time-of-use electricity and gas prices for the upstream power grid and the natural gas grid, respectively; and These represent the electricity and natural gas purchased by the multi-energy coupling aggregator from the upstream energy grid, respectively; that is, the energy interaction between the multi-energy coupling aggregator and the upstream energy grid. CHP r GT r P2G and r EB These are the unit output operation and maintenance cost coefficients for four types of equipment: combined heat and power (CHP) units, gas turbines (GT), power-to-gas (P2G) equipment, and electric boilers (EB); r flexi_RU For the four types of equipment: combined heat and power (CHP) units, gas turbines (GT), power-to-gas (P2G) equipment, and electric boilers (EB), the upward flexibility reserve cost coefficient is calculated. flexi_RD Downward flexibility reserve cost coefficient for four types of equipment: combined heat and power unit (CHP), gas turbine (GT), power-to-gas conversion equipment (P2G), and electric boiler (EB). and These are the output powers of four types of devices: CHP, GT, P2G, and EB. and This indicates the upward and downward flexibility of four types of equipment: combined heat and power unit (CHP), gas turbine (GT), electric-to-gas conversion equipment (P2G), and electric boiler (EB). This is the unit price for providing flexible ramping services to the upstream power grid; FRU t and FRD t This refers to the flexible uphill and downhill ramping services provided by multi-energy coupling aggregators to the upstream power grid; cur The penalty costs for reducing new energy sources; To reduce the amount of renewable energy.

3. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 1, characterized in that, The constraints of the optimized operating model specifically include: The equipment operation constraints of four types of energy conversion equipment, namely, combined heat and power (CHP) units, gas turbines (GT), power-to-gas (P2G) equipment, and electric boilers (EB), include energy conversion constraints, output limitation constraints, and equipment ramp-up constraints. Constraints on new energy output, including distributed wind power (WT) and distributed photovoltaic (PV); Power flow constraints for distribution network linearization; Second-order cone relaxation constraint for natural gas network; Heating network balance constraints.

4. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 3, characterized in that, The combined heat and power (CHP) unit generates electricity and heat by consuming natural gas. The relationship between its power generation and natural gas consumption is as follows: The relationship between power generation and heat production is as follows: The output and climbing constraints are as follows: In the formula, and These are the thermal output and natural gas consumption of the combined heat and power (CHP) unit, respectively. and These are the gas-to-electricity and electricity-to-heat conversion efficiencies of the combined heat and power (CHP) unit, respectively. and The upper and lower limits of the output of the combined heat and power (CHP) unit; and These are the upper and lower ramp limits for CHP of a combined heat and power unit; The gas turbine GT generates electricity by burning natural gas, and its energy conversion relationship is shown in the following formula: Its operational constraints are as follows: In the formula, This refers to the gas consumption of the gas turbine GT. The gas-to-electric conversion efficiency of the gas turbine GT; and These are the upper and lower limits of the output of the gas turbine GT, respectively; and These are the uphill and downhill ramp limits for the gas turbine GT, respectively. The P2G (Power-to-Gas) equipment converts electrical energy into natural gas. The following formula represents the relationship between the gas production and electricity consumption of the P2G equipment: Its operational constraints are as follows: In the formula: For natural gas output from P2G (Power-to-Gas) equipment; and The upper and lower limits of the output of the P2G (Power-to-Gas) equipment; and The upper and lower limits of ramp rate for P2G electro-gas conversion equipment; The electric boiler EB meets the heat load demand by consuming electrical energy. The relationship between the power consumption and heat output of the electric boiler EB is shown in the following formula: Its constraints are as follows: In the formula: The heat output of the electric boiler EB; and The minimum and maximum power of the electric boiler EB; and The upper and lower limits of the ramp rate for electric boiler EB; The operational constraints for distributed wind power (WT) and distributed photovoltaic (PV) power are as follows: In the formula: and These are the maximum output power of distributed wind power (WT) and distributed photovoltaic (PV), respectively. and This indicates the reduction in distributed wind power (WT) and distributed photovoltaic (PV).

5. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 3, characterized in that, The power flow constraint for linearization of the distribution network is shown in the following equation: In the formula, and For the minimum and maximum power purchase quantities of multi-energy coupling aggregators; P ij,t It is the active power flow on the line; P ij,max This represents the maximum active power flow of the line; For electrical load, δ elec This is a preset safety margin constant used to provide sufficient operational buffer. Rewrite the energy balance formula using chance constraints: Where pr{} represents the probability of the event occurring in}, ζ is the confidence level, and γ is a function that follows a position parameter . A random variable with a logistic distribution and a scaling parameter of s; The deterministic equivalence method is used to handle opportunity constraints, transforming them into deterministic constraints, as shown in the following equation: K ζ =F -1 (1-g) Among them, F -1 It is the inverse function of the cumulative distribution function of the Logistic distribution.

6. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 3, characterized in that, The second-order cone relaxation constraint of the natural gas network: The distribution network and the natural gas network are connected through multi-energy coupling devices. Ignoring the influence of compressors, the steady-state natural gas network model includes node flow balance constraints and steady-state pipeline operation constraints. The node flow balance constraint is shown in the following equation: The Weymouth equation is shown below: The pipeline flow constraint is shown in the following formula: -G mn,max ≤G mn ≤G mn,max The nodal pressure constraint is shown in the following equation: p m,min +d gas ≤π m ≤π m,max -d gas p m,min ≤π m ≤π m,max The constraints on the purchase of natural gas by multi-energy coupled aggregators are shown in the following equation: In the formula: This refers to the amount of natural gas exchanged between the multi-energy coupling aggregator and the upstream natural gas network; and These represent the minimum and maximum quantities of natural gas that a multi-energy coupled aggregator can purchase; G mn For the current flow in pipe mn; Let π be the natural gas load at node m; m The air pressure at node m; m mn These are Weymouth feature parameters; G mn,max For the maximum current flow in pipe mn; π m,max and π m,min δ represents the minimum and maximum air pressures at node m, respectively. gas This is a preset safety margin constant used to provide sufficient operational buffer. Transforming the Weymouth equation into a second-order cone form, the problem is rewritten as a second-order cone programming problem: Introducing new variables: Weymouth's equation can be transformed into the following: (G mn ) 2 ≤τ mn Transforming it into a second-order cone constraint form is as follows: The heat network balance constraints: In the formula, The heat load within the multi-energy coupled polymer is η{heat}, which is the heat network transport efficiency coefficient and is a constant close to 1.

7. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 3, characterized in that, The multi-energy coupling aggregator constraint is expressed in compact form as follows: In the formula, x is a variable related to energy interaction and FRP, i.e. FRU t and FRD t y1 is a controllable variable under linear constraints, y2 is a controllable variable under nonlinear constraints, and the feasible region of x is a convex hull formed by the projections of all constraints in the model onto the x-axis. The mathematical expression for the feasible region is: For any x, there exists a set of corresponding scheduling schemes y1 and y2. The scheduling variables y1 and y2 for a fixed x′ are defined as follows:

8. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 1, characterized in that, The specific steps for constructing a max-min robust optimization model to identify the feasible region boundary are as follows: Introducing slack variables, construct the following second-order cone programming problem: stAx+By1-Iv≤b In the formula: v is a slack variable, which is used for... Perform the following transformations: (2y2) 2 +[(cy2+Iv)-1] 2 ≤[(cy2+Iv)+1] 2 The compact form is: ||Ey2+Iv+q i ||≤Fy2+Iv+d i The objective function of the second-order cone programming problem is to minimize the adjustment that satisfies all constraints, where f(x) = 0 indicates that x is a point in the feasible region and... Otherwise, if f(x) > 0, x is a point outside the feasible region and For a fixed x′, Y FR Since the condition is non-empty, the optimal value of the second-order cone programming problem is 0. The model is rewritten as follows: Ax + By1 - Iv ≤ b ||Ey2+Iv+q i ||≤Fy2+Iv+d i Introducing the concept of robust optimization, the model is reformulated as follows: stAx+By1-Iv≤b:π ||Ey2+Iv+q i ||≤Fy2+Iv+d i :(oh i ,l i )。 9. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 1, characterized in that, The method for transforming the max-min model into a mixed-integer second-order cone programming (MISOCP) model through mathematical transformation is as follows: Based on duality theory, the min problem within the model is transformed into a max problem for solution. The dual problem of the min problem within the model is: ||ω i ||≤λ i -1≤π≤0 The objective function in the above equation contains a bilinear term, and the constraints on x and π are independent. Based on robust optimization, the following form is obtained: For the max problem within the model: max(-Axπ) stDx≤d:u According to strong duality theory, we get: -Axπ=d·u The KKT conditions for the max problem within the model are: A T π+D T u=0 0≤u T ⊥(d-Dx)≤0 Constraint 0≤u T ⊥(d-Dx)≤0 is equivalent to 0≤u T , (d-Dx)≤0 and u T With the constraints (d-Dx) = 0, the mixed-integer second-order cone programming (MISOCP) model is obtained as follows: |||ω i ||≤λ i -1≤π≤0 A T π+D T u=0 0≤(d-Dx)≤M(δ-1) 0≤u T ≤Mδ。 10. The method for rapid identification of feasible regions of multi-energy coupled aggregators according to claim 1, characterized in that, The method for solving the MISOCP model using a polyhedral projection algorithm based on the bisection method and the cutting plane mechanism, and iteratively calculating and outputting the boundaries of the interactive energy feasible region (IEFR) and the flexible climbing feasible region (FRFR) is as follows: Step 1: For variable x, initialize a sufficiently large region. And select a point within the feasible region, denoted as x. in ; Step 2: By solving the mixed-integer second-order cone programming (MISOCP) model, obtain the optimal solution and the optimal value f1. In the worst case, select an external point x. out If f1(x) out If x > 0, it means x out Outside the feasible region, record x out And proceed to step 3; otherwise, f1(x) out ) = 0 means x out Within the feasible region, exit and log. Step 3: When x in Let x be an interior point of the feasible region. out When a point is an external point of the feasible region, the line connecting the internal and external points must intersect the boundary of the feasible region. Therefore, the bisection method is used to find the boundary point x of the feasible region. o An auxiliary parameter γ is introduced, which is obtained by using γ = 0.5*(γ1+γ2) and x o =(1-γ)·x in +γ·x out Determine x o Solve to determine x o Is it a feasible region boundary point when x o When determined, the corresponding π within the feasible region o , λ oi and ω oi This also led to the determination of how to record π. o , λ oi and ω oi ; ||ω′ i ||≤λ′ i -1≤h≤0 If f2(x) o )≠0, meaning x o For points that are not boundary points of the feasible region, update γ and x. o If f2(x) o If ) = 0, record the optimal solution to the above equation, which is π. o , λ oi and ω oi Then proceed to step 4; Step 4: Determine x o Then, using the known π o , λ oi and ω oi Add a boundary constraint as a cutting plane (b-Ax) to the mixed integer second-order cone programming (MISOCP) model. o ) T h-∑ i (λ oi d i +ω o T i F)≤0; Update the initial region Return to step; The iterative process from steps 2 to 4 continues until one of the following two conditions is met: Accuracy convergence: The volume change or maximum displacement of the boundary point of the polyhedral region defined by the cutting plane generated in two consecutive iterations is less than the preset convergence tolerance; Iteration convergence: The number of iterations reaches the preset maximum number of iterations.

Images