Industrial park integrated energy system stochastic planning method and device considering source load uncertainty

Abstract

The invention relates to an industrial park integrated energy system stochastic programming method considering source load uncertainty. The method comprises the following steps: obtaining a representative typical operation scene; establishing a carbon transaction price prediction model based on modal experience decomposition, long and short term memory and an improved grey wolf algorithm; an industrial park two-stage planning model considering the carbon transaction price is established and solved, and an optimal stochastic planning strategy is obtained. According to the method, the source load uncertainty can be effectively processed, various possible operation states in the industrial park can be comprehensively covered, meanwhile, the calculation complexity is reduced, and reliable input data is provided for a planning model; the dynamic change rule of the carbon transaction price can be accurately captured, and powerful support is provided for energy planning of an industrial park; the industrial park comprehensive energy system planning is optimized, the investment cost and the operation cost of the system can be effectively reduced, the energy utilization efficiency is improved, and the economical efficiency and the reliability of the industrial park comprehensive energy system are enhanced.

Description

Technical Field

[0001] This invention relates to the field of integrated energy system planning technology for industrial parks, and in particular to a stochastic planning method and equipment for integrated energy systems in industrial parks that takes into account source-load uncertainty. Background Technology

[0002] In recent decades, with the continuous development and breakthroughs in distributed generation, distributed energy storage, combined cooling, heating and power (CCHP), and power-to-gas (EPG) technologies, the coupling relationship between previously independently planned and operated electricity, natural gas, and heat networks is strengthening, gradually forming integrated energy systems (IES) that coordinate the supply of multiple energy forms. To break the traditional energy system pattern of "separate design, independent operation, and divide-and-conquer," in-depth research into the planning, design, and operation control technologies of integrated energy systems can effectively improve energy utilization efficiency and the proportion of new energy consumption, control greenhouse gas emissions, and thus promote the balanced and sustainable development of a low-carbon economy. In addition to building IESs, establishing a mature and fully functional carbon emissions trading market is also an effective means of controlling further greenhouse gas emissions. Introducing carbon trading costs into optimization models and achieving low-carbon goals through economic leverage is crucial.

[0003] Existing integrated energy system planning methods for industrial parks primarily employ deterministic models. These models, assuming fixed values ​​or typical daily curves for parameters such as renewable energy output and load demand, focus on optimizing the capacity configuration and operational strategies of multi-energy coupled equipment (electricity, gas, heat, etc.) within the system, prioritizing economic and technical efficiency. Even when carbon trading mechanisms are introduced in some studies, they are typically simplified to a fixed cost parameter. This approach suffers from two major flaws: First, it fails to effectively address the inherent randomness and volatility of energy sources and loads, potentially leading to insufficient reliability or poor economic performance in real-world operation, and a lack of robustness. Second, it ignores the volatility of carbon trading market prices, resulting in inaccurate assessments of the system's low-carbon benefits, potentially misleading investment decisions and hindering the achievement of a true optimal balance between economic efficiency and low carbon emissions in the context of energy transition.

[0004] Given the above background, developing a stochastic programming method for integrated energy systems in industrial parks that can account for source-load uncertainty and carbon trading price fluctuations is of significant practical importance. Summary of the Invention

[0005] To address the significant discrepancy between existing system planning and actual operational results, the primary objective of this invention is to provide a stochastic planning method for integrated energy systems in industrial parks that takes into account source-load uncertainty. This method can effectively handle source-load uncertainty, improve the accuracy of carbon trading price prediction, and optimize the planning of integrated energy systems in industrial parks.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a stochastic programming method for an integrated energy system in an industrial park that takes into account source-load uncertainty, the method comprising the following sequential steps:

[0007] (1) The Latin hypercube sampling method is adopted to generate the output scenario of the industrial park based on the cumulative distribution function of photovoltaic power generation, electrical load and heat load inside the industrial park. The scenario reduction method based on Manhattan distance is adopted to obtain a representative typical operation scenario.

[0008] (2) A carbon trading price prediction model is established based on modal empirical decomposition, long short-term memory and improved gray wolf algorithm;

[0009] (3) Based on typical operating scenarios and carbon trading price prediction models, establish and solve a two-stage planning model for industrial parks that takes into account carbon trading prices to obtain the optimal stochastic planning strategy.

[0010] Step (1) specifically includes the following steps in sequence:

[0011] (1a) The scene is generated using the Latin hypercube sampling method. Assume there are N independent random variables X1, X2...X NS Then the nth random variable X n The cumulative distribution function value is:

[0012] ;

[0013] In the formula: Let Y be the cumulative distribution function of the nth random variable. n Let X be the nth random variable. n The cumulative distribution function value;

[0014] Suppose we need to generate K output scenarios, then we divide the vertical axis of the cumulative distribution function of N random variables into K equal parts, and then randomly select a point from each small interval as the cumulative probability value for that interval. The calculation formula is as follows:

[0015] ;

[0016] In the formula: r k It is a random number, with a value range of 0-1, corresponding to the location of the sampling point in the small interval;

[0017] By using inverse function operations, the cumulative sampling probability value corresponding to the k-th subinterval is transformed into the sampled value X. nk The calculation formula is as follows:

[0018] ;

[0019] In the formula, It represents the inverse function of the cumulative distribution function;

[0020] By randomly sorting the K sampled values, the correlation of random variable samples is reduced, and K random scenarios are finally obtained.

[0021] (1b) Perform scene reduction based on Manhattan distance, which is defined as:

[0022] ;

[0023] In the formula, d gh Let represent the Manhattan distance between the g-th and h-th data groups; Λ is the total number of elements in the data group; and Let g and h be the values ​​of the elements at the u-th position, respectively;

[0024] By calculating the Manhattan distance for each pair of scene data, a corresponding distance matrix d is constructed. ge :

[0025] ;

[0026] Based on distance matrix d ge Calculate the sum of the Manhattan distances between any scene and the other scenes, and construct the probability matrix Y. ge The calculation formula is as follows:

[0027] ;

[0028] Based on probability matrix Y ge The two scenarios with the smallest sum of Manhattan distances are selected and reduced according to the following formula. This process is repeated until the typical operating scenario is obtained, and the probability of the typical operating scenario is calculated.

[0029] ;

[0030] In the formula, e and f are the scene numbers to be reduced; p e (t se ), p f (t se ) are respectively the tth se The probabilities of scenarios e and f in the next iteration.

[0031] Step (2) specifically includes the following steps in sequence:

[0032] (2a) Select three key structured data related to carbon trading price prediction—natural gas price, crude oil price, and coal price—as input data for the carbon trading price prediction model; select unstructured data related to carbon trading price prediction—carbon trading, carbon sinks, low-carbon economy, greenhouse effect, and integrated energy system—as output data for the carbon trading price prediction model; use linear interpolation to complete the structured data sequence, and use a dimensionality reduction strategy based on manifold learning to reduce the dimensionality of the unstructured data;

