A comprehensive energy system distributed stochastic optimization method

Abstract

The application provides a kind of integrated energy system distributed random optimization method, belongs to integrated energy system optimization scheduling and uncertainty processing technical field, including: obtaining the historical data of multidimensional physical quantity in integrated energy system, and generating a plurality of typical scenes and its occurrence probability;For each subject, a two-stage random optimization model is constructed, the first stage decision variable is the point-to-point cross-zone transaction power between each subject, and the second stage decision variable is the power response and adjustment state of the internal equipment of each subject under each typical scene;A plurality of subjects are cooperatively optimized by using an alternating direction multiplier method, only the point-to-point transaction power information is exchanged between each subject, and the iteration is solved until the global clearing convergence condition is met, and the optimization scheduling scheme of each subject is output.The method quantifies uncertainty by scene generation, cooperatively anticipates decision-making by two-stage optimization, and protects privacy by ADMM distributed solution, to realize stable mutual aid operation of multiple subjects in uncertain environment.

Description

Technical Field

[0001] This invention belongs to the field of integrated energy system optimization scheduling and uncertainty handling technology, specifically relating to a distributed stochastic optimization method for integrated energy systems. Background Technology

[0002] Integrated energy systems (IES) with combined cooling, heating, and power (CCHP) as their core architecture have become an important carrier supporting the transformation of the energy structure due to their multi-energy cascade utilization and high energy efficiency. Multiple geographically adjacent microgrid entities achieve cross-entity source-load mutual assistance and multi-energy synergy through power interaction. However, the dense integration of a high proportion of distributed renewable energy sources and the spatiotemporal interweaving of diverse loads pose significant challenges to the system due to multi-dimensional and multi-grade source-load uncertainties. Traditional methods struggle to mitigate the risks of high-frequency, sudden fluctuations within the day while simultaneously protecting microgrid business privacy.

[0003] Existing short-term fine-grained scheduling and distributed collaborative methods have significant drawbacks: First, directly fitting historical environmental data ignores prior information such as weather forecasts, leading to amplified uncertainty, excessive redundancy, and a tendency to generate zero-value noise distortion of cold and hot loads during off-seasons, which violates physical norms. Second, traditional overall residual fitting severs the inherent strong spatiotemporal physical coupling and interdependence between wind, solar, temperature, and multiple loads, resulting in scenario distortion. Third, centralized optimization faces serious risks of commercial privacy leakage, while existing distributed optimization is mostly limited to deterministic frameworks, which can generate uncertain fluctuations when a large amount of renewable energy is connected, making it difficult for traditional integrated energy systems to operate stably. Summary of the Invention

[0004] In view of the above-mentioned defects or deficiencies in the prior art, the present invention provides a distributed stochastic optimization method for integrated energy systems, comprising: Acquire historical data on multidimensional physical quantities in integrated energy systems; Based on the historical data, multiple typical scenarios and the probability of occurrence of each typical scenario are generated; A two-stage stochastic optimization model is constructed for each entity. The first-stage decision variable of the two-stage stochastic optimization model is the point-to-point cross-regional electricity transaction volume between the entities, and the second-stage decision variable is the power response and regulation status of the internal equipment of each entity under various typical scenarios. The alternating direction multiplier method is used to perform distributed collaborative optimization on multiple entities. The entities only interact with each other through point-to-point electricity trading information. The solution is iteratively solved until the global clearing convergence condition is met, and the optimized scheduling scheme of each entity is output.

[0005] According to the technical solution provided by the present invention, the step of generating multiple typical scenarios and the probability of occurrence of each typical scenario based on the historical data includes: A dual-branch prediction network is constructed for the multidimensional physical quantity. The dual-branch prediction network includes a continuous value regression main network and a start-stop classification auxiliary network. The main network outputs continuous prediction values ​​of the multidimensional physical quantity, and the auxiliary network outputs the start-stop state probabilities of the cooling load and heating load and constructs a start-stop mask matrix. The continuous prediction values ​​are fused with the start-stop mask matrix to obtain the corrected baseline prediction values. Calculate the prediction residual between the true value of the multidimensional physical quantity and the baseline prediction value, and perform mode decomposition on the prediction residual to obtain multiple residual sub-modes in different frequency bands; Each residual sub-mode is fitted with a probability density, and random sampling is performed based on the fitted probability density model. The sampling results are superimposed on the baseline prediction value and the start / stop mask matrix is ​​applied to generate multiple initial prediction scenarios. High-dimensional joint feature flattening is performed on each of the initial prediction scenarios, and the flattened feature vectors are clustered and reduced to obtain a set of typical scenarios and the occurrence probability of each typical scenario.

[0006] According to the technical solution provided by the present invention, the continuous value regression main network is constructed based on a long short-term memory network, and the start-stop classification auxiliary network is constructed based on a long short-term memory network and a sigmoid activation function; the fusion process is to perform Hadamard product operation on the continuous predicted values ​​of cooling load and heating load with the corresponding start-stop state probability masks, respectively.

[0007] According to the technical solution provided by the present invention, the modal decomposition of the predicted residual includes: Variational mode decomposition is used to decompose the prediction residual into low-frequency mode, mid-frequency mode, high-frequency mode and residual component; the low-frequency mode represents the long-term trend error in the prediction deviation, the mid-frequency mode represents the intraday periodic fluctuation error in the prediction deviation, and the high-frequency mode represents the random mutation error in the prediction deviation.

[0008] According to the technical solution provided by the present invention, the probability density of each residual sub-mode is fitted by the kernel density estimation method, and the initial prediction scene is generated by the Monte Carlo sampling method; the high-dimensional joint feature flattening is to concatenate the data of all time steps and all physical dimensions contained in each initial prediction scene to form a single high-dimensional joint feature row vector, and the clustering reduction adopts the K-Medoids clustering algorithm.

[0009] According to the technical solution provided by the present invention, the construction of a two-stage stochastic optimization model for each subject includes: The optimization objective is to minimize the expected total cost across all scenarios throughout the day. The expected total cost includes the expected cost of scenario operation and the cost of point-to-point electricity trading. The expected cost of scenario operation includes the cost of purchasing electricity from the external grid, the cost of curtailment penalties, and the cost of natural gas fuel.

[0010] According to the technical solution provided by the present invention, the method of using alternating direction multipliers to perform distributed collaborative optimization of multiple entities includes: In each iteration of the alternating direction multiplier method, each entity sends its expected transaction power vector to the virtual trading center. The virtual trading center then aggregates the data and sends the sum of transaction power information of all entities except the entity itself to each entity. Each entity updates its augmented Lagrange objective function based on the received information and solves it again.

[0011] According to the technical solution provided by the present invention, the augmented Lagrange objective function includes a proximal damping term, which is used to limit the magnitude of change of the current round of trading decisions relative to the previous round of trading decisions.

[0012] According to the technical solution provided by the present invention, the method for determining the global clearing convergence condition is as follows: After collecting the trading decisions returned by all participants, the virtual trading center calculates the sum of the net trading gaps in the entire market for each time period as the trading residual, and determines whether the trading residual is less than a preset convergence threshold. If it is satisfied, it is determined that the multi-participant market has achieved distributed clearing, the algorithm terminates and outputs an optimized scheduling scheme; if it is not satisfied, the Lagrange multipliers are updated and the next round of iteration begins.

[0013] According to the technical solution provided by the present invention, the Lagrange multiplier is used as the market transaction electricity price, which is updated using the gradient descent method. An upper limit and a lower limit are set for the updated market transaction electricity price. The lower limit is 0, and the upper limit is the difference between the external grid purchase price and the preset price margin, so that the internal transaction price between entities is lower than the price of directly purchasing electricity from the external grid.

