Real-time decision-making method and device for industrial integrated energy system, and storage medium

Abstract

The invention relates to a real-time decision-making method and device for an industrial integrated energy system and a storage medium, and belongs to the technical field of industrial automation and energy system control. According to the method, the hybrid digital twinborn model formed by coupling the physical mechanism module and the data driving module is constructed, and the multi-time-scale renewable energy output and multi-element load probability distribution scene set is generated in combination with the conditional generative adversarial network. The fitness of the model to the physical law of the system is guaranteed, and the uncertainty of energy supply and demand can be accurately captured; meanwhile, a preliminary control action is generated through a pre-trained MADDPG network, rolling optimization is carried out by relying on a hybrid digital twin model, and real-time safety verification is carried out by introducing a control barrier function, so that economical efficiency optimization and safety guarantee of the control action are realized, the problems that a traditional control response is lagged, and safety constraints are easy to break through are solved, and the control effect is improved. And the autonomous operation capability and stability of the system are improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of industrial automation and energy system control, and particularly relates to an industrial comprehensive energy system real-time decision method and device and a storage medium. BACKGROUND

[0002] An industrial comprehensive energy system integrates various energy forms such as electricity, heat and gas, and its efficient, safe and economic operation is crucial to industrial production. Traditional industrial energy system control methods usually rely on scheduling based on fixed rules or simplified physical models. In the face of intermittent renewable energy output and dynamic changes in multi-element loads, it is difficult to make accurate predictions and effective scheduling, often leading to energy supply and demand imbalance, high operating costs or energy waste and other problems.

[0003] With the development of artificial intelligence technology, some more advanced control methods have been proposed. For example, some technical solutions use a digital twin model to simulate system behavior, which can couple physical mechanism models and data-driven models to improve prediction accuracy; or use multi-agent reinforcement learning algorithms to generate control strategies; some solutions use models such as autoencoders to calculate reconstruction errors to monitor abnormal states of system operation. However, these existing technical solutions still have obvious defects in actual application: first, most methods do not handle the uncertainty factors in the system (such as weather changes and load fluctuations) well, usually using deterministic predictions, resulting in control strategies that are not robust in real fluctuating scenarios; second, some methods lack effective real-time safety verification mechanisms while pursuing operational economy, which may generate dangerous instructions that break through the physical limits of devices or the safety boundaries of the system, posing a safety hazard. SUMMARY

[0004] To this end, the application provides an industrial comprehensive energy system real-time decision method, device and storage medium, aiming to solve the problem of how to balance economy, safety and adaptability in complex scenarios with high uncertainty and volatility of renewable energy output and multi-element loads, and to realize closed-loop autonomous optimization control of industrial comprehensive energy systems.

[0005] According to a first aspect of an embodiment of the application, an industrial comprehensive energy system real-time decision method is provided, the method comprising:

[0006] Obtaining real-time operation data of an industrial comprehensive energy system, and preprocessing the real-time operation data to obtain a standardized multi-source time series data set;

[0007] Based on the standardized multi-source time series data set, a pre-trained probability distribution scenario prediction model is used to generate a probability distribution scenario set of renewable energy output and multi-element loads at different prediction time scales;

[0008] based on the standardized multi-source time series dataset and the probability distribution scenario set, a pre-trained MADDPG prediction model is used to generate preliminary control actions of each controllable device in the industrial integrated energy system;

[0009] based on the standardized multi-source time series dataset, a pre-trained hybrid digital twin model is used to predict the future system state vector of the industrial integrated energy system, wherein the hybrid digital twin model couples the output of a physical mechanism module constructed by differential algebraic equations and a data-driven module constructed by a time series convolution network;

[0010] The preliminary control actions of the controllable devices are used as initial control sequences, and a rolling optimization is performed to minimize the total operating cost of the system. In the optimization process, a control barrier function is introduced, and the initial control sequence is checked in real time for safety by the system state vector output by the pre-trained hybrid digital twin model at each optimization step, to generate a final control instruction sequence.

[0011] Preferably,

[0012] The training of the probability distribution scenario prediction model includes:

[0013] A conditional generative adversarial network is constructed, which consists of a generator and a discriminator;

[0014] The historical window length and the future prediction step of three different prediction time scales, i.e., ultra-short term, short term and day-ahead, are set respectively;

[0015] The historical standardized multi-source time series dataset is obtained, and the conditional information of the historical window length is extracted therefrom. The conditional information and a random noise vector are used as inputs of the generator, and a probability distribution scenario of the future prediction step is output;

[0016] Real data is extracted from the historical standardized multi-source time series dataset, and the probability that the real data and the prediction data of the generator are real data is judged by the discriminator;

[0017] The difference between the generated prediction data and the real data is calculated by an adversarial loss function, and an adversarial loss value is obtained;

[0018] The conditional generative adversarial network is trained using the historical standardized multi-source time series dataset. The generator and the discriminator are alternately optimized during the training process, and a physical consistency loss value is introduced. The adversarial loss value and the physical consistency loss value are jointly optimized to minimize the total loss. A trained probability distribution scenario prediction model is obtained;

[0019] The generating renewable energy output and multi-element load probability distribution scenario set of different prediction time scales by using the pre-trained probability distribution scenario prediction model comprises:

[0020] For future ultra-short-term, short-term and day-ahead time scales, a plurality of probability distribution scenarios are generated by the generator of the pre-trained probability distribution scenario prediction model based on different future prediction steps set;

[0021] The K-means clustering method is used to reduce the probability distribution scenarios generated at each time scale, the Euclidean distance between scenarios is used as the similarity measure, the minimization of the within-class sum of squares is used as the optimization objective, the center scenarios of each class are selected as the representative scenario subset, and the proportion of the number of initial scenarios contained in the class to the total number of scenarios is used as the occurrence probability, forming the reduced probability distribution scenario set.

[0022] Preferably,

[0023] The training of the MADDPG prediction model comprises:

[0024] Each controllable device in the industrial integrated energy system is regarded as an agent, an independent Actor network is established for each agent, and a centralized Critic network is established to guide the policy update of each Actor network, thereby constructing a MADDPG network;

[0025] Local observation variables of each agent are extracted from a historical standardized multi-source time series data set;

[0026] Based on the local observation variables of each agent, the Actor network receives the local observation variables of a single agent, outputs the continuous control action of the agent, and constitutes the joint action space of all agents;

[0027] The energy procurement cost, the device start-stop loss cost and the carbon trading cost are obtained respectively, and the total system operation cost is calculated; the reward function is set as the negative value of the total system operation cost;

[0028] An offline reinforcement learning algorithm is used to pre-train the MADDPG network, a four-tuple containing state, action, reward and next state is constructed using historical operation data, the four-tuple is used to train the MADDPG network; the reward function is used as the benchmark of the Q value of the Critic network, the policy gradient of the Actor network is updated by minimizing the time difference error of the Critic network to optimize the network parameters of the MADDPG network, the network parameters include the connection weights and neuron biases in the Actor network and the Critic network; a pre-trained MADDPG prediction model is obtained;

[0029] The generating the preliminary control actions of the controllable devices in the industrial integrated energy system by using the pre-trained MADDPG prediction model comprises:

[0030] The standardized multi-source time series data set and the probability distribution scenario set are input into the pre-trained MADDPG prediction model, for each probability distribution scenario, a joint action space of each agent is generated by an Actor network, and a target Q value of each action is evaluated by a Critic network, and the action with the highest target Q value is selected as the preliminary action under the scenario;

[0031] The preliminary actions under all probability distribution scenarios are weighted in expectation to obtain the preliminary control actions of each controllable device.

[0032] Preferably,

[0033] The output of the physical mechanism module comprises:

[0034] Based on the standardized multi-source time series data set, the current state variable of the physical mechanism dominant unit is obtained, and the control variable related to the operation of the physical mechanism dominant unit is extracted;

[0035] Based on the current state variable and the control variable, a set of physical parameters of the physical mechanism dominant unit is introduced, and by solving a system of differential algebraic equations describing the dynamic process and static relationship of the physical mechanism dominant unit, the state variable of the physical mechanism dominant unit at the next moment is obtained;

[0036] Based on the state variable at the next moment and the set of physical parameters, the output variable at the current moment is output by calculating the algebraic constraint equation in the system of differential algebraic equations.

[0037] Preferably,

[0038] The output of the data-driven module comprises:

[0039] Based on the standardized multi-source time series data set, the historical multi-dimensional time series feature data of the data-driven dominant unit is obtained, and the output variable at the current moment provided by the physical mechanism module is fused as the input data of the data-driven module;

[0040] The data-driven module extracts features of the input data by an expansion convolution layer of a time series convolution network to obtain a high-dimensional feature sequence;

[0041] The data-driven module performs nonlinear transformation and feature fusion on the high-dimensional feature sequence by a residual connection layer of the time series convolution network to obtain a deep time series feature representation;

[0042] The data-driven module performs linear mapping on the deep time series feature representation by a fully connected layer to output a predicted value at the next moment;

[0043] The method of using a pre-trained hybrid digital twin model to predict the future system state vector of the industrial integrated energy system includes:

[0044] The predicted value for the next moment output by the data-driven module is integrated with the state variable for the next moment obtained by the physical mechanism module to output the future system state vector of the industrial integrated energy system.

[0045] Preferably,

[0046] The training of the hybrid digital twin model includes:

[0047] The objective function of the hybrid digital twin model is set as follows: the objective function is to minimize the mean square error between the output of the hybrid digital twin model and the actual measured value.

[0048] The hybrid digital twin model is iteratively trained using historical operating data, including: adjusting the physical parameters of the physical mechanism module using the maximum likelihood estimation method; adjusting the network parameters of the data-driven module using the backpropagation algorithm based on stochastic gradient descent; and obtaining the trained hybrid digital twin model by iteratively training it to satisfy the objective function.

[0049] Preferably,

[0050] The rolling optimization, which uses the initial control actions of the controllable device as the initial control sequence and aims to minimize the total system operating cost, includes:

[0051] Set the rolling optimization time domain and the control time domain, wherein the rolling optimization time domain is greater than the control time domain;

[0052] In each control cycle, an objective function is constructed with the goal of minimizing the total system operating cost within the rolling optimization time domain.

[0053] Using the actual operating state of the system reflected by the standardized multi-source time series dataset as the system state vector at the current moment, and the initial control action in the control time domain as the initial control sequence, multi-step forward simulation is performed through the pre-trained hybrid digital twin model to obtain the state prediction sequence in the rolling optimization time domain.