[0033] (2b) The normalized data of the input sequence is decomposed into n intrinsic mode function terms (IMF) and residuals (RES) by the empirical modal decomposition method (FCEEMD) to further mine the data information of the original sequence; let the original carbon trading price time series be p(t), and add white noise denoted as ω n (t), where the initial noise coefficient is represented by β0, the carbon trading price time series becomes p(t)+ω n (t)β0, the original signal is repeatedly decomposed N times using the EMD empirical mode decomposition method, and then the average value of the decomposition results is calculated according to the EEMD ensemble empirical mode decomposition method. The IMF component of p(t) is denoted as c1(t), and the remaining component of p(t) is taken as the first-order residual r1(t) as shown below:

[0034] ;

[0035] In the formula, N is the total number of IMFs or the number of repeated decompositions; E1 is the amplitude coefficient of the noise added in the first stage, E1=ε*σ0, σ0 is the standard deviation of p(t), and ε is a noise intensity coefficient with a value range of 0.1 to 0.4;

[0036] Continue to apply the signal r1(t)+E1ω n (t)β0 is decomposed N times, and the decomposed signal r k (t)+E k ω n (t)β0, decompose again and calculate c. k+1 (t) and k-th order residual r k (t), where k=2,…,K, as shown in the following equation:

[0037] (1);

[0038] In the formula, Let p(t) be the (k+1)th eigenmode function component. This is the residual signal after the k-th decomposition; β is the amplitude coefficient of the noise added in the k-th stage; k Let K be the noise figure for the k-th stage. The nth Gaussian white noise sequence is added.

[0039] Repeat formula (1) until the residuals cannot be further decomposed, at which point the entire decomposition process is terminated; K IMF components are obtained, and the final residual RES(t) is shown in the following formula:

[0040] ;

[0041] Therefore, the original carbon trading price time series p(t) is processed as shown in the following equation:

[0042] ;

[0043] (2c) Establishing a carbon trading price prediction model based on LSTM long short-term memory network: using the Sigmoid function to determine information retention and forgetting:

[0044] ;

[0045] In the formula, f t σ is the output of the forget gate at time t; σ is the sigmoid activation function; W f This is the weight matrix corresponding to the forget gate; [h t-1 , x t ] indicates that h t-1 and x t Concatenate two vectors; b f This corresponds to the bias term of the forget gate; x t Given the input at time t, h t-1 This is the hidden state from the previous moment;

[0046] h is adjusted using the Sigmoid function. t-1 and x t The output value i at time t is calculated. t Then use the tanh function to adjust h. t-1 and x t The calculation yields the output value. :

[0047] ;

[0048] In the formula, This is the weight matrix corresponding to the input gate; The bias term corresponding to the input gate; It is the hyperbolic tangent activation function; This is the weight matrix corresponding to the candidate cell states; This refers to the bias term corresponding to the candidate cell state;

[0049] Update the cell state at the current moment:

[0050] ;

[0051] h is calculated using the Sigmoid function. t-1 and x t Obtain the output value o at time t t Then use o t C after tanh operation t Multiply to obtain the hidden state h at the current time step. t :

[0052] ;

[0053] In the formula, These are the weights and biases of the output gate, respectively;

[0054] (2d) The parameters of the carbon trading price prediction model are optimized using an improved gray wolf algorithm. First, the wolf pack population hierarchy is subdivided into four levels: α1, α2, α3, and ω. The mathematical model for individual updates is as follows:

[0055] ;

[0056] In the formula, D is the absolute value of the distance vector between two individuals; A and C are both coefficient vectors; X is the individual's position vector; X b X is the position vector of the alpha wolf; p The position vector of the alpha wolf or prey used to guide position updates;

[0057] The expression for the coefficient vector in the above formula is:

[0058] ;

[0059] In the formula, R ram A random number between 0 and 1; T is the current iteration number; T max is the maximum number of iterations; e is the unit vector; a is the convergence operator;

[0060] The wolf pack position update equation is expressed as:

[0061] ;

[0062] ;

[0063] In the formula, , and These are the position vectors of the three wolves with the highest rank in the population at the Tth iteration; , and These are the coefficient vectors of the three-headed wolf, which is the highest-ranking wolf. , and X represents the coefficient vector of the three-headed wolf with the highest rank; X is the current position vector of the ordinary individual wolf being updated at the Tth iteration; X final (T) represents the final position of a normal individual; X1, X2, and X3 are the position vectors of the remaining individuals after being updated according to the position of the alpha wolf; the three alpha wolves with the highest rank refer to the three solutions with the best fitness.

[0064] The formula for calculating the introduced convergence operator 'a' is as follows:

[0065] ;

[0066] In the formula: a max a min These are the termination value and initial value of the convergence operator, respectively;

[0067] Combining improvements to the convergence operator, parameter optimization is performed using the value of |A| as a boundary;

[0068] When |A|≥1, the search direction of the subordinate population is corrected using the worst individual. The expression for the corrected individual is as follows:

[0069] ;

[0070] In the formula, X z To correct the individual; w The worst individual position; p * X is a dynamically decreasing factor; b Let alpha wolf's position vector be alpha.

[0071] The improved Grey Wolf algorithm is expressed as follows:

[0072] ;

[0073] (2e) Based on the improved Grey Wolf algorithm, the carbon trading price prediction model is optimized, each subsequence is predicted, and then the prediction components are integrated to obtain the final prediction result, that is, the predicted value of the original carbon trading price time series p(t) in the future.

[0074] In step (3), the two-stage planning model for industrial parks that takes into account carbon trading prices includes an upper-level planning model and a lower-level planning model;

[0075] (3a) The objective function of the upper-level programming model is , where C ann The comprehensive annualized cost is C. ann The formula is as follows:

[0076] ;

[0077] In the formula, C inv For equipment investment costs; C ope For operating costs; C main To maintain costs; C pun The cost of penalties for curtailing wind and solar power; C CO2 For carbon trading costs;

[0078] Equipment investment cost C inv Including the construction cost of expanded equipment, the expression is as follows:

[0079] ;

[0080] In the formula, The f-th model of newly built equipment for class a1; N g N represents the number of equipment types to be constructed. f This refers to the number of models for class q devices; The unit investment cost of the f-th model of the newly constructed equipment in category a1; r is the discount rate; N pl The planning period is in years;

[0081] Operating cost C ope The expression is as follows:

[0082] ;

[0083] In the formula, N y N is the number of days in a year, taken as 365; s N represents the number of typical days, categorized into summer, winter, and transitional season typical days; π(s) represents the probability of typical day s occurring; N t c is the number of time periods within a typical day; gas,t c represents the unit cost of natural gas at time t; elec,t P represents the unit cost of the external power grid at time t. gas,s,t P represents the input power of natural gas at time t on a typical day; elec,s,t The input power of the external power grid at time t on a typical day (s);