[0014] Compared with existing technologies, the beneficial effects of this invention are as follows: By performing predictive processing and residual decomposition to generate typical scenarios, it effectively avoids the problems of uncertainty amplification and non-operating season noise distortion caused by directly fitting environmental data in traditional methods. At the same time, it avoids the drawback of the overall residual fitting severing the coupling relationship of multiple physical quantities, ensuring the physical authenticity and fidelity of the scenario. Furthermore, by constructing a two-stage stochastic optimization model, it achieves the coordinated optimization of forward-looking planning and adaptive scheduling by making hierarchical decisions on cross-regional electricity trading and internal equipment response. At the same time, it adopts the alternating direction multiplier method for distributed collaborative optimization, with each entity only interacting with point-to-point electricity trading information. Under the premise of protecting the commercial privacy of each entity, iterative solution is used until global clearing convergence, effectively offsetting the multidimensional fluctuation risks brought about by the high proportion of renewable energy access, and realizing the stable and mutually supportive operation of multiple entities. Attached Figure Description

[0015] Other features, objects, and advantages of this application will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings: Figure 1 A flowchart illustrating the steps of the distributed stochastic optimization method for an integrated energy system provided in Example 1; Figure 2 This is a schematic diagram of the VMD decomposition of residuals in each dimension in Example 2; Figure 3 This is a schematic diagram of the probability density fitting for each dimension in Example 2; Figure 4 This is a schematic diagram illustrating the generation and clustering reduction results of the three scenarios in Example 2; Figure 5 This is a schematic diagram showing the generation and clustering reduction results of the other three scenarios in Example 2; Figure 6 This is the power dispatch balance diagram of main body 1 in Example 2; Figure 7 This is a cold energy scheduling balance diagram for main body 1 in Example 2; Figure 8 This is the thermal energy scheduling balance diagram of main body 1 in Example 2; Figure 9 This is a diagram showing the electrical energy interaction power between the three entities in Example 2. Detailed Implementation

[0016] The present application will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, only the parts relevant to the invention are shown in the accompanying drawings.

[0017] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0018] Example 1 As mentioned in the background section, this embodiment proposes a distributed stochastic optimization method for integrated energy systems to solve the optimization scheduling problem of multi-entity integrated energy systems under the impact of source-load uncertainty, while taking into account operational economy, physical authenticity and the protection of entity privacy. Figure 1 Please refer to the flowchart of the distributed stochastic optimization method for integrated energy systems provided by this invention. Figure 1 The method includes the following steps S1-S4: S1: Obtain historical data on multidimensional physical quantities in the integrated energy system.

[0019] Specifically, in step S1, the integrated energy system encompasses the coupled conversion and coordinated supply of multiple energy forms, including cooling, heating, and electricity. To accurately characterize the supply and demand characteristics of various energy sources within the system and their temporal variations, historical operational data of multidimensional physical quantities is needed as the basis for subsequent prediction and optimization. These multidimensional physical quantities include wind speed, temperature, solar radiation, cooling load, heating load, and electrical load. Wind speed and solar radiation are key environmental factors affecting the output of wind and solar power generation, while temperature is closely related to the demand levels of cooling and heating loads. Cooling load, heating load, and electrical load reflect the system's diverse energy needs. This historical data is typically collected at a certain time resolution, such as hourly intervals, covering multiple complete operating cycles, such as seasons or years, to fully reflect the temporal characteristics, periodic patterns, and extreme operating conditions of the data, providing data support for subsequent scenario generation and optimization decisions. By collecting historical data of multidimensional physical quantities, it is possible to comprehensively capture uncertainties on both the source and load sides, avoiding distortion in subsequent scenario generation due to missing information, thereby ensuring that the optimized scheduling scheme has sufficient adaptability and robustness in actual operation.

[0020] S2: Generate multiple typical scenarios and the probability of occurrence of each typical scenario based on the historical data.

[0021] Specifically, in step S2, in actual operation, renewable energy output and diverse load demands are comprehensively affected by multiple factors such as weather conditions and human activities, exhibiting significant randomness and volatility. If scheduling decisions are made based solely on a single deterministic forecast value, the system will face the risk of supply-demand imbalance if the actual source load deviates from the forecast value. To quantify the impact of this uncertainty on system operation, this step generates multiple representative typical scenarios based on the historical data obtained in step S1, and assigns a corresponding probability of occurrence to each scenario. Typical scenarios are probabilistic descriptions of possible future source load state combinations, and each scenario covers source load information at various times within the future scheduling cycle. By introducing typical scenarios and their occurrence probabilities, the uncertainty that was originally difficult to describe precisely is transformed into a finite number of discrete probabilistic expressions, enabling the subsequent optimization scheduling problem to be modeled and solved within a probabilistic framework. The specific generation method of typical scenarios will be discussed below. Figures 2-5 Please provide a detailed explanation.

[0022] S3: Construct a two-stage stochastic optimization model for each entity. The first-stage decision variable of the two-stage stochastic optimization model is the point-to-point cross-regional electricity transaction between the entities, and the second-stage decision variable is the power response and regulation status of the internal equipment of each entity under various typical scenarios.

[0023] Specifically, in step S3, under the multi-entity architecture of the integrated energy system, multiple geographically adjacent microgrid entities achieve cross-entity source-load mutual assistance and multi-energy coordination through point-to-point electricity trading. However, the renewable energy output and load demand within each entity are affected by uncertainty, which needs to be considered in the optimization decision-making. To this end, this step constructs a two-stage stochastic optimization model for each entity.

[0024] The core idea of ​​the two-stage stochastic optimization model is to divide the decision-making process into two stages to address information uncertainty at different time scales. The first stage of decision-making occurs before the actual uncertainty scenario materializes. The decision variables are point-to-point inter-regional electricity transactions between entities, including the purchase and sale of electricity with other entities. This part of the decision is forward-looking, needing to be determined before the uncertainty is fully revealed, and therefore does not change with specific scenarios; it belongs to the "here-now" decision-making. The second stage of decision-making occurs after the uncertainty scenario materializes. The decision variables are the power response and regulation status of the internal equipment of each entity under typical scenarios, including the start-stop and output of internal combustion engines, the charging and discharging power of energy storage systems, the heat generation power of gas boilers, and the cooling power of absorption chillers and electric chillers. This part of the decision can be adaptively adjusted according to the actual scenario, belonging to the "wait-and-see" decision-making. Through this hierarchical decision-making structure, the two-stage stochastic optimization model can achieve synergistic optimization of forward-looking transaction plans and adaptive scheduling strategies under uncertain environments. It locks in cross-regional mutual assistance solutions in advance through the first-stage decision to ensure synergy among multiple entities, and empowers the system to flexibly adjust its internal output according to the actual scenario through the second-stage decision, thereby effectively offsetting the impact of source-load uncertainty on system operation.

[0025] S4: The alternating direction multiplier method is used to perform distributed collaborative optimization on multiple entities. The entities only interact with each other to exchange electricity information point by point. The solution is iteratively solved until the global clearing convergence condition is met, and the optimized scheduling scheme of each entity is output.

[0026] Specifically, in step S4, a two-stage stochastic optimization model was constructed for each entity in step S3. However, the point-to-point cross-regional electricity trading between multiple entities is interdependent: the electricity sales decision of one entity necessarily corresponds to the electricity purchase decision of one or more other entities, and globally, the total electricity trading volume of all entities needs to be balanced. If a centralized optimization method is adopted, the privacy information such as equipment parameters and load data of all entities needs to be aggregated to the central computing node. This not only faces serious risks of commercial privacy leakage, but also causes the curse of dimensionality when the number of entities is large, leading to a sharp decline in computational efficiency. To solve the above problems, this step adopts the alternating direction multiplier method as the top-level distributed collaborative optimization framework.