[0054] During the optimization process, system operation constraints are introduced; these constraints include gas turbine output power constraints, energy storage device charging and discharging power constraints, energy storage device charge constraints, adjustable load power constraints, power bus voltage constraints, and thermal pipeline pressure constraints.

[0055] The sequential quadratic programming algorithm is used to solve the objective function optimization problem with system operation constraints, and the optimized control sequence is obtained.

[0056] Preferably,

[0057] The process of introducing a control barrier function during optimization, and using the system state vector of each optimization step output by the pre-trained hybrid digital twin model to perform real-time safety verification on the initial control sequence, and generating the final control command sequence includes:

[0058] Based on the optimized control sequence, the first control action to be executed in the current control cycle is extracted as a candidate control instruction to be verified.

[0059] Define the key safety variables in the industrial integrated energy system, and define control barrier functions respectively in conjunction with the system operation constraints;

[0060] Substitute the current system state vector into each control barrier function to calculate the function value of each critical safety variable at the current moment;

[0061] The system state vector at the next moment is predicted by performing single-step forward simulation using the pre-trained hybrid digital twin model, and the function value of each key safety variable at the next moment is calculated by using the control barrier function.

[0062] If the function value of each key safety variable in the next time step is not lower than a safety threshold calculated from the function value at the current time step and the corresponding relaxation coefficient, then the current candidate control instruction is determined to be safe and can be directly used as the final instruction; otherwise, it is determined to be unsafe, and the current candidate control instruction is corrected to obtain the corrected safe control instruction.

[0063] The corrected safety control command will be used as the final control command for the current moment.

[0064] Each control action in the optimized control sequence is subjected to safety judgment and correction in turn, and the resulting final control commands are sorted in turn to obtain the final control command sequence.

[0065] Preferably,

[0066] The step of correcting the current candidate control command to obtain the corrected safety control command includes:

[0067] The candidate control commands are corrected by minimizing perturbation using a sequential quadratic programming algorithm. This is done to construct and solve a quadratic programming problem with the goal of minimizing the control command correction amount and the constraint of satisfying the safety conditions of the control barrier function, so as to obtain the corrected safety control command.

[0068] Preferably, it further includes:

[0069] The final control command sequence is sent to the physical system for execution, and system operation data is collected synchronously. The reconstruction error is calculated by a pre-trained autoencoder. When the reconstruction error continues to exceed a preset error threshold, the parameters of the data-driven module and the MADDPG prediction model in the hybrid digital twin model are updated.

[0070] Preferably,

[0071] The calculation of reconstruction error using a pre-trained autoencoder includes:

[0072] An autoencoder network consisting of an encoder and a decoder is constructed. The encoder compresses the standardized multi-source time-series data of historical moments from D dimensions to a d-dimensional latent space, and the decoder obtains the reconstructed data of historical moments.

[0073] Based on historical standardized multi-source time series data, the autoencoder network is pre-trained in an unsupervised manner by minimizing the mean square error between the input data and the reconstructed data, thus obtaining a trained autoencoder.

[0074] Real-time standardized multi-source time-series data is input into the trained autoencoder to obtain the reconstruction error;

[0075] The reconstruction error distribution of statistically standardized multi-source time-series data on the trained autoencoder is used to obtain the error threshold.

[0076] If the reconstruction error is greater than the error threshold for a consecutive preset time period, the data-driven module and the MADDPG prediction model in the hybrid digital twin model will be updated.

[0077] When the model parameter update is triggered, the incremental learning process is started, including: obtaining the system operation data of the most recent preset time period, wherein the system operation data includes the standardized multi-source time series dataset of the most recent preset time period and the corresponding actual measurement values;

[0078] The input features and their corresponding actual outputs related to the data-driven module in the hybrid digital twin model are extracted from the standardized multi-source time-series dataset of the most recent preset time period. Incremental training is performed using mini-batch gradient descent to update the parameters of the temporal convolutional network, thereby minimizing the mean square error between the output of the hybrid digital twin model and the actual measured value.

[0079] An experience replay buffer is constructed to store quadruples, and the parameters of the MADDPG network are updated using offline reinforcement learning through the buffer.

[0080] According to a second aspect of the present invention, a real-time decision-making device for an integrated industrial energy system is provided, the device comprising:

[0081] Data acquisition module: used to acquire real-time operating data of the industrial integrated energy system, and preprocess the real-time operating data to obtain a standardized multi-source time-series dataset;

[0082] Probability distribution scenario prediction module: Based on the standardized multi-source time series dataset, it uses a pre-trained probability distribution scenario prediction model to generate probability distribution scenario sets of renewable energy output and multi-load at different prediction time scales;

[0083] Preliminary control action generation module: Based on the standardized multi-source time-series dataset and probability distribution scenario set, it uses a pre-trained MADDPG prediction model to generate preliminary control actions for each controllable device in the industrial integrated energy system.

[0084] Future system state vector prediction module: Based on the standardized multi-source time-series dataset, it uses a pre-trained hybrid digital twin model to predict the future system state vector of the industrial integrated energy system. The hybrid digital twin model is coupled with the output of a physical mechanism module constructed by differential algebraic equations and a data-driven module constructed by temporal convolutional networks.

[0085] Control command generation module: It is used to perform rolling optimization with the initial control action of the controllable device as the initial control sequence and the goal of minimizing the total system operating cost. During the optimization process, a control barrier function is introduced. The initial control sequence is verified in real time by the system state vector of each optimization step output by the pre-trained hybrid digital twin model, and the final control command sequence is generated.

[0086] According to a third aspect of the present invention, a storage medium is provided, the storage medium storing a computer program, which, when executed by a host controller, implements the steps of the above-described method.

[0087] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects:

[0088] This application constructs a hybrid digital twin model coupling a physical mechanism module and a data-driven module, and combines it with a conditional generative adversarial network to generate a set of renewable energy output and multi-load probability distribution scenarios across multiple time scales. Compared with traditional schemes that rely on fixed physical models or single prediction methods, this approach ensures the model's close adherence to the physical laws of the system and accurately captures the uncertainty of energy supply and demand. It solves the problems of traditional methods having limited prediction scenarios and difficulty in handling complex fluctuations, providing comprehensive and scientific scenario support for subsequent decision-making. Simultaneously, it generates preliminary control actions through a pre-trained MADDPG network, and then performs rolling optimization based on a model predictive control framework with the hybrid digital twin model as the internal prediction model. It also introduces a control barrier function for real-time safety verification. Compared with the traditional control mode that relies mainly on manual intervention and lacks dynamic optimization and safety closed-loop, this approach achieves economic optimization and safety assurance of control actions, solves the problems of delayed response and easy breach of safety constraints in traditional control, and improves the system's autonomous operation capability and stability.

[0089] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and do not limit this application. Attached Figure Description

[0090] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0091] Figure 1 This is a flowchart illustrating a real-time decision-making method for an integrated industrial energy system according to an exemplary embodiment;

[0092] Figure 2 This is a system schematic diagram of a real-time decision-making device for an integrated industrial energy system according to another exemplary embodiment;

[0093] In the attached diagram: 1-Data acquisition module, 2-Probability distribution scenario prediction module, 3-Preliminary control action generation module, 4-Future system state vector prediction module, 5-Control command generation module. Detailed Implementation

[0094] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be described in detail below. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Based on the embodiments in this application, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0095] Example 1

[0096] Figure 1 This is a flowchart illustrating a real-time decision-making method for an integrated industrial energy system according to an exemplary embodiment, such as... Figure 1 As shown, the method includes:

[0097] S1, acquire real-time operating data of the industrial integrated energy system, and preprocess the real-time operating data to obtain a standardized multi-source time-series dataset;

[0098] S2, Based on the standardized multi-source time-series dataset, a pre-trained probability distribution scenario prediction model is used to generate probability distribution scenario sets of renewable energy output and multi-load at different prediction time scales;

[0099] S3. Based on the standardized multi-source time-series dataset and probability distribution scenario set, the pre-trained MADDPG prediction model is used to generate the preliminary control actions of each controllable device in the industrial integrated energy system.

[0100] S4. Based on the standardized multi-source time-series dataset, a pre-trained hybrid digital twin model is used to predict the future system state vector of the industrial integrated energy system. The hybrid digital twin model is coupled with the output of a physical mechanism module constructed by differential algebraic equations and a data-driven module constructed by a temporal convolutional network.

[0101] S5, taking the initial control action of the controllable device as the initial control sequence, and performing rolling optimization with the goal of minimizing the total operating cost of the system, a control barrier function is introduced during the optimization process, and the initial control sequence is verified in real time by the system state vector of each optimization step output by the pre-trained hybrid digital twin model, and the final control command sequence is generated.

[0102] It is understandable that in this embodiment,

[0103] Step S1: Obtain real-time operating data of the industrial integrated energy system, and preprocess the real-time operating data to obtain a standardized multi-source time-series dataset, specifically including:

[0104] Sensor arrays are deployed at key topology nodes of power supply networks, heat transmission networks, refrigeration cycle networks, and gas supply networks. The sensor arrays include voltage transformers, current transformers, temperature sensors, pressure transmitters, flow meters, and equipment status sensors to collect real-time operating data of the integrated industrial energy system.

[0105] The key topology nodes include network structural locations such as busbars, pipeline junctions, energy storage access points, and controllable load interfaces;

[0106] The real-time operation data of the industrial integrated energy system includes data in three dimensions: electrical operation parameters, flow parameters, and equipment status.

[0107] Among them, electrical operating parameters include at least voltage, current, active power, reactive power, grid frequency, and power quality data;

[0108] Flow parameters include at least temperature, pressure, mass flow rate, heat flow rate, and the enthalpy and entropy of the medium;

[0109] The equipment status includes at least the rotational speed of rotating equipment, the charge status of energy storage devices, and the open / closed status of circuit breakers, disconnect switches, and regulating valves;

[0110] After obtaining the real-time operation data of the aforementioned integrated industrial energy system, timestamp alignment and normalization preprocessing are required to obtain a standardized multi-source time-series dataset, specifically including:

[0111] Based on real-time operational data, the original sampling frequency of each sensor data is identified. Using the highest sampling frequency as a benchmark, a fixed time window interpolation method is adopted to uniformly interpolate the data of each sensor to the same time series.

[0112] Among them, the fixed time window interpolation method refers to calculating the estimated values ​​at all intermediate reference time points by fitting a straight line between the measured values ​​at two known time points before and after the time point.

[0113] Based on the time-aligned data, the sliding window Z-score method is used to detect outliers, and the outliers are replaced by linear interpolation of the normal data before and after them.