[0084] Maintenance cost C main As shown in the following formula:

[0085] ;

[0086] In the formula, P al,s,t c is the output power of device a1 at time t on a typical day; main,al The unit maintenance cost of equipment A1;

[0087] The cost of curtailing wind and solar power (C) pun The mathematical expression for it is shown below:

[0088] ;

[0089] In the formula, P pun,s,t c is the power of light discarded at time t during a typical day; pv,pun Cost per unit of curtailment;

[0090] Carbon trading cost C CO2 The mathematical expression for it is shown below:

[0091] ;

[0092] In the formula, P ψ,s,t P represents the output power of the carbon emission equipment ψ at time t on a typical day; P2G,s,t τ is the output power of the P2G unit at time t on a typical day; ψ,emis τ is the carbon quota coefficient for carbon emission equipment ψ; ψ,rat η is the carbon quota allocation coefficient for carbon emission equipment ψ; ab For the conversion efficiency of the P2G unit; θ CO2 For carbon trading prices; A collection of potential carbon-emitting devices;

[0093] (3b) After determining the initial capacity of the equipment, the lower-level planning model optimizes the operation and scheduling scheme of the park's integrated energy system according to the load demand of different typical days. The optimization variables include continuous variables such as the output power of the energy supply equipment and energy storage equipment, the energy input power and the curtailment of wind and solar power.

[0094] The objective function of the lower-level programming model is expressed as:

[0095] ;

[0096] In the formula: u1 and u2 are weighting coefficients; C ope,s The system operating cost per typical day is C. CO2,s The carbon trading cost for a typical day; C pun,s The cost of wind and solar power curtailment on a typical day; E CO2,s s represents the carbon emissions on a typical day; c gas and c elec For the unit price of natural gas and the unit price of electricity respectively; P gas,s,t and P elec,s,t These represent the gas consumption and electrical consumption at time t on a typical day, respectively; c ope,q Let P be the unit maintenance cost of equipment q; q,s (t) represents the power of a typical daily device q at time t; P pun,wt,s,t and P pun,pv,s,t These represent the wind and solar power curtailment at time t on a typical day s;

[0097] The equipment scheduling scheme calculated by the lower-level planning model is fed back into the upper-level planning model to further optimize the equipment planning scheme, thus completing one iterative calculation process; through the n-axis between the upper-level and lower-level planning models... max After several iterations, the final planning result of the integrated energy system, i.e., the optimal stochastic programming strategy, is obtained.

[0098] In step (3), the constraints of the two-stage planning model for industrial parks that takes into account carbon trading prices include power balance constraints, curtailment power constraints, energy storage equipment constraints, and energy supply equipment operation characteristic constraints.

[0099] (5a) The power balance constraints include supply and demand balance constraints on the electric bus side, supply and demand balance constraints on the natural gas bus side and supply and demand balance constraints on the thermal bus side.

[0100] The supply and demand balance constraints on the power bus side are as follows:

[0101] ;

[0102] In the formula, P pv,s,t Let s be the output power of PV at time t on a typical day; The electrical power output of CHP at time t on a typical day s; , These represent the charging and discharging power of ES at time t on a typical day s; Let be the power consumption of EB at time t on a typical day s; P represents the power consumption of the P2G at time t on a typical day; pun,s,t P represents the power of light discarded at time t on a typical day. E,s,t P represents the electrical load at time t on a typical day; elec,s,t The output power of the upstream power grid at time t on a typical day s; ES is electric energy storage; P2G is electricity-to-gas conversion; CHP is a combined cooling, heating and power system; EB is an electric boiler; PV is photovoltaic power generation;

[0103] The supply and demand balance constraints on the natural gas bus side are:

[0104] ;

[0105] In the formula, P gas,s,t The input power of the external natural gas network at time t on a typical day s; The gas production capacity of the P2G unit at time t on a typical day s; The gas consumption power of CHP at time t on a typical day s; The gas consumption power of GB at time t on a typical day s; GB is a gas-fired boiler;

[0106] The supply and demand balance constraints on the hot bus side are:

[0107] ;

[0108] In the formula, The thermal power output of EB at time t on a typical day s; The thermal power output of GB at time t on a typical day s; , These represent the charging and discharging power of HS at time t on a typical day s; P H,s,t denoted as GB, it represents the gas consumption power at time t on a typical day s; HS represents a thermal storage system. The thermal power output of CHP at time t on a typical day s;

[0109] (5b) The power constraint for wasted light is:

[0110] ;

[0111] In the formula, P pun,s,t The discarded power is the typical solar power at time t on day s; The theoretical maximum power that the PV can generate at time t on a typical day s;

[0112] (5c) The constraints on the energy storage device include energy storage state constraints and charge / discharge power constraints, wherein the energy storage state constraints are:

[0113] ;

[0114] The charge / discharge power constraint is:

[0115] ;

[0116] In the formula, SOC sto,s,t The energy storage state of the energy storage device at time t on a typical day; ρ sto,min ρ sto,max This represents the upper and lower limits of the remaining energy of the energy storage device. sto This refers to the rated capacity of the energy storage device. , , These are the standby efficiency, charging efficiency, and discharging efficiency of energy storage devices, respectively; P st,non L is the rated power of the energy storage device. sto for Rated capacity;

[0117] SOC sto,s,t The energy storage state of the energy storage device at time t on a typical day; ρ sto,min ρ sto,max L represents the upper and lower limits of the remaining energy of the energy storage device. sto This refers to the rated capacity of the energy storage device. , , These are the standby efficiency, charging efficiency, and discharging efficiency of energy storage devices, respectively. and These represent the energy storage power and energy release power of the energy storage device at time t on a typical day s;

[0118] (5d) The operating characteristics constraints of the power supply equipment are upper and lower limits of operating power:

[0119] ;

[0120] In the formula, P sup,max and P sup,min These represent the upper and lower limits of the power output of the energy supply equipment.

[0121] Another object of the present invention is to provide an electronic device comprising:

[0122] Processor; and

[0123] A memory storing computer program instructions, which, when executed by the processor, cause the processor to perform the stochastic planning method for an integrated energy system in an industrial park that takes into account source-load uncertainty, as described above.

[0124] The present invention also provides a computer-readable storage medium having stored thereon computer program instructions, which, when executed by a processor, cause the processor to perform the stochastic planning method for an integrated energy system in an industrial park that takes into account source-load uncertainty, as described above.