[0027] The alternating direction multiplier method is a decomposition and coordination algorithm suitable for large-scale distributed optimization problems. In this invention, each entity only needs to exchange point-to-point electricity transaction information during the iteration process, without disclosing internal equipment capacity, load data, or information on generated uncertain scenarios. Through information aggregation and feedback from the virtual trading center, each entity independently solves its own local optimization problem and updates its trading decisions. Through multiple iterations, the trading decisions of each entity gradually converge to global consistency, ultimately achieving distributed clearing in a multi-entity market. This interaction mechanism, while protecting the commercial privacy of each entity, decomposes the originally high-dimensional centralized optimization problem into multiple low-dimensional local subproblems that are solved in parallel, significantly reducing the computational complexity of a single optimization.

[0028] In a preferred embodiment, step S2 generates multiple typical scenarios and the probability of occurrence of each typical scenario based on the historical data, including the following steps S21-S24: S21: Construct a dual-branch prediction network for the multidimensional physical quantity. The dual-branch prediction network includes a continuous value regression main network and a start-stop classification auxiliary network. The main network outputs continuous prediction values ​​of the multidimensional physical quantity, and the auxiliary network outputs the start-stop state probabilities of the cooling load and heating load and constructs a start-stop mask matrix. The continuous prediction values ​​are fused with the start-stop mask matrix to obtain the corrected baseline prediction values.

[0029] Specifically, in the scenario generation process, it is first necessary to obtain the basic predicted values ​​of the source load state at each moment within the future scheduling cycle, as a benchmark for subsequent uncertain scenario construction. However, the prediction difficulty and characteristics of different variables among multidimensional physical quantities vary significantly. The values ​​of wind speed, temperature, solar radiation, and electrical load have physical meaning at every moment throughout the year, and their prediction belongs to the continuous value regression problem. In contrast, cooling load and heating load have obvious seasonal characteristics. In non-operating seasons, such as the summer heating season or the winter cooling season, their true values ​​should be zero, rather than a very small positive number. If a continuous value regression model is directly used for prediction, the model will still output non-zero predicted values ​​in non-operating seasons, producing zero-value noise distortion that violates physical common sense. This noise, after being amplified by uncertainty sampling in subsequent scenario generation, will lead to excessive reserve capacity configuration and reduced economic efficiency.

[0030] To address these issues, this step constructs a dual-branch prediction network. The dual-branch prediction network comprises two branches: a continuous value regression main network and a start-stop classification auxiliary network. The continuous value regression main network is responsible for outputting preliminary continuous predictions for all multidimensional physical quantities, including wind speed, temperature, solar radiation, cooling load, heating load, and electrical load at various times within the future scheduling cycle. This main network learns the temporal dependencies in historical data to capture the evolution patterns of each physical quantity over time. Simultaneously, the start-stop classification auxiliary network specifically performs probabilistic predictions of the start-stop states of cooling and heating loads, outputting the start-stop probabilities for each load at various future times. The start-stop probability reflects the likelihood that the corresponding load is in operation at that moment.

[0031] A start-stop mask matrix is ​​constructed based on the start-stop state probabilities output by the auxiliary network. This start-stop mask matrix is ​​a binary matrix with the same dimension as the predicted cooling and heating loads. The elements within it are determined by the relationship between the start-stop state probabilities and a preset threshold, taking values ​​of 1 or 0. When the start-stop state probability is greater than or equal to the preset threshold, the mask value is 1, indicating that the corresponding load is in operation at that moment, and its predicted value should be retained. When the start-stop state probability is less than the preset threshold, the mask value is 0, indicating that the corresponding load is out of operation at that moment, and its predicted value should be masked. In this way, the start-stop mask matrix acts as a spatial filter, capable of identifying and eliminating predicted values ​​that do not conform to physical seasonal patterns.

[0032] After obtaining the continuous predicted values ​​and the start / stop mask matrix, the two are fused. Specifically, the continuous predicted values ​​of cooling load and heating load are multiplied by the corresponding start / stop state masks, i.e., the predicted value at the corresponding position is multiplied by the mask value one by one. Other multidimensional physical quantities besides cooling load and heating load, such as wind speed, temperature, solar radiation, and electrical load, do not need to be masked; their continuous predicted values ​​are directly used as the prediction results for the corresponding physical quantities. After the above fusion process, the noise floor of all non-operating seasons in the predicted values ​​of cooling load and heating load is forcibly set to zero, resulting in the corrected baseline predicted value. This baseline predicted value contains the prediction results of all physical dimensions at each time point in the future scheduling cycle, serving as the data benchmark for subsequent scenario generation. Through the construction and fusion processing of the above dual-branch prediction network, the negative impact of non-operating season noise distortion on subsequent scenario generation and optimized scheduling in traditional direct prediction methods is effectively eliminated, avoiding the problem of excessive reserve capacity redundancy caused by this. At the same time, the prediction results strictly conform to the physical operating laws of cooling and heating loads, improving the physical authenticity and reliability of the prediction data.

[0033] S22: Calculate the prediction residual between the true value of the multidimensional physical quantity and the benchmark prediction value, and perform mode decomposition on the prediction residual to obtain multiple residual sub-modes in different frequency bands.

[0034] Specifically, after obtaining the corrected baseline prediction value in step S21, it is compared with the actual values ​​of the multidimensional physical quantities collected during actual operation. The actual values ​​are derived from historical operation datasets and have the same time resolution and time range as the baseline prediction value. The difference between the actual value and the baseline prediction value for each physical dimension at each moment is calculated; this difference is the prediction residual. :

[0035] In the formula, For the true value, The baseline predicted value is used as the reference. The prediction residual reflects the degree and direction of deviation of the baseline predicted value from the actual value, and is a core element in characterizing prediction uncertainty in subsequent scene generation.

[0036] The prediction residuals contain a mixture of error components from various sources and with different characteristics. To accurately characterize the uncertainty and perform precise sampling, this step performs mode decomposition on the prediction residuals. Mode decomposition separates the original residual signal into several independent sub-mode components based on their different frequency characteristics, each representing an error component in a specific frequency band. By performing mode decomposition on the prediction residuals, the originally mixed multi-source error signals are decoupled into independent sub-modes. This allows for subsequent probability density fitting and random sampling of errors in different frequency bands, avoiding the distribution distortion problem caused by traditional overall residual fitting, which mixes errors from different sources and with different characteristics.

[0037] S23: Perform probability density fitting on each of the residual sub-modes, perform random sampling based on the fitted probability density model, superimpose the sampling results onto the baseline prediction value and apply the start / stop mask matrix to generate multiple initial prediction scenarios.

[0038] Specifically, after obtaining multiple residual sub-modes in different frequency bands in step S22, each residual sub-mode has its own statistical distribution characteristics. To generate random samples that conform to the original error distribution, this step first performs probability density fitting on each residual sub-mode. The purpose of probability density fitting is to establish a mathematical model, i.e., a probability density function, for each residual sub-mode that accurately describes its value distribution. This probability density function quantifies the likelihood of each residual sub-mode occurring within different value intervals.

[0039] After completing the probability density fitting, independent random sampling is performed based on the probability density model corresponding to each residual sub-mode. For each time step within a future scheduling period, a random sample is drawn from the probability density function of each residual sub-mode in each dimension. All residual sub-mode samples drawn at the same time step are accumulated according to their corresponding mode order to obtain a comprehensive error sequence. This comprehensive error sequence represents the overall deviation of the predicted value from the true value under this sampling.