[0114] The sliding window Z-score method calculates the mean and standard deviation of all data within a fixed-length data window. For each data point within the window, it calculates the multiple of the difference between the data point and the mean relative to the standard deviation. If the absolute value of this multiple exceeds a set threshold, the data point is considered an outlier. For outlier data points, they are replaced and filled using the adjacent normal data values ​​through linear calculation. If outliers occur consecutively, the average value of historical data from the same period is used for filling.

[0115] After outlier processing, the data is subjected to min-max normalization according to its type. The normalized data is then arranged in chronological order to obtain the standardized multi-source time series dataset {X(t)}, where t is the timestamp, expressed by the following formula:

[0116] X(t)=[U norm (t),I norm (t),P norm (t),Tnorm (t),p norm (t),

[0117]

[0118] In the formula, U norm (t) is the normalized voltage, I norm (t) represents the normalized current, P norm (t) represents the normalized active power, T norm (t) represents the normalized temperature, p norm (t) represents the normalized pressure. For normalized mass flow rate, SOC norm (t) represents the normalized energy storage state of charge, valve state {X(t)} represents the normalized valve state, etc.; {X(t)} is the final standardized multi-source time series dataset with dimension D.

[0119] Step S2: Based on the standardized multi-source time-series dataset, a pre-trained probability distribution scenario prediction model is used to generate probability distribution scenario sets for renewable energy output and diverse loads at different prediction time scales, specifically including:

[0120] A conditional generative adversarial network is constructed and trained to obtain a pre-trained probability distribution scene prediction model; wherein the conditional generative adversarial network consists of a generator G and a discriminator D, both of which are deep neural network structures;

[0121] The generator G adopts a feedforward network architecture based on fully connected layers and one-dimensional transposed convolutional layers. Its input is a concatenated random noise vector and a conditional information vector. After feature mapping through the fully connected layer, it is temporally upsampled through two levels of one-dimensional transposed convolutional layers. Finally, it generates time series data that conforms to the target dimension through a one-dimensional convolutional output layer, and uses the Tanh activation function to constrain the output value within the normalization range.

[0122] Discriminator D adopts a discriminative network architecture based on one-dimensional convolutional layers and fully connected layers. It concatenates the input data with the conditional features projected by the fully connected layer in the channel dimension, performs feature extraction and downsampling through two levels of one-dimensional convolutional layers, integrates the features through the fully connected layer, and finally outputs the probability value representing the authenticity of the data through the Sigmoid activation function.

[0123] The generator G is responsible for learning the conditional probability distribution and generating scene data, while the discriminator D is responsible for evaluating the authenticity of the data. The two are jointly optimized through adversarial training, and the specific training process includes:

[0124] Based on three different forecast time scales—ultra-short-term, short-term, and day-ahead—different historical window lengths T0 and future forecast step sizes T are set. f ;

[0125] Specifically, when used for ultra-short-term forecasting (15 minutes to 2 hours into the future), if the data sampling frequency is 15 minutes per step, then the historical window length T0 is 96 (covering the previous 24 hours of historical data), and the future forecasting step size T... f The value ranges from 1 to 8 depending on the specific prediction duration (corresponding to 15 minutes to 2 hours);

[0126] When used for short-term forecasting (6 to 24 hours from now), if the data sampling frequency is 15 minutes per step, then the historical window length T0 is 192 (covering the previous 48 hours of historical data), and the future forecast step size T... f The value is between 24 and 96 depending on the specific forecast duration (corresponding to 6 to 24 hours);

[0127] When used for day-ahead forecasting (future 24 hours), if the data sampling frequency is 1 hour / step, then the historical window length T0 is 168 (covering historical data from the previous week), and the future forecast step size T... f Take 24 (corresponding to 24 hours);

[0128] Renewable energy output data and multi-load data for historical T0 time steps are extracted from a standardized multi-source time series dataset as conditional information c. Conditional information c and random noise vector z are used as input to generator G to output the probability distribution scenario of renewable energy output and multi-load in future periods. Conditional information c includes time series data of solar irradiance, ambient temperature, wind speed, electrical load, heat load and cooling load.

[0129] Where c is the conditional information, corresponding to the standardized multi-source time series data X(t-T0:t) of the historical T0 time steps, which includes key factors affecting energy output and load such as solar irradiance, ambient temperature, and wind speed, and is used to constrain the generator G to generate scenarios that meet specific conditions; z is a random noise vector, which satisfies z~N(0,1) (a normal distribution with a mean of 0 and a variance of 1) in the scenario. Its role is to provide randomness to the generator G and ensure that the generated scenarios are diverse;

[0130] Historical renewable energy output and multi-source load data are extracted from a standardized historical multi-source time-series dataset as real data x. Combined with conditional information c, the discriminator D determines the probability that the real data x and the output data of the generator G are real data. The real data x includes historical photovoltaic output, wind power output, electrical load, heat load and cooling load data.

[0131] Where x represents real data, specifically referring to historical standardized renewable energy output (photovoltaic power generation, wind power generation) and multi-source (electricity, heat, cooling) load data;

[0132] The adversarial loss value L is obtained by calculating the degree of difference between the generated data (i.e., the probability distribution scenario) and the real data using the adversarial loss function. cgan ;

[0133] The expression for the adversarial loss function is:

[0134] L cgan =E[logD(x|c)]+E[log(1-D(G(z|c)|c))]

[0135] In the formula, L cgan To counteract the loss value, the smaller the loss value, the closer the generated data is to the true data distribution; E(*) is the mathematical expectation operator, log( * ) is the logarithmic function; G(*) is the generator's output function, where G(z|c) is the generated data obtained by the generator after inputting random noise z under given conditional information c, i.e., the probability distribution scenario of future renewable energy output and diversified loads; D(*) is the discriminator's output function, representing the discriminator's judgment result on the input data, with an output value between 0 and 1. The closer to 1, the more likely the discriminator considers the input data to be real data; D(x|c) is the probability that the discriminator judges real data x as "real data" under given conditional information c;

[0136] The conditional generative adversarial network was trained using a historically normalized multi-source time-series dataset. During training, the generator G and discriminator D were alternately optimized, and a physical consistency loss value L was introduced. phy (To improve the physical plausibility of the generated data), jointly optimize the adversarial loss value L. cgan Loss value L of physical consistency phy To minimize the total loss L total ;

[0137] L phy =‖Y sim -Y scenario || 2

[0138] In the formula, L phy Y represents the physical consistency loss value, reflecting the deviation between the generated data (i.e., the probability distribution scenario) and the simulation results of the hybrid digital twin model. The smaller the loss value, the more the results conform to the physical operating laws of the system; sim For the simulation output of the hybrid digital twin model under the same condition information c, Y scenario The scene output generated by the generator for a single noise z; ||*| is the Euclidean norm;

[0139] L total =L cgan +λ·L phy

[0140] Among them, L total The total loss value for training the conditional generative adversarial network is the core objective of optimizing the generator G and discriminator D. During training, the total loss value needs to be minimized by adjusting network parameters (such as generator weights and discriminator biases). λ is a tradeoff coefficient that is set in advance according to actual needs.

[0141] After training, for three time scales—ultra-short term (minutes), short term (hours), and near term (24 hours)—different future prediction step sizes T are used. f A generator G generates several probability distribution scenarios (e.g., 1000 initial probability distribution scenarios for each time scale); each probability distribution scenario contains the predicted photovoltaic output P. pv (t+1:t+T f ), wind power output P wind (t+1:t+T f ), electrical load L elec (t+1:t+T f ), heat load L heat (t+1:t+T f ), cooling load L cool (t+1:t+T f );

[0142] Among them, the pre-level network of generator G (including fully connected layers and transposed convolutional layers) serves as a shared feature encoder. Its role is to receive the concatenated input vector [z,c] and extract a unified, high-level shared feature tensor that contains the potential laws of the future system's operating state.

[0143] Following the shared feature encoder, the network is connected in parallel to five independent dedicated output heads, each of which is a lightweight one-dimensional convolutional layer specifically designed to decode future sequences of specific variables from the shared feature tensor. These five output heads are:

[0144] The photovoltaic output head is responsible for outputting future T f The photovoltaic power output prediction sequence P at each time step pv (t+1:t+T f );

[0145] The wind power output head is responsible for outputting future T f Wind power output prediction sequence P at each time step wind (t+1:t+T f);

[0146] The electrical load output head is responsible for outputting the future T f The electrical load prediction sequence L at each time step elec (t+1:t+T f );

[0147] The heat load output head is responsible for outputting the future T f The heat load prediction sequence L for each time step heat (t+1:t+T f );

[0148] The cold load output head is responsible for outputting the future T f The cooling load prediction sequence L at each time step cool (t+1:t+T f );

[0149] This multi-output head structure offers significant advantages over using a single network to sequentially generate variables: Since all variables are synchronously decoded from the same set of shared features, this structure naturally captures and maintains the inherent spatiotemporal coupling and synchronization characteristics between various energy outputs (electricity, heat, cooling, etc.) and loads; shared feature encoding reduces the total number of model parameters and makes learning the joint distribution of multiple related variables more efficient; this structure allows for adjustments to the time step T of the final output layer for different time scales, such as ultra-short-term, short-term, and day-ahead. f This allows for the flexible generation of probabilistic scenarios of appropriate length without altering the core network structure.

[0150] To reduce computational complexity, K-means clustering is used to reduce the probability distribution scenes generated at each time scale. Euclidean distance between scenes is used as a similarity measure, and minimizing the sum of squares within each cluster is the optimization objective. The central scene of each category is selected as a representative subset of scenes, and the probability of occurrence is determined by the proportion of the initial scenes contained within that cluster to the total number of scenes, thus forming a simplified probability distribution scene set {Y}. scenario (k)};where k is the scene index, k=1,...,K, and K is the total number of scenes;

[0151] The K-means clustering method calculates the distance between samples, automatically groups a large number of scenes into K representative categories, and extracts the central scene of each category and its probability of occurrence to achieve dimensionality reduction and simplification of the scene set. The total number of scenes K is determined according to the variance percentage criterion, that is, the minimum number of clusters that makes the ratio of the sum of squares within a class to the total sum of squares reach a preset percentage (e.g., 85%) is selected.

[0152] Step S3: Based on the standardized multi-source time-series dataset and probability distribution scenario set, the pre-trained MADDPG prediction model is used to generate preliminary control actions for each controllable device in the industrial integrated energy system, specifically including:

[0153] Each controllable device in the industrial integrated energy system is set as an intelligent agent, an independent Actor network is established for each intelligent agent, and a centralized Critic network is established to guide the policy updates of each Actor network, thereby constructing the MADDPG network; among them, controllable devices (i.e. key devices in the system that can be directly controlled to execute optimization strategies) include gas turbines, energy storage devices, and adjustable loads, etc.