[0125] As can be seen from the above technical solution, the beneficial effects of the present invention are as follows: First, the present invention can effectively handle source-load uncertainty: through Latin hypercube sampling and Manhattan distance scenario reduction methods, it can comprehensively cover all possible operating states within the industrial park, while reducing computational complexity and providing reliable input data for the planning model; Second, the present invention improves the accuracy of carbon trading price prediction: the prediction model based on modal empirical decomposition, LSTM, and the improved gray wolf algorithm can accurately capture the dynamic change law of carbon trading prices, providing strong support for energy planning in industrial parks; Third, the present invention optimizes the planning of integrated energy systems in industrial parks: the two-stage planning model comprehensively considers long-term planning and short-term operational optimization, which can effectively reduce the investment and operating costs of the system, improve energy utilization efficiency, and enhance the economy and reliability of the integrated energy system in industrial parks. Attached Figure Description

[0126] Figure 1 This is a flowchart of the method of the present invention;

[0127] Figure 2 This is a flowchart of the carbon trading price prediction method in this invention. Detailed Implementation

[0128] like Figure 1 As shown, a stochastic programming method for an integrated energy system in an industrial park that takes into account source-load uncertainty is presented. The method includes the following sequential steps:

[0129] (1) The Latin hypercube sampling method is adopted to generate the output scenario of the industrial park based on the cumulative distribution function of photovoltaic power generation, electrical load and heat load inside the industrial park. The scenario reduction method based on Manhattan distance is adopted to obtain a representative typical operation scenario.

[0130] (2) A carbon trading price prediction model is established based on modal empirical decomposition, long short-term memory and improved gray wolf algorithm;

[0131] (3) Based on typical operating scenarios and carbon trading price prediction models, establish and solve a two-stage planning model for industrial parks that takes into account carbon trading prices to obtain the optimal stochastic planning strategy.

[0132] Step (1) specifically includes the following steps in sequence:

[0133] (1a) The scene is generated using the Latin hypercube sampling method. Assume there are N independent random variables X1, X2...X NS Then the nth random variable X n The cumulative distribution function value is:

[0134] ;

[0135] In the formula: Let Y be the cumulative distribution function of the nth random variable. n Let X be the nth random variable. n The cumulative distribution function value;

[0136] Suppose we need to generate K output scenarios, then we divide the vertical axis of the cumulative distribution function of N random variables into K equal parts, and then randomly select a point from each small interval as the cumulative probability value for that interval. The calculation formula is as follows:

[0137] ;

[0138] In the formula: r k It is a random number, with a value range of 0-1, corresponding to the location of the sampling point in the small interval;

[0139] By using inverse function operations, the cumulative sampling probability value corresponding to the k-th subinterval is transformed into the sampled value X. nk The calculation formula is as follows:

[0140] ;

[0141] In the formula, It represents the inverse function of the cumulative distribution function;

[0142] By randomly sorting the K sampled values, the correlation of random variable samples is reduced, and K random scenarios are finally obtained.

[0143] (1b) Perform scene reduction based on Manhattan distance, which is defined as:

[0144] ;

[0145] In the formula, d gh Let represent the Manhattan distance between the g-th and h-th data groups; Λ is the total number of elements in the data group; and Let g and h be the values ​​of the elements at the u-th position, respectively;

[0146] By calculating the Manhattan distance for each pair of scene data, a corresponding distance matrix d is constructed. ge :

[0147] ;

[0148] Based on distance matrix d ge Calculate the sum of the Manhattan distances between any scene and the other scenes, and construct the probability matrix Y. ge The calculation formula is as follows:

[0149] ;

[0150] Based on probability matrix Y ge The two scenarios with the smallest sum of Manhattan distances are selected and reduced according to the following formula. This process is repeated until the typical operating scenario is obtained, and the probability of the typical operating scenario is calculated.

[0151] ;

[0152] In the formula, e and f are the scene numbers to be reduced; p e (t se ), p f (t se ) are respectively the tth se The probabilities of scenarios e and f in the next iteration.

[0153] like Figure 2 As shown, step (2) specifically includes the following steps in sequence:

[0154] (2a) Select three key structured data related to carbon trading price prediction—natural gas price, crude oil price, and coal price—as input data for the carbon trading price prediction model; select unstructured data related to carbon trading price prediction—carbon trading, carbon sinks, low-carbon economy, greenhouse effect, and integrated energy system—as output data for the carbon trading price prediction model; use linear interpolation to complete the structured data sequence, and use a dimensionality reduction strategy based on manifold learning to reduce the dimensionality of the unstructured data;

[0155] (2b) The normalized data of the input sequence is decomposed into n intrinsic mode function terms (IMF) and residuals (RES) by the empirical modal decomposition method (FCEEMD) to further mine the data information of the original sequence; let the original carbon trading price time series be p(t), and add white noise denoted as ω n (t), where the initial noise coefficient is represented by β0, the carbon trading price time series becomes p(t)+ω n (t)β0, the original signal is repeatedly decomposed N times using the EMD empirical mode decomposition method, and then the average value of the decomposition results is calculated according to the EEMD ensemble empirical mode decomposition method. The IMF component of p(t) is denoted as c1(t), and the remaining component of p(t) is taken as the first-order residual r1(t) as shown below:

[0156] ;

[0157] In the formula, N is the total number of IMFs or the number of repeated decompositions; E1 is the amplitude coefficient of the noise added in the first stage, E1=ε*σ0, σ0 is the standard deviation of p(t), and ε is a noise intensity coefficient with a value range of 0.1 to 0.4;

[0158] Continue to apply the signal r1(t)+E1ω n (t)β0 is decomposed N times, and the decomposed signal r k (t)+E k ω n (t)β0, decompose again and calculate c. k+1 (t) and k-th order residual r k (t), where k=2,…,K, as shown in the following equation:

[0159] (1);

[0160] In the formula, Let p(t) be the (k+1)th eigenmode function component. This is the residual signal after the k-th decomposition; β is the amplitude coefficient of the noise added in the k-th stage; k Let K be the noise figure for the k-th stage. The nth Gaussian white noise sequence is added.

[0161] Repeat formula (1) until the residuals cannot be further decomposed, at which point the entire decomposition process is terminated; K IMF components are obtained, and the final residual RES(t) is shown in the following formula:

[0162] ;

[0163] Therefore, the original carbon trading price time series p(t) is processed as shown in the following equation:

[0164] ;

[0165] (2c) Establishing a carbon trading price prediction model based on LSTM long short-term memory network: using the Sigmoid function to determine information retention and forgetting:

[0166] ;

[0167] In the formula, f t σ is the output of the forget gate at time t; σ is the sigmoid activation function; W f This is the weight matrix corresponding to the forget gate; [h t-1 , x t ] indicates that h t-1 and x t Concatenate two vectors; b f This corresponds to the bias term of the forget gate; x t Given the input at time t, h t-1 This is the hidden state from the previous moment;

[0168] h is adjusted using the Sigmoid function. t-1 and x t The output value i at time t is calculated. t Then use the tanh function to adjust h. t-1 and x t The calculation yields the output value. :

[0169] ;