[0040] The comprehensive error sequence obtained from sampling is superimposed onto the baseline prediction curve output in step S21, and the start-stop mask matrix is ​​applied again to the cooling load and heating load. This superposition operation introduces random errors into the baseline prediction value, ensuring that the generated scenario reflects various potential source load deviations in the future. The start-stop mask matrix is ​​applied again to the cooling and heating loads after superposition to ensure that their scenario data is always forced to zero during the off-season, preventing the reappearance of non-zero noise floor due to error superposition. Through these operations, each set of random sampling results generates an initial prediction scenario. The number of initial scenarios generated is set to a preset first number, such as 1000. Each initial prediction scenario contains source load data for all physical dimensions at each moment within a future scheduling cycle, such as 24 hours. By combining probability density fitting and random sampling, the generated initial prediction scenario fully reflects the distribution characteristics and randomness of errors in each frequency band, while preserving the statistical correlation information between different physical dimensions.

[0041] S24: Flatten the high-dimensional joint features of each of the initial prediction scenarios, and then perform clustering reduction on the flattened feature vectors to obtain a set of typical scenarios and the probability of occurrence of each typical scenario.

[0042] Specifically, step S23 generates a large number of initial prediction scenarios. While these scenarios statistically reflect the uncertainty of the source load, directly substituting all initial scenarios into the subsequent optimization model would lead to a sharp increase in the size of the optimization problem and an excessive computational burden. Therefore, it is necessary to reduce the number of initial scenarios, replacing the massive number with fewer typical scenarios, while preserving as much of the uncertainty information contained in the initial scenarios as possible.

[0043] However, a significant challenge arises in reducing the number of physical quantities in a scenario: the six physical dimensions—wind speed, temperature, solar radiation, cooling load, heating load, and electrical load—naturally exhibit spatiotemporal physical coupling relationships. For example, solar radiation is positively correlated with temperature, temperature is negatively correlated with cooling and heating loads, and wind speed is positively correlated with wind power output. Performing independent dimensionality reduction or clustering on each physical dimension artificially severs these inherent physical coupling relationships, causing the reduced typical scenario to lose its physical realism and fail to accurately reflect the interconnected changes between physical quantities in the actual system.

[0044] To address this issue, this step performs high-dimensional joint feature flattening on each initial prediction scene. This involves concatenating the data from all time steps and all physical dimensions within each scene to form a high-dimensional joint feature row vector. Each row vector fully preserves the values ​​of all physical dimensions and their inherent relationships at all times within the corresponding scene. A clustering algorithm is then used to reduce the number of high-dimensional joint feature row vectors, dividing the initial scenes into a predetermined number of categories. Cluster centers are extracted from each category as typical scenes, and the proportion of initial scenes in each category relative to the total number of scenes is calculated as the probability of that typical scene occurring.

[0045] In a preferred embodiment, in step S21, the continuous value regression main network is constructed based on a long short-term memory network, and the start-stop classification auxiliary network is constructed based on a long short-term memory network and a sigmoid activation function; the fusion process involves performing Hadamard product operations on the continuous predicted values ​​of cooling load and heating load with the corresponding start-stop state probability masks, respectively.

[0046] Specifically, in step S21, when constructing the dual-branch prediction network, the continuous value regression main network is built based on a Long Short-Term Memory (LSTM) network. It takes multi-dimensional physical quantity data from historical moments as input, extracts temporal features through the LSTM network layers, and finally outputs preliminary continuous predictions of wind speed, temperature, solar radiation, cooling load, heating load, and electrical load for future moments through a fully connected layer. These predictions are denoted as... During training, the network aims to minimize the mean squared error between the predicted and actual values, and continuously adjusts the network weight parameters through backpropagation.

[0047] The start-stop classification auxiliary network is also built on a Long Short-Term Memory (LSTM) network, with a Sigmoid activation function added after the output layer. Unlike the continuous value regression main network, the start-stop classification auxiliary network's task is to perform binary classification probability prediction of the start-stop status of cold and hot loads. This network takes multi-dimensional physical quantity data from historical moments as input, and the output layer uses a Sigmoid activation function to map the network output to probability values ​​between 0 and 1. The closer the output value is to 1, the higher the probability that the corresponding load is in an operating state at that moment; the closer it is to 0, the higher the probability that it is in a shutdown state. During training, the network optimizes by minimizing the cross-entropy loss between the predicted probability and the true start-stop label. After training, it can output the start-stop status probabilities of cold and hot loads, constructing a start-stop mask matrix. :

[0048] In the formula, This represents the predicted probability value for the classification state.

[0049] In obtaining continuous predicted values and start / stop mask matrix Next, the two are fused. Specifically, the fusion process involves performing a Hadamard product operation on the continuous predicted values ​​of cooling and heating loads with their corresponding start / stop mask vectors to obtain the corrected baseline predicted values. :

[0050] Start / Stop Mask Matrix When the value is 1, the predicted value at that moment is retained; when the value is 0, the predicted value at that moment is forcibly set to zero. Through this calculation, the prediction noise of cooling and heating loads during the non-operating season is completely filtered out at once, while the effective prediction values ​​for the operating season are fully retained, achieving strict consistency between the prediction results and physical laws. Other multidimensional physical quantities do not require masking processing; their continuous prediction values ​​are directly used as the prediction results for the corresponding dimensions.

[0051] In a preferred embodiment, step S22 involves performing mode decomposition on the predicted residuals, including: Variational mode decomposition is used to decompose the prediction residual into low-frequency mode, mid-frequency mode, high-frequency mode and residual component; the low-frequency mode represents the long-term trend error in the prediction deviation, the mid-frequency mode represents the intraday periodic fluctuation error in the prediction deviation, and the high-frequency mode represents the random mutation error in the prediction deviation.

[0052] Specifically, for the predicted residuals When performing mode decomposition, variational mode decomposition is used to adaptively decompose the predicted residual signal into a predetermined number of eigenmode functions with finite bandwidth. Variational mode decomposition constructs and solves a variational problem to separate the original signal into multiple independent sub-mode components based on different frequency characteristics. Each sub-mode is distinguishable from the others in the frequency domain, so that error components in different frequency bands can be effectively decoupled.

[0053] Predicting residuals Decomposed into 3 intrinsic mode functions ( ) and a residual component ( ),Right now:

[0054] In the formula, For the first One eigenmode function .

[0055] The three intrinsic mode functions (EMFs) correspond to error components in three different frequency bands: low-frequency mode, mid-frequency mode, and high-frequency mode. The low-frequency mode is the first EMF (EMF). The mid-frequency mode (MFM) characterizes the slowly changing long-term trend error in prediction bias, reflecting systematic prediction shifts caused by factors such as equipment performance degradation and seasonal environmental changes. This error persists over a long timescale and changes gradually. The mid-frequency mode is the second intrinsic mode function (IMF). This characterizes the intraday periodic fluctuations in prediction bias related to human daily activity patterns or the diurnal cycle of the environment. Examples include differences in energy consumption patterns between weekdays and rest days, and load fluctuations caused by diurnal temperature variations. These errors exhibit an approximately regular oscillation characteristic with a 24-hour period. The high-frequency mode is the third intrinsic mode function (...). This characterizes the random abrupt error in prediction bias caused by unforeseen factors such as sudden gusts, brief cloud cover, or random equipment start-ups and shutdowns. This type of error is characterized by random occurrence, drastic amplitude changes, and short duration. Residual component. This represents the remaining terms in the residuals that were not captured by the three eigenmode functions mentioned above.