[0154] The Actor network receives local observation information from a single agent and outputs the agent's continuous control actions (such as valve opening and power setpoint); the Critic network receives joint observation information and joint actions from all agents and evaluates the long-term expected reward (Q value) of the current joint action in the global state to guide the policy update of each Actor network.

[0155] Local observation variables of each agent are extracted from the standardized multi-source time series dataset {X(t)}; among which, the local observation variables of each agent include the speed and exhaust temperature of the gas turbine, the state of charge and power value of the energy storage device, and the power level of the adjustable load;

[0156] Based on the local observation variables of each agent, the local observation variables of a single agent are received through the Actor network, and the continuous control actions of the agent are output to form the joint action space a(t) of all agents; wherein, the joint action space includes the fuel valve opening adjustment of the gas turbine, the power setpoint of the energy storage device, and the power reduction of the adjustable load.

[0157] a(t) = [a1(t), a2(t), a3(t)]

[0158] In the formula, a(t) is the joint action space of all intelligent agents, a1(t) represents the adjustment amount of the fuel valve opening of the gas turbine, a2(t) represents the power setting value of the energy storage device, and a3(t) represents the power reduction amount of the adjustable load.

[0159] By statistically analyzing the total amount of electricity and gas purchased within and outside the control cycle at time t, and multiplying it by the corresponding real-time electricity and gas prices, the energy procurement cost C is obtained. purchase (t);

[0160] C purchase (t)=[P grid (t)·π grid (t)+P gas (t)·π gas(t)]·Δt

[0161] In the formula, C purchase (t) represents the energy procurement cost (unit: yuan); P grid (t) represents the system's purchased power (unit: kW) at time t, i.e., the power purchased from the external power grid. A positive value indicates external purchase, and a negative value indicates power sold to the grid; π grid (t) represents the real-time electricity price at time t (unit: yuan / kWh), which changes dynamically over time (e.g., peak-valley price difference); P gas (t) represents the gas consumption of the system at time t (unit: m³). 3 / h), mainly from equipment such as gas turbines and gas boilers, is the basis for calculating gas procurement costs; π gas (t) represents the gas price at time t (unit: yuan / m³). 3 The duration of the control period may vary depending on the region and season; Δt is the control period duration (unit: h), for example Δt = 0.25 h (i.e. 15 min);

[0162] The equipment start-stop loss cost C is obtained by multiplying the number of times each controllable device changes its start-stop status within adjacent time periods by the preset loss cost coefficient for each start-stop operation. switch (t);

[0163]

[0164] In the formula, C switch (t) represents the equipment start-up and shutdown loss cost (unit: yuan); κ device The start-stop loss coefficient (unit: yuan / time) for a specific device is determined by the device type and model. For example, the start-stop loss coefficient for a gas turbine is much higher than that for a small energy storage device, and needs to be calibrated using the equipment manual or historical operation and maintenance data. device (t) represents the start / stop state of the device at time t, and is a binary variable. Typically, 1 indicates that the device is running, and 0 indicates that the device is stopped. ONOFF device (t-1) represents the start / stop state of the device at time t-1, used to compare with the current state to determine whether a change in start / stop status has occurred; ∑ device (*) is used to sum the start-up and shutdown loss costs of all controllable equipment (such as gas turbines, energy storage devices, and adjustable loads) in the system to ensure coverage of all equipment with high start-up and shutdown losses;

[0165] Multiplying the total carbon emissions generated during system operation by the unit carbon allowance price in the carbon emission trading market yields the carbon trading cost C. carbon (t);

[0166] C carbon (t)=γ·E carbon (t)

[0167] Among them, C carbon (t) represents the carbon trading cost (unit: yuan); γ is the unit carbon emission price (unit: yuan / tCO2), determined by the regional carbon trading market; E carbon (t) represents the total carbon emissions of the system at time t (unit: tCO2), calculated from the gas consumption of the gas turbine (e.g., a fixed mass of CO2 emitted per cubic meter of gas combustion);

[0168] Energy procurement cost C purchase (t), Equipment start-up and shutdown loss cost C switch (t), carbon trading cost C carbon The total operating cost of the system is obtained by adding the values ​​of (t) together, and the reward function R(t) is set to be the negative of the total operating cost of the system.

[0169] R(t) = -[C purchase (t)+C switch (t)+C carbon (t)]

[0170] Where R(t) is the reward function (unit: yuan), the larger its value, the lower the total operating cost of the system, thereby driving the agent to make decisions and regulate in the direction of economic optimization;

[0171] The training of the MADDPG prediction model includes:

[0172] The MADDPG network is pre-trained using an offline reinforcement learning algorithm. A quadruple containing state, action, reward, and next state is constructed using historical running data, and the MADDPG network is trained using this quadruple. The reward function R(t) is used as the benchmark for the target Q value of the Critic network. The network parameters of the MADDPG network are optimized by minimizing the temporal difference error of the Critic network and simultaneously updating the policy gradient of the Actor network. The network parameters are the connection weights and neuron biases in the Actor network and the Critic network.

[0173] Specifically, the loss function (i.e., time-series difference error) of the Critic network is:

[0174] L critic =E[(Q(s(t),a(t))-y(t)) 2 ]

[0175] In the formula, y(t)=R(t)+γ discount ·Q target (s(t+1),a′(t+1)) represents the objective Q-value, which is the optimization objective of the Critic network, combining the current reward R(t) and the expected future reward; γ discountThis is a discount factor (ranging from 0 to 1), used to weigh the importance of current rewards against future rewards. A value closer to 1 indicates a greater emphasis on long-term future returns; a typical value in industrial scenarios is 0.9 to 0.99. target (*) represents the target Critic network; Q target (s(t+1),a′(t+1)) represents the future action value (i.e., the target Q-value) output by the target Critic network, indicating the expected long-term reward for performing the action at time t+1; s(t+1) is the global state vector of the system at time t+1; a′(t+1) is the action generated by the target Actor network based on the next state, which is the predicted value of the optimal future action; L critic is the loss value of the Critic network. The smaller the value, the more accurate the Critic's evaluation of the action value, and it is the core basis for updating the Critic network parameters; E(*) is the mathematical expectation operator; Q(s(t),a(t)) is the action value (Q value) output by the current Critic network, representing the expected long-term reward for performing the action in the state at time t; s(t) is the global state vector of the system at time t, containing the local observation variables of all agents (such as gas turbine speed, energy storage SOC, load demand); a(t) is the joint action space of all agents;

[0176] After the MADDPG prediction model is trained, in the real-time decision-making stage, the standardized multi-source time series dataset {X(t)} and the probability distribution scene set {Y} generated in step S2 are used. scenatio (k)} Input into the pre-trained MADDPG network, for each probability distribution scenario k, the joint action space of each agent is generated through the Actor network, and the Q value of each action is evaluated through the Critic network. The action with the highest Q value is selected as the initial action a in that scenario. k (t);

[0177] For the initial action a in all probability distribution scenarios k (t) Perform expected weighting to obtain the preliminary control actions a of each controllable device. preliminary (t);

[0178]

[0179] Where prob(k) is the probability of scenario k occurring, generated by a conditional generative adversarial network, and the sum of the probabilities is 1, reflecting the likelihood of the scenario occurring; a k (t) represents the initial action in scenario k; a preliminary (t) represents the initial control actions of each controllable unit at time t, including the initial control commands of all intelligent agents, such as the target power of the gas turbine, the charging and discharging power of the energy storage device, the power adjustment amount of the adjustable load, and the opening commands of various valves.

[0180] Step S4: Based on the standardized multi-source time-series dataset, a pre-trained hybrid digital twin model is used to predict the future system state vector of the industrial integrated energy system. The hybrid digital twin model couples the outputs of a physical mechanism module constructed using differential-algebraic equations and a data-driven module constructed using a temporal convolutional network. Specifically, it includes:

[0181] Based on the variable types in the standardized multi-source time series dataset {X(t)}, the system modeling objects are divided into physical mechanism-dominated units and data-driven dominant units; among them, physical mechanism-dominated units include gas turbines, steam turbines and heat exchangers, and data-driven dominant units include photovoltaic power generation units, wind power generation units and complex fluid pipelines.

[0182] During the system simulation, for each simulation time step, the physical mechanism module is first run to calculate the state and output of the gas turbine, steam turbine and heat exchanger. Then, these outputs are used as inputs to the data-driven module to further calculate the dynamic behavior of the photovoltaic power generation unit, wind power generation unit and complex fluid pipeline network, thereby realizing the coupling of the physical mechanism module and the data-driven module and constructing a hybrid digital twin model.

[0183] The working steps of the physical mechanism module are as follows:

[0184] Based on the standardized multi-source time-series dataset {X(t)}, the current state variable X of the physical mechanism-dominant unit is obtained. phys And extract the control variable U related to the operation of the physical mechanism-dominant unit. phys Among them, the main physical mechanism units include gas turbines, steam turbines, and heat exchangers;

[0185] Among them, X phys State variables representing the dominant physical mechanism unit include, for a gas turbine, rotor speed, turbine inlet temperature, and exhaust temperature; for a steam turbine, rotor speed, steam pressure at each stage, and steam temperature at each stage; and for a heat exchanger, shell-side fluid temperature distribution and metal wall temperature. phys The control variables representing the dominant unit of the physical mechanism include, for gas turbines, fuel valve opening and guide vane angle; for steam turbines, inlet regulating valve opening and extraction valve opening; and for heat exchangers, cold / hot medium flow regulating valve opening.

[0186] Based on the current state variables and control variables, the set of physical parameters Θ of the physical mechanism-dominant unit is introduced. phys By solving the system of differential-algebraic equations describing the dynamic process and static relationship of the dominant physical mechanism unit, the state variable X of the dominant physical mechanism unit at the next time step is obtained. phys (t+1);

[0187] The differential equation part of the system of differential algebraic equations is expressed as:

[0188]

[0189] In the formula, X phys U represents a state variable; phys Indicates a control variable; Θ phys This refers to the set of physical parameters representing the dominant unit of the physical mechanism. For gas turbines, this includes heat capacity ratio, mechanical efficiency, and combustion efficiency; for steam turbines, it includes isentropic efficiency, mechanical loss coefficient, and rated capacity; and for heat exchangers, it includes heat transfer coefficient, heat transfer area, and flow resistance coefficient. phys (*) represents the state equation function describing the dynamic changes of the state variable; the updated state variable X is obtained by solving this equation using numerical integration. phys (t+1);