[0170] In the formula, This is the weight matrix corresponding to the input gate; The bias term corresponding to the input gate; It is the hyperbolic tangent activation function; This is the weight matrix corresponding to the candidate cell states; This refers to the bias term corresponding to the candidate cell state;

[0171] The cell state C at the previous moment t-1 The vector f output by the forget gatet Perform point-by-point multiplication, if f t If the value of f is close to 0, the corresponding information will be discarded; if f t If the value of f is close to 1, the information is preserved. t Between (0, 1), information is not fully preserved, but only partially retained. Then, the current cell state is updated by adding the cell state value from the previous time step to the output of the input gate point by point:

[0172] ;

[0173] h is calculated using the Sigmoid function. t-1 and x t Obtain the output value o at time t t Then use o t C after tanh operation t Multiply to obtain the hidden state h at the current time step. t :

[0174] ;

[0175] In the formula, These are the weights and biases of the output gate, respectively;

[0176] (2d) The parameters of the carbon trading price prediction model are optimized using an improved gray wolf algorithm. First, the wolf pack population hierarchy is subdivided into four levels: α1, α2, α3, and ω. The mathematical model for individual updates is as follows:

[0177] ;

[0178] In the formula, D is the absolute value of the distance vector between two individuals; A and C are both coefficient vectors; X is the individual's position vector; X b X is the position vector of the alpha wolf; p The position vector of the alpha wolf or prey used to guide position updates;

[0179] The expression for the coefficient vector in the above formula is:

[0180] ;

[0181] In the formula, R ram A random number between 0 and 1; T is the current iteration number; T max is the maximum number of iterations; e is the unit vector; a is the convergence operator;

[0182] The wolf pack position update equation is expressed as:

[0183] ;

[0184] ;

[0185] In the formula, , and These are the position vectors of the three wolves with the highest rank in the population at the Tth iteration; , and These are the coefficient vectors of the three-headed wolf, which is the highest-ranking wolf. , and X represents the coefficient vector of the three-headed wolf with the highest rank; X is the current position vector of the ordinary individual wolf being updated at the Tth iteration; X final (T) represents the final position of a normal individual; X1, X2, and X3 are the position vectors of the remaining individuals after being updated according to the position of the alpha wolf; the three alpha wolves with the highest rank refer to the three solutions with the best fitness.

[0186] The formula for calculating the introduced convergence operator 'a' is as follows:

[0187] ;

[0188] In the formula: a max a min These are the termination value and initial value of the convergence operator, respectively;

[0189] Combining improvements to the convergence operator, parameter optimization is performed using the value of |A| as a boundary;

[0190] When |A|≥1, the search direction of the subordinate population is corrected using the worst individual. The expression for the corrected individual is as follows:

[0191] ;

[0192] In the formula, X z To correct the individual; w The worst individual position; p * X is a dynamically decreasing factor; b Let alpha wolf's position vector be alpha.

[0193] The improved Grey Wolf algorithm is expressed as follows:

[0194] ;

[0195] (2e) Based on the improved Grey Wolf algorithm, the carbon trading price prediction model is optimized, each subsequence is predicted, and then the prediction components are integrated to obtain the final prediction result, that is, the predicted value of the original carbon trading price time series p(t) in the future.

[0196] In step (3), the two-stage planning model for industrial parks that takes into account carbon trading prices includes an upper-level planning model and a lower-level planning model;

[0197] (3a) The objective function of the upper-level programming model is , where C ann The comprehensive annualized cost is C. ann The formula is as follows:

[0198] ;

[0199] In the formula, C inv For equipment investment costs; C ope For operating costs; C main To maintain costs; C pun The cost of penalties for curtailing wind and solar power; C CO2 For carbon trading costs;

[0200] Equipment investment cost C inv Including the construction cost of expanded equipment, the expression is as follows:

[0201] ;

[0202] In the formula, The f-th model of newly built equipment for class a1; N g N represents the number of equipment types to be constructed. f This refers to the number of models for class q devices; The unit investment cost of the f-th model of the newly constructed equipment in category a1; r is the discount rate; N pl The planning period is in years;

[0203] Operating cost C ope The expression is as follows:

[0204] ;

[0205] In the formula, N y N is the number of days in a year, taken as 365; s N represents the number of typical days, categorized into summer, winter, and transitional season typical days; π(s) represents the probability of typical day s occurring; N t c is the number of time periods within a typical day; gas,t c represents the unit cost of natural gas at time t; elec,t P represents the unit cost of the external power grid at time t. gas,s,t P represents the input power of natural gas at time t on a typical day; elec,s,t The input power of the external power grid at time t on a typical day (s);

[0206] Maintenance cost C main As shown in the following formula:

[0207] ;

[0208] In the formula, P al,s,t c is the output power of device a1 at time t on a typical day; main,al The unit maintenance cost of equipment A1;

[0209] The cost of curtailing wind and solar power (C) pun The mathematical expression for it is shown below:

[0210] ;

[0211] In the formula, P pun,s,t c is the power of light discarded at time t during a typical day; pv,pun Cost per unit of curtailment;

[0212] Carbon trading cost C CO2 The mathematical expression for it is shown below:

[0213] ;

[0214] In the formula, P ψ,s,t P represents the output power of the carbon emission equipment ψ at time t on a typical day; P2G,s,t τ is the output power of the P2G unit at time t on a typical day; ψ,emis τ is the carbon quota coefficient for carbon emission equipment ψ; ψ,rat η is the carbon quota allocation coefficient for carbon emission equipment ψ; ab For the conversion efficiency of the P2G unit; θ CO2 For carbon trading prices; A collection of potential carbon-emitting devices;

[0215] (3b) After determining the initial capacity of the equipment, the lower-level planning model optimizes the operation and scheduling scheme of the park's integrated energy system according to the load demand of different typical days. The optimization variables include continuous variables such as the output power of the energy supply equipment and energy storage equipment, the energy input power and the curtailment of wind and solar power.

[0216] The objective function of the lower-level programming model is expressed as:

[0217] ;

[0218] In the formula: u1 and u2 are weighting coefficients; C ope,s The system operating cost per typical day is C. CO2,s The carbon trading cost for a typical day; C pun,s The cost of wind and solar power curtailment on a typical day; E CO2,s s represents the carbon emissions on a typical day; c gas and celec For the unit price of natural gas and the unit price of electricity respectively; P gas,s,t and P elec,s,t These represent the gas consumption and electrical consumption at time t on a typical day, respectively; c ope,q Let P be the unit maintenance cost of equipment q; q,s (t) represents the power of a typical daily device q at time t; P pun,wt,s,t and P pun,pv,s,t These represent the wind and solar power curtailment at time t on a typical day s;