[0056] Through the variational mode decomposition described above, the originally superimposed multi-source error signals are decoupled into three independent sub-modes with clear frequency bands and distinct physical meanings, and a residual component. This allows for subsequent probability density fitting and random sampling of errors in different frequency bands, avoiding the distribution distortion problem caused by the traditional overall residual fitting that mixes errors from different sources and with different characteristics.

[0057] In a preferred embodiment, in step S23, the kernel density estimation method is used to fit the probability density of each residual sub-mode, and the Monte Carlo sampling method is used to generate the initial prediction scene; in step S24, the high-dimensional joint feature flattening is to concatenate the data of all time steps and all physical dimensions contained in each initial prediction scene to form a single high-dimensional joint feature row vector, and the clustering reduction adopts the K-Medoids clustering algorithm.

[0058] Specifically, kernel density estimation is a nonparametric probability density estimation method that does not require a predefined data distribution and can adaptively fit a probability density function of arbitrary shape based solely on the sample data itself. Specifically, a standard Gaussian kernel function is used as the kernel function to perform kernel density estimation on the sample data of each residual submode, establishing a probability density function model (PDF) for the error in each frequency band.

[0059] Based on the fitted probability density functions, an initial prediction scenario is generated using Monte Carlo sampling. The number of initial scenarios is set to 1000. For a 24-hour prediction window, at each time point, a random sample is independently drawn from each of the four modal probability density function models across all dimensions. All residual submodal samples drawn at the same time point are summed to obtain the comprehensive error sequence for that time point. After iterating through all 24 time points, a complete comprehensive error sequence is obtained, covering all time points within the entire prediction window. The obtained comprehensive error sequence is then superimposed on the baseline prediction value output in step S21. Then, apply the start / stop status mask again:

[0060] In the formula, To represent the physical dimension (such as wind, light, cold, heat, electricity, etc.) generating the first The values ​​of the initial scene, For the output regarding physical dimensions The baseline prediction curve value, This is a comprehensive error sequence obtained by independently and randomly sampling from the probability density function models of each mode for this physical dimension and then summing the results. For physical dimension The corresponding start / stop status mask variable.

[0061] Each complete sampling and overlay process generates a massive initial prediction scene containing all 6 physical dimensions and covering 24 time points. The above process is repeated until a preset 1000 initial prediction scenes are generated, resulting in a 6-dimensional massive initial prediction scene matrix.

[0062] In step S24, the high-dimensional joint feature flattening is implemented as follows: All 24 time steps and all 6 physical dimensions of the data in each initial prediction scene are concatenated end-to-end to form a single high-dimensional joint feature row vector. After flattening, each scene forms a high-dimensional joint feature row vector with dimensions 24 x 6, i.e., 144. Through this flattening operation, each row vector fully preserves the values ​​of all physical dimensions in the corresponding scene at all times and their inherent relationships, ensuring that the coupling relationship between physical dimensions is not disrupted during subsequent clustering.

[0063] After flattening the high-dimensional joint features, the K-Medoids clustering algorithm is used to cluster and reduce all high-dimensional joint feature row vectors. The K-Medoids algorithm is a centroid-based clustering algorithm where cluster centers must be selected from the sample data rather than calculated from the mean. During clustering, the Manhattan distance is used as a measure of similarity between samples, measuring the degree of difference between row vectors corresponding to different initial scenes. The K-Medoids algorithm divides all high-dimensional joint feature row vectors into a predetermined number of categories, iteratively searching for representative cluster centers that minimize the sum of distances between samples and centroids within each category. This ensures that scenes within the same category are as similar as possible in terms of spatiotemporal joint features, while scenes between different categories are as dissimilar as possible. Finally, the cluster center of each category is extracted as a typical scene. The sample size percentage of each category is then used as the probability of occurrence of this typical scenario. .

[0064] In a preferred embodiment, step S3 involves constructing a two-stage stochastic optimization model for each subject, including: The optimization objective is to minimize the expected total cost across all scenarios throughout the day. The expected total cost includes the expected cost of scenario operation and the cost of point-to-point electricity trading. The expected cost of scenario operation includes the cost of purchasing electricity from the external grid, the cost of curtailment penalties, and the cost of natural gas fuel.

[0065] Specifically, in step S3, when constructing a two-stage stochastic optimization model for each subject, the optimization objective is to minimize the expected total cost across all scenarios throughout the day. Since the probability of occurrence of each typical scenario is different, the contribution of its corresponding operating cost to the total cost should also be weighted according to its probability.

[0066] Expected total cost This is the weighted sum of the expected operating costs for each scenario and the costs of peer-to-peer electricity trading, i.e.:

[0067] in, This represents the total number of typical scenarios. For the scene The probability of occurrence The cost of point-to-point electricity transactions reflects the net expenditure generated by the buying and selling of electricity between entities through cross-regional transactions. Positive values ​​represent electricity purchase expenditures, while negative values ​​represent electricity sales revenues. This cost is calculated based on the amount of electricity traded between the entities and the transaction price formed through negotiation or market clearing. The expected cost of running the scenario. The cost of purchasing electricity from the external power grid reflects the economic cost that the entity needs to pay when its own power generation is insufficient. Its size is related to the amount of electricity purchased and the time-of-use price of the external power grid. The cost of purchasing electricity is higher during peak hours and lower during off-peak hours. The introduction of this cost makes the optimization model prioritize the use of internal power generation resources or increase energy storage discharge to reduce the purchase of electricity when the electricity price is high. The curtailment penalty cost reflects the economic losses caused by the forced curtailment of wind and solar power due to the inability of renewable energy output to be fully absorbed. The introduction of this cost means that when faced with a surplus of renewable energy output, the optimization model will prioritize improving the renewable energy absorption rate through methods such as energy storage charging and increasing the electricity consumption of electric chillers, so as to reduce the curtailment penalty. Natural gas fuel cost reflects the fuel expenses incurred by internal combustion engines and gas boilers consuming natural gas. This cost is an important indicator for measuring the economic efficiency of multi-energy conversion equipment.

[0068] In a preferred embodiment, the two-stage distributed stochastic optimization in step S3 must strictly adhere to the following physical boundaries and coupling constraints during the solution process: I. Market Transaction Auxiliary Constraints: To address the issue of network access fee calculation, a linearization transformation constraint for the absolute value of transaction power is introduced:

[0069] In the formula, For actual network transaction power, It is an auxiliary variable for the absolute value of trading power.

[0070] II. Physical boundary constraints for equipment operation: The equations of state for the energy storage system and thermal storage tank limit the upper and lower limits of charging and discharging power, and ensure the conservation balance of the charged / thermal state throughout the day:

[0071] In the formula, In the scene Underground energy storage system or thermal storage tank The heat storage state during a certain period of time, , In the scene Underground energy storage system The actual charging / discharging (or heating / heat dissipation) power during the time period. , This indicates the charging / discharging (or energy charging / discharging) efficiency of the energy storage system. This is the rated maximum capacity of the energy storage system.

[0072] Internal combustion engine gradeability and output constraints, introducing dual-phase start-stop state variables. Limit its minimum safe output and maximum permissible gradeability:

[0073] In the formula, In the scene The internal combustion engine in the current The actual active power output during the time period. For internal combustion engines The binary start / stop status variable for the time period (a value of 1 represents startup and a value of 0 represents shutdown). , This indicates the upper limit of the maximum rated power and the lower limit of the minimum safe output of an internal combustion engine. This is the maximum allowable power ramp-up limit (i.e., the maximum allowable power increase or decrease) for an internal combustion engine between adjacent time steps.

[0074] III. Multi-energy balance coupling constraint Cooling balance: The sum of the cooling output of the absorption chiller and the electric chiller must equal the predicted cooling load of the scenario. .