[0190] Specifically, taking a gas turbine as an example, by defining its state variables, control variables, and output variables, differential-algebraic equations are used to describe its dynamic processes and static relationships. Its physical mechanism module is described by the following set of equations:

[0191]

[0192] 0 = g gt (X gt ,Y gt ,Θ gt )

[0193] In the formula, X gt This is the state variable vector of the gas turbine, including rotor speed, turbine inlet temperature, and exhaust temperature. These variables are core indicators describing the real-time operating state of the gas turbine; U gt The control variable vector includes fuel valve opening and guide vane angle; adjusting these variables can change the operating state and output of the gas turbine. gt The output variable vector includes output electrical power, exhaust flow rate, and exhaust enthalpy, which are key parameters for measuring the operating performance and external energy supply capacity of a gas turbine; Θ gt This is a set of physical parameters for a gas turbine, including heat capacity ratio, mechanical efficiency, and combustion efficiency. These parameters are determined by the gas turbine's inherent hardware characteristics and serve as the fundamental input for the physical mechanism model; f gt (*) represents the differential equation for the state variables, describing a dynamic process, such as how the rotor speed changes with the adjustment of the fuel valve opening; g gt (*) represents an algebraic constraint equation that describes static relationships, such as the fixed calculation relationship between exhaust enthalpy and exhaust temperature and pressure.

[0194] Based on the state variable X at the next time step phys (t+1) and the set of physical parameters are used to calculate the algebraic constraint equations in the system of differential-algebraic equations, thereby obtaining the output variable Y of the physical mechanism-dominant unit at the current moment. phys ;

[0195] The algebraic constraint equation is expressed as:

[0196] 0 = g phys (X phys (t+1),Y phys ,Θ phys )

[0197] In the formula, X phys (t+1) represents the state variable at the next time step; g phys (*) represents an algebraic equation function describing the static relationship between variables; Y phys The output variables of the physical mechanism dominant unit are: for gas turbines, they include output electric power, exhaust flow rate and exhaust enthalpy; for steam turbines, they include output electric power, extraction steam flow rate and extraction steam enthalpy; and for heat exchangers, they include outlet temperature, heat exchange capacity and pressure drop.

[0198] The workflow of the data-driven module is as follows:

[0199] Based on the standardized multi-source time series dataset {X(t)}, historical multi-dimensional time series feature data of the data-driven dominant unit are extracted and fused with the output variable Y provided by the physical mechanism module. phys As internal boundary conditions of the system, the input data S constitutes the data-driven module. in Among them, the data-driven main units include photovoltaic power generation units, wind power generation units, and complex fluid pipeline networks;

[0200] Among them, S in The input data for the data-driven model includes, for photovoltaic power generation units, solar irradiance, ambient temperature, cloud cover coefficient, photovoltaic panel surface temperature, historical photovoltaic output data, and exhaust flow rate and exhaust enthalpy provided by the gas turbine model; for wind power generation units, wind speed, wind direction, air density, ambient temperature, atmospheric pressure, historical wind power output data, and extraction steam flow rate and extraction steam enthalpy provided by the steam turbine model; and for complex fluid networks, pipeline inlet pressure, pipeline inlet temperature, medium flow rate, ambient temperature, network topology data, historical load data, and outlet temperature and heat exchange rate provided by the heat exchanger model.

[0201] Input data S based on the data-driven module in Feature extraction is performed using dilated convolutional layers of a temporal convolutional network to obtain a high-dimensional feature sequence C. out ;

[0202] The calculation process of dilated convolutional layers is represented as follows:

[0203] C out =σ(W conv,d ·S in +b conv )

[0204] Among them, W conv,d b represents the weight matrix of the convolutional layer. conv The vector represents the bias vector of the convolutional layer; d represents the dilation factor of the convolution operation, which expands the receptive field by adjusting the dilation factor d (which takes values ​​of 1, 2, 4, 8, ...); σ(*) represents the ReLU activation function; C out This represents the high-dimensional feature sequence of the output;

[0205] Based on high-dimensional feature sequence C out By performing nonlinear transformation and feature fusion through the residual connection layer of a temporal convolutional network, a deep temporal feature representation is obtained.

[0206] The calculation process of the residual connection layer is expressed as follows:

[0207] Output residual =C out +F(C out )

[0208] Among them, C out F(C) represents the input feature sequence from the dilated convolutional layer; out () represents the nonlinear transformation performed by subsequent convolutional layers, weight normalization layers, and activation functions within the same residual block; Output residual This represents the deep temporal feature representation after feature fusion.

[0209] Based on deep temporal feature representation, a linear mapping is performed through a fully connected layer to obtain the predicted value P of the data-driven dominant unit at the next time step. out (t+1);

[0210] The computation process of a fully connected layer is represented as follows:

[0211] P out (t+1)=W fc ·H final +b fc

[0212] In the formula, H final W represents the final feature vector obtained after processing through multiple dilated convolutional layers and residual connection layers. fc Let b represent the weight matrix of the fully connected layer. fc P represents the bias vector of the fully connected layer.out (t+1) represents the predicted value of the data-driven dominant unit at time t+1. For photovoltaic power generation unit, it represents the predicted value of photovoltaic output at the next time. For wind power generation unit, it represents the predicted value of wind power output at the next time. For complex fluid pipeline network, it represents the predicted value of pipeline pressure distribution at the next time.

[0213] Specifically, taking a photovoltaic power generation unit as an example, relevant time-series data (such as solar irradiance, ambient temperature, and cloud coefficient) are extracted from the standardized multi-source time-series dataset {X(t)} and input together with the output variables of the physical mechanism-driven unit into a time-series convolutional network. Through the network's dilated convolution and residual connection structure, the predicted photovoltaic output value is output. The expression of its data-driven module is as follows:

[0214] P pv (t)=TCN(W pv ·[X solar (t),T amb (t),H cloud (t)]+b pv )

[0215] In the formula, P pv (t) represents the predicted photovoltaic power output at time t; X solar (t) represents the time-series data of solar irradiance at time t, which is a key input variable affecting photovoltaic (PV) output. Higher solar irradiance generally corresponds to higher PV output. amb (t) represents ambient temperature, a crucial factor affecting photovoltaic (PV) output; excessively high or low temperatures will impact the power generation efficiency of PV modules; H cloud (t) is the cloud cover coefficient, used to quantify the degree to which clouds reduce solar irradiance; W pv b represents the weights of a temporal convolutional network. pv The biases of the temporal convolutional network are determined through model training; TCN(*) represents the forward computation process of the temporal convolutional network, which uses dilated convolution and residual connection structures to capture long-term dependencies.

[0216] The output of the hybrid digital twin module is:

[0217] A hybrid digital twin model is constructed by coupling the data-driven module and the physical mechanism module through a state integration function; the predicted value P output by the data-driven module is then used to construct the model. out (t+1) and the updated state variable X from the physical mechanism model phys Integrating at (t+1) yields the system state vector X of the hybrid digital twin model at the next time step. sys (t+1);

[0218] The state integration function is expressed as:

[0219] X sys (t+1)=F hybrid (X sys (t),U sys (t),Θ phy ,Θ nn )

[0220] In the formula, X sys (t+1) represents the system state vector of the hybrid digital twin model at time t+1, which is the output of the simulation calculation and reflects the new state of the system after one time step; X sys (t) represents the system state vector of the hybrid digital twin model at time t, including the state variables of all units (such as gas turbines, photovoltaic power generation units, etc.), serving as the input for the current time step; U sys (t) represents the system control vector, such as adjustment commands for gas turbine control variables and photovoltaic power generation unit operating parameters, used to guide the system state to change in the target direction; Θ phy This is the set of physical parameters (such as heat capacity ratio, mechanical efficiency, etc.) of all the dominant physical mechanisms; Θ nn F is the set of network parameters for all data-driven dominant units (such as the weights and biases of a temporal convolutional network). hybrid (*) is the state integration function, which enables collaborative solving between the physical mechanism module and the data-driven module;

[0221] It should be noted that the pre-training steps for the hybrid digital twin model are as follows:

[0222] The physical mechanism module is trained using the maximum likelihood estimation method to adjust the physical parameters. Specifically, it iteratively solves the problem using a quasi-Newton method to obtain the parameter estimates with the highest probability of occurrence in the historical data sequence. The data-driven module is trained using a backpropagation algorithm based on stochastic gradient descent to train the network parameters. Specifically, it uses the Adam optimizer to iteratively update the network weights to minimize the prediction error. The objective function is to minimize the mean squared error between the output of the hybrid digital twin model and the actual measured value.

[0223]

[0224] In the formula, L hybrid Y represents the objective function value of the hybrid digital twin model; the smaller this value, the higher the simulation accuracy of the model. meas (t) represents the measured output vector of the actual industrial integrated energy system at time t, obtained by field sensors and other equipment, and serves as the benchmark data for measuring the model's accuracy; Y sim (t) represents the simulation output vector of the hybrid digital twin model under the same conditions.

[0225] Step S5: Using the initial control actions of the controllable device as the initial control sequence, rolling optimization is performed with the goal of minimizing the total system operating cost. A control barrier function is introduced during the optimization process. The initial control sequence is then subjected to real-time safety verification using the system state vector output by the pre-trained hybrid digital twin model for each optimization step, generating the final control command sequence. Specifically, this includes:

[0226] In this embodiment, the rolling optimization time domain of model predictive control is pre-set to be T. mpc The control time domain is T ctrl Furthermore, the optimization time domain is larger than the control time domain;

[0227] Specifically, in the real-time control of an integrated industrial energy system, rolling optimization of the time domain T mpc Typically set to the next hour (corresponding to four 15-minute time steps), controlling the time domain T. ctrl Set to the next 15 minutes (corresponding to 1 time step) to adapt to the system's dynamic response speed and real-time computing requirements;

[0228] In each control cycle, the optimization objective of model predictive control is to achieve the following in the prediction time domain T: mpc To minimize the total system operating cost, construct the objective function:

[0229]

[0230] In the formula, C purchase (τ) represents the energy procurement cost at time τ, C switch (τ) represents the equipment start-up and shutdown loss cost, C carbon (τ) represents the carbon trading cost, and the calculation method for each cost item is consistent with the formula used for the reward function R(t) in step S4; τ is the time step index in the prediction time domain, and t is the start time of the current control cycle; J mpc The objective function value of the model predictive control (unit: yuan) is the total operating cost of the system in the predicted time domain. The smaller the value, the better the control sequence. It is the core minimization object of the optimization algorithm.