[0219] The equipment scheduling scheme calculated by the lower-level planning model is fed back into the upper-level planning model to further optimize the equipment planning scheme, thus completing one iterative calculation process; through the n-axis between the upper-level and lower-level planning models... max After several iterations, the final planning result of the integrated energy system, i.e., the optimal stochastic programming strategy, is obtained.

[0220] In step (3), the constraints of the two-stage planning model for industrial parks that takes into account carbon trading prices include power balance constraints, curtailment power constraints, energy storage equipment constraints, and energy supply equipment operation characteristic constraints.

[0221] (5a) The power balance constraints include supply and demand balance constraints on the electric bus side, supply and demand balance constraints on the natural gas bus side and supply and demand balance constraints on the thermal bus side.

[0222] The supply and demand balance constraints on the power bus side are as follows:

[0223] ;

[0224] In the formula, P pv,s,t Let s be the output power of PV at time t on a typical day; The electrical power output of CHP at time t on a typical day s; , These represent the charging and discharging power of ES at time t on a typical day s; Let be the power consumption of EB at time t on a typical day s; P represents the power consumption of the P2G at time t on a typical day; pun,s,t P represents the power of light discarded at time t on a typical day. E,s,t P represents the electrical load at time t on a typical day; elec,s,t The output power of the upstream power grid at time t on a typical day s; ES is electric energy storage; P2G is electricity-to-gas conversion; CHP is a combined cooling, heating and power system; EB is an electric boiler; PV is photovoltaic power generation;

[0225] The primary energy source for the natural gas busbar is the external gas network, with the main energy-consuming equipment being GB and P2G units. Additionally, due to the operating characteristics of the P2G units, the produced natural gas is transmitted back to the main busbar for redistribution. The supply and demand balance constraints on the natural gas busbar side are:

[0226] ;

[0227] In the formula, P gas,s,t The input power of the external natural gas network at time t on a typical day s; The gas production capacity of the P2G unit at time t on a typical day s; The gas consumption power of CHP at time t on a typical day s; The gas consumption power of GB at time t on a typical day s; GB is a gas-fired boiler;

[0228] The supply and demand balance constraints on the hot bus side are:

[0229] ;

[0230] In the formula, The thermal power output of EB at time t on a typical day s; The thermal power output of GB at time t on a typical day s; , These represent the charging and discharging power of HS at time t on a typical day s; P H,s,t denoted as GB, it represents the gas consumption power at time t on a typical day s; HS represents a thermal storage system. The thermal power output of CHP at time t on a typical day s;

[0231] (5b) The power constraint for wasted light is:

[0232] ;

[0233] In the formula, P pun,s,t The discarded power is the typical solar power at time t on day s; The theoretical maximum power that the PV can generate at time t on a typical day s;

[0234] (5c) The constraints on the energy storage device include energy storage state constraints and charge / discharge power constraints, wherein the energy storage state constraints are:

[0235] ;

[0236] The charge / discharge power constraint is:

[0237] ;

[0238] In the formula, SOC sto,s,t The energy storage state of the energy storage device at time t on a typical day; ρsto,min ρ sto,max This represents the upper and lower limits of the remaining energy of the energy storage device. sto This refers to the rated capacity of the energy storage device. , , These are the standby efficiency, charging efficiency, and discharging efficiency of energy storage devices, respectively; P st,non L is the rated power of the energy storage device. sto for Rated capacity;

[0239] SOC sto,s,t The energy storage state of the energy storage device at time t on a typical day; ρ sto,min ρ sto,max L represents the upper and lower limits of the remaining energy of the energy storage device. sto This refers to the rated capacity of the energy storage device. , , These are the standby efficiency, charging efficiency, and discharging efficiency of energy storage devices, respectively. and These represent the energy storage power and energy release power of the energy storage device at time t on a typical day s;

[0240] (5d) The operating characteristics constraints of the power supply equipment are upper and lower limits of operating power:

[0241] ;

[0242] In the formula, P sup,max and P sup,min These represent the upper and lower limits of the power output of the energy supply equipment.

[0243] In summary, this invention effectively addresses source-load uncertainty: by employing Latin hypercube sampling and Manhattan distance scenario reduction methods, it comprehensively covers various possible operating states within industrial parks while reducing computational complexity, providing reliable input data for planning models. This invention also improves the accuracy of carbon trading price prediction: the prediction model based on modal empirical decomposition, LSTM, and the improved Grey Wolf algorithm accurately captures the dynamic changes in carbon trading prices, providing strong support for energy planning in industrial parks. Furthermore, this invention optimizes the planning of integrated energy systems in industrial parks: the two-stage planning model comprehensively considers long-term planning and short-term operational optimization, effectively reducing system investment and operating costs, improving energy efficiency, and enhancing the economic viability and reliability of integrated energy systems in industrial parks.

[0244] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the claimed invention. The scope of protection claimed by the appended claims and their equivalents is defined.

Claims

1. A stochastic programming method for an integrated energy system in an industrial park that takes into account source-load uncertainty, characterized in that: The method includes the following steps in sequence: (1) The Latin hypercube sampling method is adopted to generate the output scenario of the industrial park based on the cumulative distribution function of photovoltaic power generation, electrical load and heat load inside the industrial park. The scenario reduction method based on Manhattan distance is adopted to obtain a representative typical operation scenario. (2) A carbon trading price prediction model is established based on modal empirical decomposition, long short-term memory and improved gray wolf algorithm; (3) Based on typical operating scenarios and carbon trading price prediction models, establish and solve a two-stage planning model for industrial parks that takes into account carbon trading prices to obtain the optimal stochastic planning strategy.

2. The stochastic programming method for integrated energy systems in industrial parks that takes into account source-load uncertainty, as described in claim 1, is characterized in that: Step (1) specifically includes the following steps in sequence: (1a) The scene is generated using the Latin hypercube sampling method. Assume there are N independent random variables X1, X2...X NS Then the nth random variable X n The cumulative distribution function value is: ； In the formula: Let Y be the cumulative distribution function of the nth random variable. n Let X be the nth random variable. n The cumulative distribution function value; Suppose we need to generate K output scenarios, then we divide the vertical axis of the cumulative distribution function of N random variables into K equal parts, and then randomly select a point from each small interval as the cumulative probability value for that interval. The calculation formula is as follows: ； In the formula: r k It is a random number, with a value range of 0-1, corresponding to the location of the sampling point in the small interval; By using inverse function operations, the cumulative sampling probability value corresponding to the k-th subinterval is transformed into the sampled value X. nk The calculation formula is as follows: ； In the formula, It represents the inverse function of the cumulative distribution function; By randomly sorting the K sampled values, the correlation of random variable samples is reduced, and K random scenarios are finally obtained. (1b) Perform scene reduction based on Manhattan distance, which is defined as: ； In the formula, d gh Let represent the Manhattan distance between the g-th and h-th data groups; Λ is the total number of elements in the data group; and Let g and h be the values ​​of the elements at the u-th position, respectively; By calculating the Manhattan distance for each pair of scene data, a corresponding distance matrix d is constructed. ge : ； Based on distance matrix d ge Calculate the sum of the Manhattan distances between any scene and the other scenes, and construct the probability matrix Y. ge The calculation formula is as follows: ； Based on probability matrix Y ge The two scenarios with the smallest sum of Manhattan distances are selected and reduced according to the following formula. This process is repeated until the typical operating scenario is obtained, and the probability of the typical operating scenario is calculated. ； In the formula, e and f are the scene numbers to be reduced; p e (t se ), p f (t se ) are respectively the tth se The probabilities of scenarios e and f in the next iteration.