[0075] Thermal energy balance: After deducting the heat consumption for charging and cooling from the heat generated by the gas boiler, the waste heat recovery from the internal combustion engine, and the heat release from the thermal storage tank, the heat load of the predicted scenario must be met. need.

[0076] Power balance (core coupling constraint):

[0077] In the formula, For the scene Down The actual active power output of the wind turbine during the time period For the scene Down The actual active power output of photovoltaic power generation during the time period. For the scene Down The actual electrical power generated by the internal combustion engine during the specified time period. For the scene Down The actual discharge power of the time-limited energy storage system For the scene Down The electrical power purchased by the time-of-use system from the external power grid. In predicting scenarios Down The system's own standard basic electrical load requirements during the time period Indicates in the scene Down The electrical power consumed by the electric refrigeration unit during a given period. For the scene Down The actual charging power of the time-of-use energy storage system For the scene Down The net transaction power of the system and other microgrid entities during the period of cross-regional point-to-point interaction (Note: it has no net transaction power). (Subscript, because it is a forward-looking decision made in the first phase and does not change with uncertain scenarios). For the scene Down The electrical power that the system sells (or abandons) to the external power grid when it has a surplus of electricity during a certain period.

[0078] This equation realizes the first-stage trading variables. Decoupling and coordinated balance adjustment of various scenario-adaptive variables (actual wind and solar power output, inter-grid interaction power, load demand, etc.) in the second phase in the dimension of electricity.

[0079] In a preferred embodiment, step S4 employs the alternating direction multiplier method to perform distributed collaborative optimization on multiple agents, including: In each iteration of the alternating direction multiplier method, each entity sends its expected transaction power vector to the virtual trading center. The virtual trading center then aggregates the data and sends the sum of transaction power information of all entities except the entity itself to each entity. Each entity updates its augmented Lagrange objective function based on the received information and solves it again.

[0080] Specifically, in step S4, the alternating direction multiplier method is used to perform distributed collaborative optimization on multiple entities. The core of this method is to gradually bring the transaction decisions of each entity towards global consistency through iterative interaction. During the iterative process of multi-entity distributed optimization, the local optimization objective of each entity cannot only consider its own economic efficiency but also needs to introduce global transaction balance constraints to promote the cross-regional electricity transaction between multiple entities towards global supply and demand balance. At the same time, since multiple entities make independent decisions, without restrictions on the magnitude of decision changes, each entity may experience drastic strategy oscillations due to game theory behavior during iteration, making it difficult for the algorithm to converge. Therefore, each entity... In the The local optimization objective function in the next iteration changes from the original expected total cost. Extended to augmented Lagrangian objective function :

[0081] In the formula, This is the secondary penalty coefficient. In the previous iteration, excluding the main body The sum of the transaction power of all other entities; For proximal damping parameters, This is the proximal damping term.

[0082] The augmented Lagrange objective function adds two terms to the original cost function: the first term is a quadratic penalty term, which introduces a quadratic penalty parameter. This term acts on the global transaction balance constraint, and its mathematical meaning is a quadratic penalty on the degree to which the transaction decisions of each subject in the current round deviate from the global balance, which is used to accelerate the convergence speed of the algorithm to the global supply and demand balance; the second term is a proximal damping term, which introduces a proximal damping parameter and constructs the square term of the difference between the transaction decision in the current round and the transaction decision in the previous round. This term is used to limit the magnitude of the deviation of the current round decision from the previous round decision and avoid the drastic oscillation of strategy caused by multi-subject game.

[0083] Before the iteration begins, each entity has completed its two-stage stochastic optimization model initialization based on the current market electricity price information. In each iteration of the alternating direction multiplier method, each entity first sets its expected transaction electricity vector for the current round. Send to the virtual trading center. The virtual trading center will collect the responses from all entities on the [number]th [day]. After each round of trading decision vectors, global trading information is aggregated and calculated. Specifically, for each entity, the virtual trading center calculates the sum of the trading power of all other entities in the current round. And then send the information to the corresponding entity.

[0084] After receiving the sum of the trading power of all other entities from the virtual trading center, each entity substitutes this information into its own local optimization problem, updating the quadratic penalty term and the proximal damping term in the augmented Lagrange objective function to values ​​based on the current round information. Each entity then utilizes its local computing resources, combining the updated penalty term with the current market electricity price information, to invoke a two-stage stochastic optimization solver in parallel, obtaining a solution that satisfies the current augmented Lagrange objective function. Minimize the optimal physical operation strategy for this round and the updated trading decision This process continues until the next iteration begins. Throughout this process, each entity only sends the expected transaction volume vector to the virtual trading center, without disclosing its internal equipment capacity, load data, or information about generated uncertain scenarios. This ensures distributed collaborative optimization while protecting the commercial privacy of each entity.

[0085] In a preferred embodiment, the method for determining the global clearing convergence condition in step S4 is as follows: After collecting the trading decisions returned by all participants, the virtual trading center calculates the sum of the net trading gaps in the entire market for each time period as the trading residual, and determines whether the trading residual is less than a preset convergence threshold. If it is satisfied, it is determined that the multi-participant market has achieved distributed clearing, the algorithm terminates and outputs an optimized scheduling scheme; if it is not satisfied, the Lagrange multipliers are updated and the next round of iteration begins.

[0086] Specifically, during the iterative process of the alternating direction multiplier method in step S4, after collecting the transaction decisions returned by all entities in the k-th round, the virtual trading center needs to determine whether the transaction decisions of each entity in the current round have achieved global consistency. The transaction decisions proposed by each entity after independent optimization may not be perfectly matched globally; that is, there may be a discrepancy in the total amount of electricity planned to be sold by some entities and the total amount planned to be purchased by others. To quantify the degree of this global supply-demand imbalance, this step first calculates the sum of the net transaction gap in the entire market at each time period as the transaction residual. :

[0087] Transaction residual is an overall measure of the difference between planned electricity sales and purchases by all entities at all times. It reflects the degree of global supply and demand imbalance in the electricity traded among multiple entities in the current round and is a core quantitative indicator for judging whether the market has achieved clearing.

[0088] After calculating the transaction residual, compare the transaction residual with the preset convergence threshold. The comparison is performed. The preset convergence threshold is a very small positive number, for example, it can be set to 10. -3 The value of this threshold determines the precision requirement for algorithm termination. When the transaction residual is less than the preset convergence threshold, it indicates that the transaction volume among the entities has achieved global supply and demand balance within an acceptable error range. The transaction scheme of the current round can be regarded as a distributed clearing result jointly recognized by all entities. At this time, it is determined that the multi-entity market has achieved distributed clearing, the algorithm terminates the iteration, and outputs the optimized scheduling scheme of each entity, including the cross-regional transaction volume decision in the first stage and the internal equipment response strategy in the second stage.

[0089] If the transaction residual is greater than or equal to the preset convergence threshold, it indicates that there is still a significant global supply-demand imbalance in the current round of market transactions, and the participants have not yet reached a consensus on transactions. In this case, the Lagrange multipliers for this round need to be updated, and the updated multiplier information will be used in the next iteration. The specific implementation of updating the Lagrange multipliers will be explained below in conjunction with specific implementation methods. After the multiplier update is completed, the next iteration begins, and the above process is repeated until the convergence condition is met. Through the closed-loop mechanism of residual calculation and threshold determination, the algorithm can adaptively determine the termination time while ensuring solution accuracy, avoiding the problems of premature process termination or wasted computational resources that may result from a fixed number of iterations.