[0231] The hybrid digital twin model is used as the internal prediction model, and the actual operating state of the system reflected by the standardized multi-source time series dataset {X(t)} is the system state vector X at the current time t. sys (t), so as to control the time domain T ctrl Initial control action a preliminary (t) is the initial control sequence U seq By performing multi-step forward simulation using a hybrid digital twin model, the system's performance in the optimized time domain T is obtained. mpc The state prediction sequence X within pred ;

[0232] The state update formula for forward simulation is:

[0233] X sys (τ+1)=F hybrid (X sys (τ),u(τ),Θ phy ,Θ nn )

[0234] In the formula, U seq =[u(t),u(t+1),...,u(t+T)] ctrl -1)] is the initial control sequence, u(τ) is the time τ (τ=t,t+1,…,t+T) ctrl -1) is the control action, and the initial value is taken from the initial control action a. preliminary The components at the corresponding time in (t); X sys (τ+1) is the system state vector at time τ+1, which is the output of the simulation prediction and reflects the new state of the system after the control action at time τ; X sys (τ) is the system state vector at time τ, containing the state variables of all controllable devices (such as rotational speed, state of charge, etc.); Θ phy For the set of physical parameters, Θ nn This is the set of network parameters, consistent with the definition in step S2;

[0235] X pred =[X sys (t+1),X sys (t+2),...,X sys (t+T mpc )]

[0236] In the formula, X pred For the state prediction sequence, X sys (τ) represents time τ (τ=t+1,t+2,…,t+T) mpc The system state vector;

[0237] During the optimization process, system operation constraints are introduced; these constraints include:

[0238] The gas turbine output power constraint requires that the actual output power of the gas turbine must be kept between the minimum and maximum power allowed by the equipment;

[0239] The charging and discharging power constraint requirement for energy storage devices is that the charging and discharging power of energy storage devices must not exceed the limits specified in their technical specifications.

[0240] Energy storage device charge constraint requires that the state of charge of the energy storage device must be maintained within the safe operating range, and adjustable load power constraint requires that the power adjustment of the adjustable load does not exceed its preset adjustable range.

[0241] The voltage constraint requirement for power busbars stipulates that the voltage deviation of each busbar node in a power system must not exceed the maximum allowable range.

[0242] Pressure constraints on thermal pipelines require that the operating pressure of each pipeline in the thermal system must be maintained between the minimum allowable pressure and the maximum allowable pressure;

[0243] Specifically, the gas turbine output power constraint:

[0244]

[0245] In the formula, P gt (τ) represents the actual output power of the gas turbine at time τ, and is the system state vector X. sys The components of (τ); This is the minimum output power limit for a gas turbine. The maximum output power limit of a gas turbine is determined by the equipment hardware.

[0246] Energy storage device charge and discharge power constraints:

[0247]

[0248] In the formula, P ess (τ) represents the charging and discharging power of the energy storage device at time τ, and is the control action output of the energy storage device; The minimum charge and discharge power limit for the energy storage device (usually a negative value, such as -500kW indicating a maximum charging power of 500kW) is determined by the charging rate of the energy storage battery. The maximum charge and discharge power limit for energy storage devices (usually a positive value, such as 500kW indicating a maximum discharge power of 500kW) is determined by the battery discharge rate and the inverter capacity.

[0249] Energy storage device charge constraint:

[0250] SOC min ≤SOC(τ)≤SOC max

[0251] In the formula, SOC(τ) represents the state of charge of the energy storage device at time τ, and is the system state vector X. sys (τ) components; SOC min The lower limit of the state of charge (SOC) of energy storage (e.g., 20%) max This represents the upper limit of the energy storage state of charge (e.g., 90%).

[0252] Adjustable load power constraints:

[0253]

[0254] In the formula, L dr(τ) represents the actual power consumption of the adjustable load at time τ (such as the power of industrial air conditioners and water pumps), and is the control action output of the adjustable load; The minimum power limit for adjustable loads is determined by the load process requirements (such as the minimum power required for air conditioning to maintain the indoor base temperature). The maximum power limit for adjustable loads is determined by the rated power of the load equipment and the power distribution capacity (such as the rated power of a water pump motor).

[0255] Power bus voltage constraints:

[0256] |U bus (τ)-U nom |≤ΔU max

[0257] In the formula, U bus (τ) represents the actual voltage of the power bus at time τ (e.g., 380V, 10kV), and is the system state vector X. sys (τ) is a component of; U nom The rated voltage of the power busbar (e.g., 380V for low-voltage busbars, 10kV for medium-voltage busbars) is the standard voltage for system design; ΔU max The maximum permissible voltage deviation (e.g., ±5% of rated voltage) is specified by power system specifications (e.g., GB / T 12325-2008);

[0258] Pressure constraints on thermal pipelines:

[0259]

[0260] In the formula, p pipe (τ) represents the actual pressure inside the thermal pipe at time τ (e.g., 0.8 MPa, 1.6 MPa), and is the system state vector X. sys The components of (τ); The minimum pressure limit for heat pipelines is determined by the heat delivery requirements (e.g., if the pressure is too low, hot water / steam cannot be delivered to distant users). The maximum pressure limit for heating pipelines is determined by the pipeline material and wall thickness (e.g., the maximum withstand pressure of carbon steel pipelines).

[0261] The sequential quadratic programming algorithm is used to solve the objective function optimization problem with system operation constraints, and the optimized control sequence U is obtained. opt During the optimization process, each iteration calls the pre-trained hybrid digital twin model to predict the state and calculates the gradient between the objective function and the system operation constraints.

[0262] U opt =[u opt (t),u opt(t+1),...,u opt (t+T ctrl -1)]

[0263] In the formula, U opt For the optimized control sequence, u opt (τ) represents time τ (τ=t,t+1,…,t+T) ctrl -1) Optimized control actions;

[0264] Based on the optimized control sequence U opt Extract the first control action u to be executed in the current control cycle. candidate , as candidate control instructions to be verified;

[0265] u candidate =u opt (t)

[0266] Among them, u candidate Candidate control commands;

[0267] For the key safety variables in the industrial integrated energy system, control barrier functions are defined respectively based on the above system operation constraints; among them, the key safety variables include gas turbine output power, energy storage device state of charge, power bus voltage and thermal pipeline pressure;

[0268]

[0269] h soc (SOC(τ))=(SOC(τ)-SOC min (SOC) max -SOC(τ))

[0270] h bus (U bus (τ))=(U bus (τ)-(U nom -ΔU max ))((U nom +ΔU max )-U bus (τ))

[0271]

[0272] In the formula, h gt (*) represents the control barrier function for the gas turbine output power, h soc (*) represents the control barrier function for the state of charge of the energy storage device, h bus (*) represents the control barrier function for the power bus voltage, h pipe (*) represents the control barrier function for the pressure in the thermal pipeline;