3. The stochastic programming method for integrated energy systems in industrial parks that takes into account source-load uncertainty, as described in claim 1, is characterized in that: Step (2) specifically includes the following steps in sequence: (2a) Select three key structured data related to carbon trading price prediction—natural gas price, crude oil price, and coal price—as input data for the carbon trading price prediction model; select unstructured data related to carbon trading price prediction—carbon trading, carbon sinks, low-carbon economy, greenhouse effect, and integrated energy system—as output data for the carbon trading price prediction model; use linear interpolation to complete the structured data sequence, and use a dimensionality reduction strategy based on manifold learning to reduce the dimensionality of the unstructured data; (2b) The normalized data of the input sequence is decomposed into n intrinsic mode function terms (IMF) and residuals (RES) by the empirical modal decomposition method (FCEEMD) to further mine the data information of the original sequence; let the original carbon trading price time series be p(t), and add white noise denoted as ω n (t), where the initial noise coefficient is represented by β0, the carbon trading price time series becomes p(t)+ω n (t)β0, the original signal is repeatedly decomposed N times using the EMD empirical mode decomposition method, and then the average value of the decomposition results is calculated according to the EEMD ensemble empirical mode decomposition method. The IMF component of p(t) is denoted as c1(t), and the remaining component of p(t) is taken as the first-order residual r1(t) as shown below: ； In the formula, N is the total number of IMFs or the number of repeated decompositions; E1 is the amplitude coefficient of the noise added in the first stage, E1=ε*σ0, σ0 is the standard deviation of p(t), and ε is a noise intensity coefficient with a value range of 0.1 to 0.4; Continue with the signal r1(t)+E1ω n (t)β0 is decomposed N times, and the decomposed signal r k (t)+E k ω n (t)β0, decompose again and calculate c. k+1 (t) and k-th order residual r k (t), where k=2,…,K, as shown in the following equation: （1）； In the formula, Let p(t) be the (k+1)th eigenmode function component. This is the residual signal after the k-th decomposition; β is the amplitude coefficient of the noise added in the k-th stage; k Let K be the noise figure for the k-th stage. The nth Gaussian white noise sequence added; Repeat formula (1) until the residuals cannot be further decomposed, at which point the entire decomposition process is terminated; K IMF components are obtained, and the final residual RES(t) is shown in the following formula: ； Therefore, the original carbon trading price time series p(t) is processed as shown in the following equation: ； (2c) Establishing a carbon trading price prediction model based on LSTM long short-term memory network: using the Sigmoid function to determine information retention and forgetting: ； In the formula, f t σ is the output of the forget gate at time t; σ is the sigmoid activation function; W f This is the weight matrix corresponding to the forget gate; [h t-1 , x t ] indicates that h t-1 and x t Concatenate two vectors; b f This corresponds to the bias term of the forget gate; x t Given the input at time t, h t-1 This is the hidden state from the previous moment; h is adjusted using the Sigmoid function. t-1 and x t The output value i at time t is calculated. t Then use the tanh function to adjust h. t-1 and x t The calculation yields the output value. : ； In the formula, This is the weight matrix corresponding to the input gate; The bias term corresponding to the input gate; It is the hyperbolic tangent activation function; This is the weight matrix corresponding to the candidate cell states; This refers to the bias term corresponding to the candidate cell state; Update the cell state at the current moment: ； h is calculated using the Sigmoid function. t-1 and x t Obtain the output value o at time t t Then use o t C after tanh operation t Multiply to obtain the hidden state h at the current time step. t : ； In the formula, These are the weights and biases of the output gate, respectively; (2d) The parameters of the carbon trading price prediction model are optimized using an improved gray wolf algorithm. First, the wolf pack population hierarchy is subdivided into four levels: α1, α2, α3, and ω. The mathematical model for individual updates is as follows: ； In the formula, D is the absolute value of the distance vector between two individuals; A and C are both coefficient vectors; X is the individual's position vector; X b X is the position vector of the alpha wolf; p The position vector of the alpha wolf or prey used to guide position updates; The expression for the coefficient vector in the above formula is: ； In the formula, R ram A random number between 0 and 1; T is the current iteration number; T max is the maximum number of iterations; e is the unit vector; a is the convergence operator; The wolf pack position update equation is expressed as: ； ； In the formula, , and These are the position vectors of the three wolves with the highest rank in the population at the Tth iteration; , and These are the coefficient vectors of the three-headed wolf, which is the highest-ranking wolf. , and X represents the coefficient vector of the three-headed wolf with the highest rank; X is the current position vector of the ordinary individual wolf being updated at the Tth iteration; X final (T) represents the final position of a normal individual; X1, X2, and X3 are the position vectors of the remaining individuals after being updated according to the position of the alpha wolf; the three alpha wolves with the highest rank refer to the three solutions with the best fitness. The formula for calculating the introduced convergence operator 'a' is as follows: ； In the formula: a max a min These are the termination value and initial value of the convergence operator, respectively; Combining improvements to the convergence operator, parameter optimization is performed using the value of |A| as a boundary; When |A|≥1, the search direction of the subordinate population is corrected using the worst individual. The expression for the corrected individual is as follows: ； In the formula, X z To correct the individual; w The worst individual position; p * X is a dynamically decreasing factor; b Let alpha wolf's position vector be alpha. The improved Grey Wolf algorithm is expressed as follows: ； (2e) Based on the improved Grey Wolf algorithm, the carbon trading price prediction model is optimized, each subsequence is predicted, and then the prediction components are integrated to obtain the final prediction result, that is, the predicted value of the original carbon trading price time series p(t) in the future.