[0090] In a preferred embodiment, the Lagrange multiplier serves as the market transaction electricity price, which is updated using a gradient descent method. An upper limit and a lower limit are set for the updated market transaction electricity price, with the lower limit being 0 and the upper limit being the difference between the external grid purchase price and a preset price margin, so that the internal transaction price between entities is lower than the price of directly purchasing electricity from the external grid.

[0091] In this invention, the Lagrange multiplier is given the actual physical meaning of the market-traded electricity price, serving as the core price signal for the virtual trading center to guide various participants in adjusting their trading strategies. Specifically, the virtual trading center uses a gradient descent method to analyze the Lagrange multiplier (i.e., the market-traded electricity price). Perform adaptive updates:

[0092] In the formula, Update the price step size. For large power grids Time-of-use electricity pricing for different time periods This is a very small margin for price difference (e.g., 0.01 yuan).

[0093] The gradient descent method uses the transaction residual of the current round as the gradient direction, adjusting the Lagrange multipliers in the direction that reduces the global supply-demand imbalance. The price update step size is used to control the magnitude of each electricity price adjustment: a larger step size results in faster convergence but may cause oscillations, while a smaller step size leads to a smoother convergence but slower speed. In practical applications, the step size can be flexibly set according to the system size and the number of participants.

[0094] After each update, to ensure the physical and economic rationality of market-traded electricity prices, it is necessary to set upper and lower limits for the updated market-traded electricity prices. The lower limit is set to 0; the physical meaning of this setting is to ensure that market electricity prices are always non-negative, avoiding the illogical situation of negative electricity prices. The upper limit is set to... The physical and economic significance of this setting lies in ensuring that the internal point-to-point electricity transaction price between entities is always lower than the price of directly purchasing electricity from the external power grid. Because the internal transaction price always has a price advantage, driven by economic considerations, entities will prioritize internal power exchange rather than directly purchasing electricity from the external grid. This mechanism prioritizes incentivizing the local consumption and power exchange of new energy within multiple entities, improving the utilization efficiency of renewable energy and the overall economic benefits of the system. Through the combination of the gradient descent update and upper and lower limit cutoff mechanisms, the electricity price signal can both adaptively adjust to the system's supply and demand status and always remain within a reasonable range with clear physical and economic significance.

[0095] After the price update is completed, the virtual trading center will distribute the updated market transaction electricity price to each participant. Each participant will then perform a local two-stage stochastic optimization solution based on the new electricity price and enter the next iteration until complete convergence.

[0096] Example 2 This embodiment uses a microgrid group consisting of three geographically adjacent integrated energy systems with complementary load characteristics as the experimental object to verify the distributed stochastic optimization method for integrated energy systems described in this application.

[0097] Experimental conditions settings: The microgrid cluster comprises three main entities: Entity 1 is a photovoltaic-dominated integrated energy system, while Entities 2 and 3 are wind-dominated integrated energy systems. Each entity is equipped with renewable energy generation units, multi-energy conversion units, multi-energy storage units, and external interaction interfaces, providing the physical foundation for participating in multi-energy coordinated optimization scheduling and cross-entity power exchange. The system's day-ahead scheduling cycle is set to 24 hours, with a time step of 1 hour, resulting in 24 scheduling periods throughout the day.

[0098] To address source load uncertainty, typical scenarios were generated based on a 90-day historical operational dataset. First, probability distribution features were extracted from the historical data using kernel density estimation. Then, Monte Carlo sampling was used to generate an initial scenario set. Finally, the K-Medoids clustering algorithm was employed to jointly reduce the initial scenarios, ultimately generating five typical scenarios, each with a definite probability of occurrence, for use in the subsequent two-stage stochastic optimization model calculation.

[0099] Typical scenario generation results: Taking subject 1 as an example, the prediction residuals are first subjected to variational mode decomposition. For example... Figure 2As shown, the prediction residuals are decomposed into four sub-modes: low-frequency mode, mid-frequency mode, high-frequency mode, and residual components. The low-frequency mode exhibits gradual changes, reflecting systematic prediction shifts caused by factors such as seasonal environmental changes; the mid-frequency mode shows fluctuations with a 24-hour period, consistent with the day-night cycle and daily human activity patterns; the high-frequency mode exhibits drastic amplitude changes and occurs at random times, corresponding to random errors caused by occasional factors such as sudden gusts of wind and short-term cloud cover. Figure 3 As shown, after kernel density estimation for each sub-mode, the probability density curves of each mode exhibit their own unique distribution patterns, indicating that the errors in different frequency bands have different statistical characteristics, thus verifying the necessity of performing mode decomposition on the residuals to fit the probability densities separately.

[0100] Monte Carlo sampling was performed based on the above probability density model to generate 1000 initial predicted scenarios. Then, high-dimensional joint feature flattening and K-Medoids clustering reduction were used to obtain 5 typical scenarios. Figure 4 and Figure 5 As shown, the reduced typical scenarios cover the main distribution areas of the source-load uncertainty space. The differences between the scenarios are mainly reflected in the fluctuation range of new energy output and the combination mode of cooling, heating and power loads. This verifies that the proposed scenario generation method can significantly reduce the number of scenarios while retaining the main statistical characteristics of the original uncertainty space.

[0101] Optimization scheduling result analysis: The five generated typical scenarios and their occurrence probabilities are substituted into a two-stage stochastic optimization model for solution. The scheduling results are analyzed using subject 1 as an example.

[0102] like Figure 6 As shown, during the period from 00:00 to 06:00, the electrical load is at a low level, wind power output is sufficient, and the energy storage system is charging. Simultaneously, entity 1 sells excess electricity to other entities through point-to-point transactions, achieving cross-entity consumption of new energy. From 08:00 to 16:00, photovoltaic power output reaches its peak, but due to the simultaneous increase in electrical load, the output of new energy is still insufficient to fully meet the load demand. The energy storage system switches to a discharging state, and entity 1 purchases some electricity from the external grid to make up the gap. During the evening peak period from 18:00 to 22:00, photovoltaic output is zero, and although wind power output increases, it still cannot cover all load demand. At this time, the internal combustion engine starts generating electricity, the energy storage system further discharges, and electricity is purchased from other entities through point-to-point transactions to jointly support the electrical load demand. Figure 6The power dispatch balance shown indicates that the two-stage stochastic optimization model can flexibly allocate various adjustment methods such as internal unit output, energy storage charging and discharging, external grid interaction, and cross-entity point-to-point transactions at different times according to the source and load characteristics of each typical scenario, thereby achieving power supply and demand balance and economic optimization.

[0103] Figure 7 This is a cold energy scheduling balance diagram for main body 1 in this embodiment. For example... Figure 7 As shown, during the daytime hours of the summer cooling season, the demand for cooling load is high, and absorption chillers and electric chillers operate in tandem to meet the cooling load. The absorption chiller utilizes waste heat from the internal combustion engine to drive refrigeration, achieving cascaded utilization of thermal energy; the electric chiller increases its output during periods of low electricity prices or when renewable energy supply is abundant, thereby reducing overall cooling costs. The output ratio of the two cooling methods adaptively adjusts according to different typical scenarios, verifying the adaptability of the two-stage stochastic optimization model to the uncertainty of the cooling load.

[0104] Figure 8 This is the thermal energy scheduling balance diagram of main body 1 in this embodiment. For example... Figure 8 As shown, waste heat recovery from the internal combustion engine is one of the main sources of thermal energy for the system, while the gas boiler serves as a supplementary heat source, starting when waste heat is insufficient. The thermal storage tank is charged during periods of abundant waste heat and released heat during peak heat load periods, playing a role in peak shaving and valley filling. Under various typical scenarios, the thermal balance constraints are strictly satisfied, indicating that the optimization model can effectively coordinate multiple heat source devices and thermal storage devices to achieve a continuous and stable supply of thermal energy.