[0273] Substitute the system state vector at the current time t into each control barrier function to calculate the function value of each critical safety variable at the current time t, including the function value h of the gas turbine output power at the current time. gt (P gt (t) and the function value h of the current state of charge of the energy storage device. soc (SOC(t)) and the function value h of the current power bus voltage. bus (U bus (t) and the function value h of the current thermal pipeline pressure. pipe (p pipe (t));

[0274] Single-step forward simulation using a hybrid digital twin model to predict and execute candidate control commands u candidate The system state vector X of the system at the next time step. sys (t+1), and calculate the function value of each critical safety variable at the next time step by controlling the barrier function;

[0275] The state update formula for forward simulation is:

[0276] X sys (t+1)=F hybrid (X sys (t),u candidate ,Θ phy ,Θ nn )

[0277] In the formula, the state update formula is the same as the state update formula for the forward simulation described above;

[0278] If the function value of each key safety variable in the next time step is not lower than a safety threshold calculated from the function value at the current time step and the corresponding relaxation coefficient, then the current candidate control instruction is determined to be safe and can be directly used as the final instruction; otherwise, it is determined to be unsafe and needs to proceed to the next step for correction.

[0279] Safety conditions include the following inequalities:

[0280] h gt (P gt (t))≥-γ gt ·h gt (P gt (t+1))

[0281] h soc (SOC(t+1))≥-γ soc ·h soc (SOC(t))

[0282] h bus (U bus(t))≥-γ bus ·h bus (U bus (t+1))

[0283] h pipe (p plpe (t))≥-γ pipe ·h pipe (p pipe (t+1))

[0284] In the formula, γ gt γ is the relaxation factor for the output power of the gas turbine. soc γ is the relaxation factor for the state of charge of the energy storage device. bus γ is the relaxation factor for the power bus voltage. pipe The relaxation factor for the pressure in the thermal pipeline is typically taken as 0.1 to 0.3; h gt (P gt (t+1)) is the function value of the gas turbine output power at the next moment, h soc (SOC(t+1)) is the function value of the energy storage device's state of charge at the next moment, h bus (U bus (t+1)) is the function value of the power bus voltage at the next moment, h pipe (p pipe (t+1)) is the function value of the thermal pipeline pressure at the next moment; if the safety inequality is also satisfied, then u final (t)=u candidate As the final control command, otherwise the following corrections are performed;

[0285] When a candidate control command is determined to be unsafe, the candidate control command u is evaluated using a sequential quadratic programming algorithm. candidate Minimum disturbance correction is performed to construct and solve a quadratic programming problem with the goal of minimizing the control command correction amount and the constraint of satisfying the safety condition of the control barrier function. The corrected control command u obtained from the solution is then used. corrected The above security check is re-executed as input for state update, resulting in the corrected security control command u. safe ;

[0286] The optimization objective of the quadratic programming problem is:

[0287] min‖u corrected -u candidate || 2

[0288] In the formula, u corrected The corrected control command;

[0289] X s ′ys (t+1)=F hybrid (X sys (t),u corrected ,Θ phy ,Θ nn )

[0290] In the formula, X′ sys (t+1) is the corrected system state vector, including the function value h′ of the gas turbine output power at the next corrected moment. gt (P gt (t+1));

[0291] The constraints of the quadratic programming problem are inequalities that simultaneously satisfy safety conditions, where the function value (e.g., h′) of each critical safety variable in the safety condition inequality is the value of that variable at the next time step. gt (P gt (t+1)) is updated to the corrected system state vector X′ sys (t+1) is calculated;

[0292] The corrected safety control command u safe As the final control command at the current moment u final (t);

[0293] For the optimized control sequence U opt Each control action u in opt (τ)(τ=t,t+1,…,t+T ctrl -1) Perform real-time safety checks and necessary corrections sequentially, and then convert the final control command u final (t) are sorted sequentially to obtain the final control command sequence U. final And distribute it to the physical system for execution;

[0294] U final =[u final (t),u final (t+1),...,u final (t+T ctrl -1)]

[0295] Among them, U final For the final control instruction sequence, t final (τ) represents time τ (τ=t,t+1,…,t+T) ctrl -1) final control command.

[0296] This embodiment also provides a self-updating mechanism for model parameters. By monitoring reconstruction errors through an autoencoder, incremental model updates are triggered. Compared to traditional control modes that rely heavily on manual intervention and lack dynamic optimization and safety loops, this mechanism achieves economic optimization and safety assurance of control actions. Furthermore, it can adaptively update the model to reflect changes in system state, solving the problems of traditional control's delayed response, susceptibility to safety constraints, and difficulty in long-term adaptation to system dynamics. This improves the system's autonomous operation capability and stability. In simple terms, after obtaining the final control command sequence through step S5, the final control command sequence is sent to the physical system for execution, and system operation data is collected simultaneously. The autoencoder calculates the reconstruction error. When the reconstruction error continuously exceeds a preset threshold, parameter updates are triggered in the data-driven module and the MADDPG prediction model within the hybrid digital twin model. Specifically, this includes:

[0297] In this embodiment, the training of the autoencoder network includes:

[0298] An autoencoder network consisting of an encoder and a decoder is constructed. The encoder compresses the standardized multi-source time-series data of historical moments from D dimensions to a d-dimensional latent space, and the decoder obtains the reconstructed data of historical moments.

[0299] Specifically, the encoder consists of two fully connected layers. The first layer has an input dimension of D (i.e., the dimension of the standardized multi-source time series dataset) and an output dimension of... The activation function is ReLU; the input dimension of the second layer is... The output dimension is d (the dimension of the latent space, The activation function is linear; the decoder is symmetric to the encoder, the first layer has an input dimension of d and an output dimension of d. The activation function is ReLU; the input dimension of the second layer is... The output dimension is D, and the activation function is Sigmoid;

[0300] Based on standardized multi-source time-series data from historical moments, an autoencoder network is unsupervised pre-trained by minimizing the mean squared error between the input and reconstructed data. During training, the Adam optimizer is used with a learning rate of r (e.g., 0.001) and a batch size of Batch. size (e.g., 128), training period is T. train (e.g., 200); After training is complete, save the network parameters of the encoder and decoder;

[0301] In real-time operation, the standardized multi-source time-series data X(t) at the current moment is input into the autoencoder to calculate the reconstruction error E. rec (t);

[0302]

[0303] In the formula, X is the reconstructed data output by the autoencoder network at time t and dimension i; i (t) represents standardized multi-source time-series data at time t and dimension i; E rec (t) represents the reconstruction error at time t;

[0304] The error threshold E is determined by the statistical distribution of reconstruction errors on the autoencoder based on historical normal data. threshold (Take the upper 95th percentile of the distribution);

[0305] Monitor the reconstruction error within a continuous time window, if the continuous θ n Reconstruction error E per control cycle rec (t) all exceed the error threshold E threshold If the system's operating state deviates significantly from the expectations of the hybrid digital twin model, a model update flag is triggered; otherwise, the system continues to operate according to the current model and control strategy.

[0306] Where, θ n For a preset periodic threshold, such as θ n =5, which corresponds to approximately 1.25 hours in an industrial setting, assuming a control cycle of 15 minutes;

[0307] When the model update flag is triggered, the incremental learning process is initiated, collecting system operation data over a recent period (e.g., the past 24 hours), including the latest standardized multi-source time-series dataset X. new and the corresponding actual measured value

[0308] Extract the latest standardized multi-source time series dataset X new The input features related to the data-driven dominant unit (such as solar irradiance, ambient temperature, etc.) and their corresponding actual outputs (such as photovoltaic output, wind power output) are used for incremental training with mini-batch gradient descent to update the parameters of the temporal convolutional network and minimize the output of the hybrid digital twin model. Compared with actual measured value Mean square error between:

[0309]

[0310] In the formula, The actual measured values ​​of the data-driven unit at time t, such as the measured values ​​of photovoltaic power output and the instrument readings of wind power output, come directly from the system sensors; L represents the simulation output of the data-driven component at time t, i.e., the output value predicted by the model under the same input (such as solar irradiance and wind speed), which is the output result to be optimized; nnThe mean squared error loss value is the target to be minimized during incremental training, covering the entire time range of the incremental training dataset (e.g., the past 24 hours).

[0311] The learning rate in the incremental training phase is set to 1 / 10 of that in the pre-training phase (e.g., 0.0001), the batch size is 64, and the training cycle is 50.

[0312] Construct a recent experience replay buffer to store the above quadruplets, and use this buffer to fine-tune the policy of the MADDPG network using offline reinforcement learning, and update the parameters of the Actor network and Critic network.

[0313] After completing incremental learning and parameter updates, the original model is replaced with the updated hybrid digital twin model and MADDPG prediction module. The model update flag is cleared, the system continues to run, and returns to step S1 to start the next round of data acquisition and control cycle.

[0314] Example 2

[0315] Figure 2 This is a schematic diagram of a real-time decision-making device for an integrated industrial energy system according to another exemplary embodiment, the device comprising:

[0316] Data acquisition module 1: used to acquire real-time operating data of the industrial integrated energy system, and preprocess the real-time operating data to obtain a standardized multi-source time series dataset;

[0317] Probability distribution scenario prediction module 2: Based on the standardized multi-source time series dataset, it uses a pre-trained probability distribution scenario prediction model to generate probability distribution scenario sets of renewable energy output and multi-load at different prediction time scales;

[0318] Preliminary control action generation module 3: Based on the standardized multi-source time series dataset and probability distribution scenario set, it uses a pre-trained MADDPG prediction model to generate preliminary control actions for each controllable device in the industrial integrated energy system;

[0319] Future system state vector prediction module 4: Based on the standardized multi-source time series dataset, it uses a pre-trained hybrid digital twin model to predict the future system state vector of the industrial integrated energy system. The hybrid digital twin model is coupled with the output of a physical mechanism module constructed by differential algebraic equations and a data-driven module constructed by temporal convolutional networks.

[0320] Control command generation module 5: It is used to perform rolling optimization with the initial control action of the controllable device as the initial control sequence and the goal of minimizing the total operating cost of the system. During the optimization process, a control barrier function is introduced. The initial control sequence is verified in real time by the system state vector of each optimization step output by the pre-trained hybrid digital twin model, and the final control command sequence is generated.

[0321] Example 3:

[0322] This embodiment provides a storage medium storing a computer program, which, when executed by a host controller, implements the various steps in the above method.

[0323] It is understood that the storage medium mentioned above can be a read-only memory, a hard disk, or an optical disk, etc.

[0324] It is understood that the same or similar parts in the above embodiments can be referred to each other, and the contents not described in detail in some embodiments can be referred to the same or similar contents in other embodiments.

[0325] It should be noted that in the description of this invention, the terms "first," "second," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance. Furthermore, in the description of this invention, unless otherwise stated, "a plurality of" means at least two.

[0326] Any process or method description in the flowchart or otherwise herein can be understood as representing a module, segment, or portion of code comprising one or more executable instructions for implementing a particular logical function or process, and the scope of the preferred embodiments of the invention includes additional implementations in which functions may be performed not in the order shown or discussed, including substantially simultaneously or in reverse order depending on the functions involved, as will be understood by those skilled in the art to which embodiments of the invention pertain.

[0327] It should be understood that various parts of the present invention can be implemented in hardware, software, firmware, or a combination thereof. In the above embodiments, multiple steps or methods can be implemented in software or firmware stored in memory and executed by a suitable instruction execution system. For example, if implemented in hardware, as in another embodiment, it can be implemented using any one or a combination of the following techniques known in the art: discrete logic circuits having logic gates for implementing logical functions on data signals, application-specific integrated circuits (ASICs) having suitable combinational logic gates, programmable gate arrays (PGAs), field-programmable gate arrays (FPGAs), etc.

[0328] Those skilled in the art will understand that all or part of the steps of the methods in the above embodiments can be implemented by a program instructing related hardware. The program can be stored in a computer-readable storage medium, and when executed, the program includes one or a combination of the steps of the method embodiments.

[0329] Furthermore, the functional units in the various embodiments of the present invention can be integrated into a processing module, or each unit can exist physically separately, or two or more units can be integrated into a module. The integrated module can be implemented in hardware or as a software functional module. If the integrated module is implemented as a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium.

[0330] The storage media mentioned above can be read-only memory, disk, or optical disk, etc.

[0331] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0332] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. An industrial integrated energy system real-time decision-making method, characterized in that, The method comprises: obtaining real-time operation data of an industrial comprehensive energy system, and preprocessing the real-time operation data to obtain a standardized multi-source time series data set; based on the standardized multi-source time series data set, using a pre-trained probability distribution scenario prediction model to generate a probability distribution scenario set of renewable energy output and multi-element load at different prediction time scales; based on the standardized multi-source time series data set and the probability distribution scenario set, using a pre-trained MADDPG prediction model to generate preliminary control actions of each controllable device in the industrial comprehensive energy system; based on the standardized multi-source time series data set, using a pre-trained hybrid digital twin model to predict a future system state vector of the industrial comprehensive energy system, wherein the hybrid digital twin model couples the output of a physical mechanism module constructed by a differential algebraic equation and a data-driven module constructed by a time series convolution network; using the preliminary control actions of the controllable devices as an initial control sequence, performing rolling optimization with the goal of minimizing the total operation cost of the system, introducing a control barrier function in the optimization process, and performing real-time safety checking on the initial control sequence by the system state vector output by the pre-trained hybrid digital twin model at each optimization step to generate a final control instruction sequence.

2. The method of claim 1, wherein the training of the probability distribution scenario prediction model comprises: constructing a conditional generative adversarial network composed of a generator and a discriminator; setting different historical window lengths and future prediction steps for three different prediction time scales of ultra-short term, short term and day-ahead; obtaining historical standardized multi-source time series data sets and extracting condition information of the historical window length therefrom, taking the condition information and a random noise vector as inputs of the generator, and outputting a probability distribution scenario of the future prediction step; extracting real data from the historical standardized multi-source time series data set, and determining the probability of real data and generator prediction data being real data by the discriminator; calculating the difference between the generated prediction data and the real data by an adversarial loss function to obtain an adversarial loss value; training the conditional generative adversarial network using the historical standardized multi-source time series data set, alternately optimizing the generator and the discriminator during the training process, and introducing a physical consistency loss value to jointly optimize the adversarial loss value and the physical consistency loss value to minimize the total loss; and obtaining the trained probability distribution scenario prediction model; the use of the pre-trained probability distribution scenario prediction model to generate a probability distribution scenario set of renewable energy output and multi-element load at different prediction time scales comprises: for future ultra-short term, short term and day-ahead time scales, generating a number of probability distribution scenarios based on the different future prediction steps set by the generator of the pre-trained probability distribution scenario prediction model. The K-means clustering method is used to reduce the generated probability distribution scenarios at each time scale, the Euclidean distance between scenarios is used as the similarity measure, the minimum within-class sum of squares is used as the optimization objective, the center scenarios of each class are selected as the representative scenario subset, and the proportion of the number of initial scenarios contained in the class to the total number of scenarios is used as the occurrence probability to form the reduced probability distribution scenario set.

3. The method of claim 2, wherein, The training of the MADDPG prediction model comprises: Each controllable device in the industrial integrated energy system is regarded as an agent, an independent Actor network is established for each agent, and a centralized Critic network is established to guide the policy update of each Actor network, thereby constructing a MADDPG network; Local observation variables of each agent are extracted from a historical standardized multi-source time series data set; Based on the local observation variables of each agent, the Actor network receives the local observation variables of a single agent, outputs the continuous control action of the agent, and forms a joint action space of all agents; The energy procurement cost, the device start-stop loss cost, and the carbon trading cost are obtained respectively, and the total operating cost of the system is calculated; a reward function is set as the negative value of the total operating cost of the system; An offline reinforcement learning algorithm is used to pre-train the MADDPG network, a four-tuple containing state, action, reward, and next state is constructed using historical operating data, the four-tuple is used to train the MADDPG network; the reward function is used as the benchmark of the target Q value of the Critic network, the policy gradient of the Actor network is updated synchronously by minimizing the time difference error of the Critic network, so as to optimize the network parameters of the MADDPG network, including the connection weights and neuron biases in the Actor network and the Critic network; and a pre-trained MADDPG prediction model is obtained; The preliminary control action of the controllable device in the industrial integrated energy system is generated by using the pre-trained MADDPG prediction model, which comprises: The standardized multi-source time series data set and the probability distribution scenario set are input into the pre-trained MADDPG prediction model, for each probability distribution scenario, the joint action space of each agent is generated through the Actor network, and the target Q value of each action is evaluated through the Critic network, and the action with the highest target Q value is selected as the preliminary action under the scenario; The preliminary actions under all probability distribution scenarios are weighted to obtain the preliminary control action of each controllable device.

4. The method of claim 3, wherein, The output of the physical mechanism module comprises: Based on the standardized multi-source time series data set, the current state variables of the physical mechanism dominant unit are obtained, and the control variables related to the operation of the physical mechanism dominant unit are extracted; Based on the current state variable and the control variable, a set of physical parameters of a physical mechanism dominant unit is introduced, and by solving a differential algebraic equation set describing the dynamic process and static relationship of the physical mechanism dominant unit, a state variable of the physical mechanism dominant unit at the next moment is obtained; Based on the state variable at the next moment and the set of physical parameters, an output variable at the current moment is output by calculating an algebraic constraint equation in the differential algebraic equation set.

5. The method of claim 4, wherein, the output of the data-driven module includes: based on the standardized multi-source time series data set, historical multi-dimensional time series feature data of a data-driven dominant unit is obtained, and an output variable at the current moment provided by the physical mechanism module is fused as input data of the data-driven module; the data-driven module extracts features of the input data through an expansion convolution layer of a time series convolution network to obtain a high-dimensional feature sequence; the data-driven module performs nonlinear transformation and feature fusion on the high-dimensional feature sequence through a residual connection layer of the time series convolution network to obtain a deep time series feature representation; the data-driven module performs linear mapping on the deep time series feature representation through a fully connected layer to output a predicted value at the next moment; the use of the pre-trained hybrid digital twin model to predict the future system state vector of the industrial comprehensive energy system includes: the predicted value at the next moment output by the data-driven module is integrated with the state variable at the next moment obtained by the physical mechanism module to output the future system state vector of the industrial comprehensive energy system.

6. The method of claim 5, wherein, the training of the hybrid digital twin model includes: setting an objective function of the hybrid digital twin model, the objective function being to minimize the mean square error between the output of the hybrid digital twin model and the actual measured value; iteratively training the hybrid digital twin model through historical operation data, including: adjusting the physical parameters of the physical mechanism module using a maximum likelihood estimation method; adjusting the network parameters of the data-driven module using a back propagation algorithm based on a stochastic gradient descent; and obtaining a trained hybrid digital twin model when the objective function is satisfied through iterative training.

7. The method of claim 6, wherein, the use of the preliminary control action of the controllable device as an initial control sequence to perform rolling optimization with the goal of minimizing the total system operation cost includes: setting a rolling optimization time domain and a control time domain, the rolling optimization time domain being greater than the control time domain; in each control period, constructing an objective function with the goal of minimizing the total system operation cost within the rolling optimization time domain; using the actual system state reflected by the standardized multi-source time series data set as the system state vector at the current moment, using the preliminary control action within the control time domain as the initial control sequence, and performing multi-step forward simulation through the pre-trained hybrid digital twin model to obtain a state prediction sequence within the rolling optimization time domain. In the optimization process, system operation constraints are introduced; the system operation constraints include gas turbine output power constraints, energy storage device charging and discharging power constraints, energy storage device charge constraints, adjustable load power constraints, power bus voltage constraints, and heat pipe pressure constraints; A sequential quadratic programming algorithm is used to solve the optimization problem of the objective function with system operation constraints to obtain the optimized control sequence.

8. The method of claim 7, wherein, The control barrier function is introduced in the optimization process, and the system state vector output by the pre-trained hybrid digital twin model at each optimization step is used to perform real-time safety verification on the initial control sequence to generate a final control instruction sequence, which includes: Based on the optimized control sequence, the first control action to be executed in the current control period is extracted as a candidate control instruction to be verified; Key safety variables in the industrial comprehensive energy system are defined, and control barrier functions are defined based on the system operation constraints; The system state vector at the current time is substituted into each control barrier function to calculate the function value of each key safety variable at the current time; A single-step forward simulation is performed through the pre-trained hybrid digital twin model to predict the system state vector at the next time after executing the candidate control instruction, and the function value of each key safety variable at the next time is calculated through the control barrier function; If the function value of each key safety variable at the next time is not less than a safety threshold calculated based on the corresponding current time function value and the corresponding relaxation coefficient, it is determined that the current candidate control instruction is safe and can be directly used as the final instruction; otherwise, it is determined that the current candidate control instruction is unsafe, and the current candidate control instruction is modified to obtain a corrected safe control instruction; The corrected safe control instruction is used as the final control instruction at the current time; The safety of each control action in the optimized control sequence is determined and modified in sequence, and the final control instructions obtained are sorted in sequence to obtain a final control instruction sequence.

9. The method of claim 8, wherein, The current candidate control instruction is modified to obtain a corrected safe control instruction, which includes: Based on the sequential quadratic programming algorithm, the candidate control instruction is corrected with minimum disturbance to construct and solve a quadratic programming problem with the minimum control instruction correction amount as the optimization objective and the satisfaction of the control barrier function safety condition as the constraint to obtain the corrected safe control instruction.

10. The method of claim 3, wherein, Further comprising: The final control instruction sequence is sent to the physical system for execution, and system operation data is synchronously collected to calculate the reconstruction error through the pre-trained autoencoder; when the reconstruction error continuously exceeds the preset error threshold, the parameter update of the data-driven module and the MADDPG prediction model in the hybrid digital twin model is triggered.

11. The method of claim 10, wherein, The reconstruction error is calculated through the pre-trained autoencoder, which includes: An autoencoder network including an encoder and a decoder is constructed, the standardized multi-source time series data at a historical time is compressed from D dimensions to d dimensions of a latent space by the encoder, and the reconstructed data at the historical time is obtained by the decoder; Based on the historical standardized multi-source time series data, the autoencoder network is pre-trained unsupervisedly by minimizing the mean square error between the input data and the reconstructed data, and a trained autoencoder is obtained; The real-time standardized multi-source time series data is input into the trained autoencoder to obtain a reconstruction error; The reconstruction error distribution of the historical standardized multi-source time series data on the trained autoencoder is counted to obtain an error threshold; If the reconstruction errors of a continuous preset time period are all greater than the error threshold, the parameter update of the data-driven module and the MADDPG prediction model in the hybrid digital twin model is triggered; When the model parameter update is triggered, an incremental learning process is started, including: obtaining system operation data of a recent preset time period, the system operation data including a standardized multi-source time series data set and corresponding actual measurement values of the recent preset time period; The input features related to the data-driven module in the hybrid digital twin model and the corresponding actual output in the standardized multi-source time series data set of the recent preset time period are extracted, incremental training is performed by using a small batch gradient descent method, the parameters of the time series convolution network are updated, and the mean square error between the output of the hybrid digital twin model and the actual measurement value is minimized. An experience replay buffer for storing four-tuples is constructed, and the MADDPG network is updated in parameters by using offline reinforcement learning through the buffer.

12. An industrial integrated energy system real-time decision device, characterized in that, The device comprises: A data acquisition module is configured to acquire real-time operation data of an industrial comprehensive energy system, and to preprocess the real-time operation data to obtain a standardized multi-source time series data set; A probability distribution scenario prediction module is configured to generate a probability distribution scenario set of renewable energy output and multi-element load at different prediction time scales based on the standardized multi-source time series data set by using a pre-trained probability distribution scenario prediction model; A preliminary control action generation module is configured to generate preliminary control actions of controllable devices in the industrial comprehensive energy system based on the standardized multi-source time series data set and the probability distribution scenario set by using a pre-trained MADDPG prediction model; A future system state vector prediction module is configured to predict a future system state vector of the industrial comprehensive energy system based on the standardized multi-source time series data set by using a pre-trained hybrid digital twin model, wherein the hybrid digital twin model couples the outputs of a physical mechanism module constructed by a differential algebraic equation and a data-driven module constructed by a time series convolution network; A control instruction generation module is configured to take the preliminary control actions of the controllable devices as an initial control sequence, to perform rolling optimization with the objective of minimizing the total operation cost of the system, to introduce a control barrier function in the optimization process, to perform real-time safety checking on the initial control sequence by using the system state vector of each optimization step output by the pre-trained hybrid digital twin model, and to generate a final control instruction sequence.

13. A storage medium, characterized by The storage medium stores a computer program, and the computer program is executed by the processor to implement each step in the industrial comprehensive energy system real-time decision method according to any one of claims 1-11.

Images