4. The stochastic programming method for integrated energy systems in industrial parks that takes into account source-load uncertainty, as described in claim 1, is characterized in that: In step (3), the two-stage planning model for industrial parks that takes into account carbon trading prices includes an upper-level planning model and a lower-level planning model; (3a) The objective function of the upper-level programming model is , where C ann The comprehensive annualized cost is C. ann The formula is as follows: ； In the formula, C inv For equipment investment costs; C ope For operating costs; C main To maintain costs; C pun The cost of penalties for curtailing wind and solar power; C CO2 For carbon trading costs; Equipment investment cost C inv Including the construction cost of expanded equipment, the expression is as follows: ； In the formula, The f-th model of newly built equipment for class a1; N g N represents the number of equipment types to be constructed. f This refers to the number of models for class q devices; The unit investment cost of the f-th model of the newly constructed equipment in category a1; r is the discount rate; N pl The planning period is in years; Operating cost C ope The expression is as follows: ； In the formula, N y N is the number of days in a year, taken as 365; s N represents the number of typical days, categorized into summer, winter, and transitional season typical days; π(s) represents the probability of typical day s occurring; N t c is the number of time periods within a typical day; gas,t c represents the unit cost of natural gas at time t; elec,t P represents the unit cost of the external power grid at time t. gas,s,t P represents the input power of natural gas at time t on a typical day; elec,s,t The input power of the external power grid at time t on a typical day (s); Maintenance cost C main As shown in the following formula: ； In the formula, P al,s,t c is the output power of device a1 at time t on a typical day; main,al The unit maintenance cost of equipment A1; The cost of curtailing wind and solar power (C) pun The mathematical expression for it is shown below: ； In the formula, P pun,s,t c is the power of light discarded at time t during a typical day; pv,pun Cost per unit of curtailment; Carbon trading cost C CO2 The mathematical expression for it is shown below: ； In the formula, P ψ,s,t P represents the output power of the carbon emission equipment ψ at time t on a typical day; P2G,s,t τ is the output power of the P2G unit at time t on a typical day; ψ,emis τ is the carbon quota coefficient for carbon emission equipment ψ; ψ,rat η is the carbon quota allocation coefficient for carbon emission equipment ψ; ab For the conversion efficiency of the P2G unit; θ CO2 For carbon trading prices; A collection of potential carbon-emitting devices; (3b) After determining the initial capacity of the equipment, the lower-level planning model optimizes the operation and scheduling scheme of the park's integrated energy system according to the load demand of different typical days. The optimization variables include continuous variables such as the output power of the energy supply equipment and energy storage equipment, the energy input power and the curtailment of wind and solar power. The objective function of the lower-level programming model is expressed as: ； In the formula: u1 and u2 are weighting coefficients; C ope,s The system operating cost per typical day is C. CO2,s The carbon trading cost for a typical day; C pun,s The cost of wind and solar power curtailment on a typical day; E CO2,s s represents the carbon emissions on a typical day; c gas and c elec For the unit price of natural gas and the unit price of electricity respectively; P gas,s,t and P elec,s,t These represent the gas consumption and electrical consumption at time t on a typical day, respectively; c ope,q Let P be the unit maintenance cost of equipment q; q,s (t) represents the power of a typical daily device q at time t; P pun,wt,s,t and P pun,pv,s,t These represent the wind and solar power curtailment at time t on a typical day s; The equipment scheduling scheme calculated by the lower-level planning model is fed back into the upper-level planning model to further optimize the equipment planning scheme, thus completing one iterative calculation process; through the n-axis between the upper-level and lower-level planning models... max After several iterations, the final planning result of the integrated energy system, i.e., the optimal stochastic programming strategy, is obtained.

5. The stochastic programming method for integrated energy systems in industrial parks considering source-load uncertainty as described in claim 1, characterized in that: In step (3), the constraints of the two-stage planning model for industrial parks that takes into account carbon trading prices include power balance constraints, curtailment power constraints, energy storage equipment constraints, and energy supply equipment operation characteristic constraints. (5a) The power balance constraints include supply and demand balance constraints on the electric bus side, supply and demand balance constraints on the natural gas bus side and supply and demand balance constraints on the thermal bus side. The supply and demand balance constraints on the power bus side are as follows: ； In the formula, P pv,s,t Let s be the output power of PV at time t on a typical day; The electrical power output of CHP at time t on a typical day s; , These represent the charging and discharging power of ES at time t on a typical day s; Let be the power consumption of EB at time t on a typical day s; P represents the power consumption of the P2G at time t on a typical day; pun,s,t P represents the power of light discarded at time t on a typical day. E,s,t P represents the electrical load at time t on a typical day; elec,s,t The output power of the upstream power grid at time t on a typical day s; ES is electric energy storage; P2G is electricity-to-gas conversion; CHP is a combined cooling, heating and power system; EB is an electric boiler; PV is photovoltaic power generation; The supply and demand balance constraints on the natural gas bus side are: ； In the formula, P gas,s,t The input power of the external natural gas network at time t on a typical day s; The gas production capacity of the P2G unit at time t on a typical day s; The gas consumption power of CHP at time t on a typical day s; The gas consumption of GB is defined as the gas power consumption of a typical day at time t; GB refers to a gas-fired boiler. The supply and demand balance constraints on the hot bus side are: ； In the formula, The thermal power output of EB at time t on a typical day s; The thermal power output of GB at time t on a typical day s; , These represent the charging and dissipating heat power of HS at time t on a typical day s; P H,s,t is the gas consumption power of GB at time t on a typical day s; HS is a thermal storage system; The thermal power output of CHP at time t on a typical day s; (5b) The power constraint for wasted light is: ； In the formula, P pun,s,t The discarded power is the typical solar power at time t on day s; The theoretical maximum power that the PV can generate at time t on a typical day s; (5c) The constraints on the energy storage device include energy storage state constraints and charge / discharge power constraints, wherein the energy storage state constraints are: ； The charge / discharge power constraint is: ； In the formula, SOC sto,s,t The energy storage state of the energy storage device at time t on a typical day; ρ sto,min ρ sto,max This represents the upper and lower limits of the remaining energy of the energy storage device. sto This refers to the rated capacity of the energy storage device. , , These are the standby efficiency, charging efficiency, and discharging efficiency of energy storage devices, respectively; P st,non L is the rated power of the energy storage device. sto for Rated capacity; SOC sto,s,t The energy storage state of the energy storage device at time t on a typical day; ρ sto,min ρ sto,max L represents the upper and lower limits of the remaining energy of the energy storage device. sto This refers to the rated capacity of the energy storage device. , , These are the standby efficiency, charging efficiency, and discharging efficiency of energy storage devices, respectively. and These represent the energy storage power and energy release power of the energy storage device at time t on a typical day s; (5d) The operating characteristics constraints of the power supply equipment are upper and lower limits of operating power: ； In the formula, P sup,max and P sup,min These represent the upper and lower limits of the power output of the energy supply equipment.

6. An electronic device, comprising: processor; as well as A memory storing computer program instructions, which, when executed by the processor, cause the processor to perform the stochastic planning method for an integrated energy system in an industrial park that takes into account source-load uncertainty, as described in any one of claims 1-5.

7. A computer-readable storage medium having stored thereon computer program instructions, which, when executed by a processor, cause the processor to perform a stochastic planning method for an integrated energy system in an industrial park that takes into account source-load uncertainty, as described in any one of claims 1-5.

Images