[0105] Figure 9 This is a diagram showing the electrical energy interaction power between the three main entities in this embodiment. For example... Figure 9 As shown, from 00:00 to 06:00, entities 2 and 3 have high wind power output and sell electricity to entity 1. From 10:00 to 16:00, entity 1 has sufficient photovoltaic output and sells electricity to entities 2 and 3. During the evening peak period, the renewable energy output of all entities is insufficient. At this time, entities mainly meet their needs through internal energy storage discharge and purchasing electricity from the external grid, resulting in a decrease in point-to-point transactions. Throughout the 24 time periods of the day, the point-to-point transaction price between entities is lower than the purchase price from the external grid, verifying the effectiveness of the price cap mechanism—this mechanism ensures the economic attractiveness of internal transactions and prioritizes the local consumption of renewable energy and the mutual assistance of electricity among entities. The sum of the power purchased and sold by the three entities in each time period is always zero, verifying the convergence of the distributed collaborative optimization method of alternating direction multiplier and the effectiveness of the global balance constraint. This indicates that by exchanging electricity information only, entities have successfully achieved distributed clearing of the multi-entity market while fully protecting the privacy of internal equipment parameters and load data.

[0106] The experimental results above show that the distributed stochastic optimization method for integrated energy systems described in this application can generate physically realistic, high-fidelity typical scenarios and effectively quantify source-load uncertainty; the constructed two-stage stochastic optimization model can achieve collaborative optimization of forward-looking trading plans and adaptive internal scheduling strategies under uncertain environments; and the distributed solution framework using the alternating direction multiplier method can achieve effective collaboration and stable mutual support among multiple entities while protecting the privacy of each entity.

[0107] The above description is merely a preferred embodiment of this application and an explanation of the technical principles employed. Those skilled in the art should understand that the scope of the invention involved in this application is not limited to technical solutions formed by specific combinations of the above-described technical features, but should also cover other technical solutions formed by arbitrary combinations of the above-described technical features or their equivalents without departing from the inventive concept. For example, technical solutions formed by substituting the above features with (but not limited to) technical features with similar functions disclosed in this application.

Claims

1. A distributed stochastic optimization method for an integrated energy system, characterized in that, include: Acquire historical data on multidimensional physical quantities in integrated energy systems; Based on the historical data, multiple typical scenarios and the probability of occurrence of each typical scenario are generated; A two-stage stochastic optimization model is constructed for each entity. The first-stage decision variable of the two-stage stochastic optimization model is the point-to-point cross-regional electricity transaction volume between the entities, and the second-stage decision variable is the power response and regulation status of the internal equipment of each entity under various typical scenarios. The alternating direction multiplier method is used to perform distributed collaborative optimization on multiple entities. The entities only interact with each other through point-to-point electricity trading information. The solution is iteratively solved until the global clearing convergence condition is met, and the optimized scheduling scheme of each entity is output.

2. The distributed stochastic optimization method for integrated energy systems according to claim 1, characterized in that, The process of generating multiple typical scenarios and the probability of occurrence of each typical scenario based on the historical data includes: A dual-branch prediction network is constructed for the multidimensional physical quantity. The dual-branch prediction network includes a continuous value regression main network and a start-stop classification auxiliary network. The main network outputs continuous prediction values ​​of the multidimensional physical quantity, and the auxiliary network outputs the start-stop state probabilities of the cooling load and heating load and constructs a start-stop mask matrix. The continuous prediction values ​​are fused with the start-stop mask matrix to obtain the corrected baseline prediction values. Calculate the prediction residual between the true value of the multidimensional physical quantity and the baseline prediction value, and perform mode decomposition on the prediction residual to obtain multiple residual sub-modes in different frequency bands; Each residual sub-mode is fitted with a probability density, and random sampling is performed based on the fitted probability density model. The sampling results are superimposed on the baseline prediction value and the start / stop mask matrix is ​​applied to generate multiple initial prediction scenarios. High-dimensional joint feature flattening is performed on each of the initial prediction scenarios, and the flattened feature vectors are clustered and reduced to obtain a set of typical scenarios and the occurrence probability of each typical scenario.

3. The distributed stochastic optimization method for integrated energy systems according to claim 2, characterized in that, The continuous value regression main network is constructed based on a long short-term memory network, and the start-stop classification auxiliary network is constructed based on a long short-term memory network and a sigmoid activation function; the fusion process involves performing Hadamard product operations on the continuous predicted values ​​of cooling load and heating load with the corresponding start-stop state probability masks, respectively.

4. The distributed stochastic optimization method for integrated energy systems according to claim 2, characterized in that, The modal decomposition of the predicted residual includes: Variational mode decomposition is used to decompose the prediction residual into low-frequency mode, mid-frequency mode, high-frequency mode and residual component; the low-frequency mode represents the long-term trend error in the prediction deviation, the mid-frequency mode represents the intraday periodic fluctuation error in the prediction deviation, and the high-frequency mode represents the random mutation error in the prediction deviation.

5. The distributed stochastic optimization method for integrated energy systems according to claim 2, characterized in that, The probability density of each residual sub-mode is fitted using the kernel density estimation method, and the initial prediction scene is generated using the Monte Carlo sampling method. The high-dimensional joint feature flattening is to concatenate the data of all time steps and all physical dimensions contained in each initial prediction scene to form a single high-dimensional joint feature row vector. The clustering reduction adopts the K-Medoids clustering algorithm.

6. The distributed stochastic optimization method for integrated energy systems according to claim 1, characterized in that, The construction of a two-stage stochastic optimization model for each subject includes: The optimization objective is to minimize the expected total cost across all scenarios throughout the day. The expected total cost includes the expected cost of scenario operation and the cost of point-to-point electricity trading. The expected cost of scenario operation includes the cost of purchasing electricity from the external grid, the cost of curtailment penalties, and the cost of natural gas fuel.

7. The distributed stochastic optimization method for integrated energy systems according to claim 1, characterized in that, The method of using alternating direction multipliers to perform distributed collaborative optimization of multiple agents includes: In each iteration of the alternating direction multiplier method, each entity sends its expected transaction power vector to the virtual trading center. The virtual trading center then aggregates the data and sends the sum of transaction power information of all entities except the entity itself to each entity. Each entity updates its augmented Lagrange objective function based on the received information and solves it again.

8. The distributed stochastic optimization method for integrated energy systems according to claim 7, characterized in that, The augmented Lagrange objective function includes a proximal damping term, which is used to limit the magnitude of change in the current round of trading decisions relative to the previous round of trading decisions.

9. The distributed stochastic optimization method for integrated energy systems according to claim 7, characterized in that, The method for determining the global clearing convergence condition is as follows: After collecting the trading decisions returned by all participants, the virtual trading center calculates the sum of the net trading gaps in the entire market for each time period as the trading residual, and determines whether the trading residual is less than a preset convergence threshold. If it is satisfied, it is determined that the multi-partner market has achieved distributed clearing, the algorithm terminates and outputs an optimized scheduling scheme; if it is not satisfied, the Lagrange multipliers are updated and the next round of iteration begins.

10. The distributed stochastic optimization method for integrated energy systems according to claim 9, characterized in that, The Lagrange multiplier is used as the market transaction price, and it is updated using the gradient descent method. An upper limit and a lower limit are set for the updated market transaction price. The lower limit is 0, and the upper limit is the difference between the external grid purchase price and the preset price margin, so that the internal transaction price between entities is lower than the price of directly purchasing electricity from the external grid.

Images