Energy complex flexible dynamic aggregation method and system based on reinforcement learning

Abstract

The invention provides an energy complex flexible dynamic aggregation method and system based on reinforcement learning, and the method comprises the steps: collecting the real-time operation parameters of electric, gas, heat and multi-energy coupling equipment in an energy complex through an industrial Ethernet protocol, and constructing a continuous time dynamic system model through employing a Shenchang differential equation process modeling technology; and the nonlinear dynamic characteristics of various energy resources are accurately described. Modeling an aggregation process into a Markov decision process, and training and optimizing by adopting an improved twin delay depth deterministic strategy gradient algorithm to generate a reinforcement learning aggregation strategy model. The system realizes self-adaptive dynamic optimization through an online learning updating mechanism, and the renewable energy consumption rate and the system operation efficiency are remarkably improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of energy management, in particular to a dynamic and real-time aggregation method and system for energy complex based on reinforcement learning. BACKGROUND

[0002] As a complex system integrating multiple energy forms, energy complex has become an important development direction of modern energy management. With the increasing penetration of renewable energy and the diversification of energy demand, how to effectively aggregate and dispatch various flexible resources within the energy complex to achieve the coordination of economy, reliability and environmental protection has become a core technical challenge in the current energy field.

[0003] Existing energy aggregation technologies mainly include static aggregation methods based on optimization algorithms and dynamic aggregation methods based on predictive control. Static aggregation methods usually use linear programming or quadratic programming algorithms to determine the optimal combination of resources according to the preset objective function and constraint conditions, but they are difficult to adapt to real-time changes in system operating state. Dynamic aggregation methods based on predictive control can consider the changes in system state in future time periods, but their prediction accuracy is limited by the accuracy of the model, and the computational complexity is high, making it difficult to meet the requirements of real-time control.

[0004] In addition, so far, the energy aggregation method based on model predictive control (MPC) has been mainly used. This technology establishes a mathematical model of the energy system, predicts the system state in future time periods, and solves an optimization problem in each control period to determine the optimal control strategy. Its basic principle is to minimize the system operating cost or maximize the benefit within a limited time window, while satisfying the operating constraints of various devices and system balance constraints, and achieving dynamic scheduling through rolling optimization.

[0005] However, there are two key problems in the existing technology: first, it is difficult to accurately model the dynamic characteristics of complex multi-element heterogeneous energy systems, especially when facing non-linear, time-varying and uncertain factors, the accuracy and adaptability of traditional modeling methods are insufficient; second, there is a lack of effective processing capability for high-dimensional, non-uniform time series data, resulting in inaccurate evaluation of resource response characteristics and affecting the optimization effect of aggregation decisions. SUMMARY

[0006] The purpose of the present application is to overcome the shortcomings of the prior art and provide a dynamic and real-time aggregation method and system for energy complex based on reinforcement learning, which can accurately model the dynamic characteristics of complex multi-element heterogeneous energy systems, effectively process high-dimensional, non-uniform time series data, and realize dynamic aggregation and optimal dispatch of flexible resources within the energy complex.

[0007] To achieve the above-mentioned purpose, the technical solutions provided by the present application are as follows: The application discloses a flexible energy complex real-time aggregation method based on reinforcement learning, which comprises the following steps: Obtaining operation data of flexible resources in an energy complex, collecting real-time operation parameters of electricity, gas, heat and multi-energy coupling equipment, and using a neural ordinary differential equation process modeling technology to construct a continuous-time dynamic system model; Based on the continuous-time dynamic system model, a static pulse neural network is used to preliminarily identify an energy aggregation problem, and after identifying a potential aggregation demand, an adaptive dynamic pulse neural network classifier is activated to accurately classify a specific energy aggregation problem; Obtaining different aggregation demand types in the accurate classification result of the specific energy aggregation problem, using a multi-marginal stochastic flow matching method to align and analyze energy data measured at non-equidistant time points, and quantitatively evaluating the response characteristics of each flexible resource through a measured value spline enhancement technology and fractional matching, wherein the response characteristics of each flexible resource include the randomness of state transition and the sequence of decision-making in the energy aggregation process; Modeling the energy aggregation process as a Markov decision process, and using a twin-delayed deep deterministic policy gradient algorithm for training and optimization to generate a reinforcement learning aggregation strategy optimization model; Using the reinforcement learning aggregation strategy optimization model to receive state data in real time and output an optimal aggregation strategy, using the optimal aggregation strategy to send adjustment instructions to each flexible resource through a distributed control system, and collecting actual response data for online learning and updating to realize adaptive dynamic optimization.

[0008] Preferably, obtaining operation data of flexible resources in an energy complex, collecting real-time operation parameters of electricity, gas, heat and multi-energy coupling equipment, and using a neural ordinary differential equation process modeling technology to construct a continuous-time dynamic system model, comprises the following steps: Obtaining operation data of each flexible resource, removing data points beyond the normal range by using a three-sigma rule and eliminating dimension differences by Z-score standardization to generate preprocessed operation data; Based on the preprocessed operation data, the state evolution process of an energy system is represented in the form of a differential equation, a neural network with a residual network architecture is used to approximate a differential equation function and a gradient is calculated by using an adjoint sensitivity method, and a battery energy storage model, a gas turbine model and a heat energy storage model are established; A variable step ordinary differential equation solver is used to train the battery energy storage model, the gas turbine model and the heat energy storage model respectively, and the continuous-time dynamic system model is obtained.

[0009] Preferably, based on the continuous-time dynamic system model, a static spiking neural network is employed to preliminarily identify the energy aggregation problem, and after identifying the potential aggregation demand, an adaptive dynamic spiking neural network classifier is activated to accurately classify the specific energy aggregation problem, including: After obtaining the normalized state data, the static spiking neural network employs a leaky integrate-and-fire neuron model for pulse coding to identify four types of aggregation demands, including load peak-valley regulation, frequency regulation, voltage support, and emergency response; Based on the four types of aggregation demands, an adaptive dynamic spiking neural network classifier is activated, and the network structure is dynamically adjusted according to the novelty of the input data through the adaptive dynamic spiking neural network classifier, automatically increasing neurons to expand the network capacity; The connection weights of the dynamic spiking neural network classifier are updated through the adaptive spike-timing-dependent plasticity rule, and the specific energy aggregation problem is accurately classified through the dynamic spiking neural network classifier.

[0010] Preferably, based on the four types of aggregation demands, an adaptive dynamic spiking neural network classifier is activated, and the network structure is dynamically adjusted according to the novelty of the input data through the adaptive dynamic spiking neural network classifier, automatically increasing neurons to expand the network capacity, including: The feature vector of the four types of aggregation demands is obtained, the Euclidean distance between the input data and the existing neuron prototype vector is calculated, and when the minimum distance exceeds the preset novelty threshold, it is determined as a new mode, triggering the network structure adjustment mechanism; Based on the network structure adjustment mechanism, a new hidden layer neuron is created in the corresponding output category, the feature vector of the current input data is taken as the prototype vector of the new neuron, and the connection weights of the new neuron with the input layer and the output layer are initialized; The activation strength of the new neuron in the network is updated through the competitive learning mechanism to determine the optimal response neuron, and finally the expanded network capacity is determined.

[0011] Preferably, a multi-marginal random flow matching method is used to align and analyze the energy data measured at non-equidistant time points, and through measurement value spline enhancement technology and fractional matching, the response characteristics of each flexible resource are quantitatively evaluated, including: The data distribution characteristics of different types of energy equipment are obtained, and multiple marginal distributions of electric energy storage resources, thermal energy storage resources, and gas energy storage resources are constructed, and non-parametric modeling is performed through kernel density estimation method to generate a smoothed distribution model; Based on the smoothed distribution model, a cubic B-spline basis function is used to process the observation values at irregular time points through measurement value spline technology to generate a continuous time series; The response characteristics of the flexible resources are obtained by calculating the difference between the minimized real data distribution of the continuous time sequence and the fractional function of the distribution model using a fractional matching algorithm.

[0012] Preferably, the response characteristics of the flexible resources are obtained by calculating the difference between the minimized real data distribution of the continuous time sequence and the fractional function of the distribution model using a fractional matching algorithm, including: The real data distribution of the continuous time sequence is obtained, the logarithmic density gradient of the real data distribution is calculated and taken as the real fractional function, the logarithmic density gradient of the distribution model is calculated and taken as the fractional function of the distribution model; Based on the real fractional function and the fractional function of the distribution model, a fractional matching objective function is constructed, the least square method is used to calculate the square error between the real fractional function and the fractional function of the distribution model, and the calculation result is integrated and summed; The model parameters of the fractional matching objective function are optimized by a gradient descent algorithm to make the fractional function of the distribution model approach the real fractional function, and the response characteristics of the flexible resources including the response time constant, the adjustment accuracy and the stability index are obtained.

[0013] Preferably, based on the smoothed distribution model, a measured value spline technology is used to construct a cubic B-spline basis function to process the observation values of irregular time points, and a continuous time sequence is generated, including: The observation values of irregular time points in the smoothed distribution model are obtained, the node sequence of the spline is determined, the node density is set according to the time interval distribution of the observation data, and the number of nodes is increased in the data-intensive area; Based on the node sequence, a cubic B-spline basis function is constructed, each B-spline basis function is a cubic polynomial in four consecutive node intervals and zero in other intervals; The coefficients of the cubic B-spline basis function are determined by the least square fitting method to make the spline curve pass through or approach all observation points, and the continuous time sequence with time continuity and smoothness is obtained.

[0014] Preferably, a twin delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregated policy optimization model, including: The energy aggregation process is modeled as a Markov decision process, a state space including load demand, renewable energy output and real-time state of resources, and an action space including resource combination and output allocation ratio are constructed, a reward function integrating response speed, regulation accuracy and operation cost is designed, and based on the Markov decision process framework, state space, action space and reward function, a twin-delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregation strategy optimization model.

[0015] Preferably, based on the Markov decision process framework, state space, action space and reward function, a twin-delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregation strategy optimization model, comprising: Multi-dimensional characteristic data is acquired, a state space including load demand dimension, renewable energy dimension, energy storage state dimension, power quality dimension and economic signal dimension, and an action space including resource participation degree allocation vector and aggregation mode selection are constructed, and a Markov decision process framework is established; In the Markov decision process framework, a multi-objective reward function including response speed reward item, regulation accuracy reward item, operation cost reward item, system stability reward item and environmental benefit reward item is designed to generate a comprehensive performance evaluation index; Based on the comprehensive performance evaluation index, an improved twin-delayed deep deterministic policy gradient algorithm is used to construct actor network and critic network to obtain the reinforcement learning aggregation strategy optimization model.

[0016] Preferably, the reinforcement learning aggregation strategy optimization model is used to receive state data in real time and output an optimal aggregation strategy, and the optimal aggregation strategy is used to send regulation instructions to each flexible resource through a distributed control system, and actual response data is collected for online learning and updating to realize adaptive dynamic optimization, comprising: The reinforcement learning aggregation strategy optimization model is used to acquire real-time state data, and a constraint checking module is used to verify the feasibility of the strategy, if the strategy violates the constraint, a projection algorithm is used to project the action vector into the feasible region to generate the optimal aggregation strategy; Based on the optimal aggregation strategy, regulation instructions are sent to each flexible resource through a communication protocol and actual execution effect data is collected to obtain actual response data; The actual response data is used to calculate an actual reward value and store it in an experience replay buffer, and online learning and updating are started to realize adaptive dynamic optimization of the system.

[0017] The application also provides an energy complex flexible dynamic aggregation system based on reinforcement learning, comprising: The acquisition module is used to acquire operational data of flexible resources within the energy complex. It collects real-time operational parameters of electrical, gas, thermal, and multi-energy coupling equipment via the industrial Ethernet protocol and constructs a continuous-time dynamic system model using the neural network differential equation process modeling technique. The identification module is used to perform preliminary identification of energy aggregation problems based on the continuous-time dynamic system model using a static spiking neural network. After identifying potential aggregation needs, it activates an adaptive dynamic spiking neural network classifier to accurately classify specific energy aggregation problems. The analysis module is used to determine the corresponding data processing strategy and time window parameters by utilizing the different aggregation demand types in the accurate classification results of the specific energy aggregation problem, and to perform alignment analysis on energy data measured at non-equidistant time points using the multilateral stochastic flow matching method. Through measurement value spline enhancement technology and fraction matching, the response characteristics of each flexibility resource are quantitatively evaluated. The response characteristics of each flexibility resource include the randomness of state transitions and the sequentiality of decision-making in the energy aggregation process. The module is used to model the energy aggregation process as a Markov decision process, and to train and optimize it using the twin-delay deep deterministic policy gradient algorithm to generate a reinforcement learning aggregation policy optimization model. The output module is used to receive state data in real time and output the optimal aggregation strategy using the reinforcement learning aggregation strategy optimization model. Using the optimal aggregation strategy, it sends adjustment instructions to each flexible resource through the distributed control system and collects actual response data for online learning and updating to achieve adaptive dynamic optimization.

[0018] The technical solution provided by this invention has the following beneficial effects: 1. This invention uses the process modeling technique of the divine constant differential equation to construct a continuous-time dynamic system model, which can accurately describe the dynamic characteristics of complex multi-element heterogeneous energy systems and effectively solve the problem of insufficient accuracy and adaptability of traditional modeling methods when facing nonlinear, time-varying and uncertain factors.

[0019] 2. This invention employs a two-layer spiking neural network for the identification and classification of aggregation problems. The static spiking neural network is responsible for the initial identification, while the adaptive dynamic spiking neural network is responsible for the precise classification. This structure features high computational efficiency and strong adaptability, and can accurately identify various energy aggregation demands.

[0020] 3. This invention employs a multilateral random flow matching method to perform alignment analysis on high-dimensional energy data measured at non-equidistant time points. Through measurement value spline enhancement technology and fractional matching algorithm, it effectively solves the problems of high-dimensional data processing and non-equidistant time series analysis, and improves the accuracy of resource response characteristic assessment.

[0021] 4. The present application models the energy aggregation process as a Markov decision process, and uses an improved twin-delayed deep deterministic policy gradient algorithm for training optimization, which effectively improves the stability and convergence speed of policy learning through techniques such as double-path network design, priority sampling and adaptive soft update.

[0022] 5. The present application enables the reinforcement learning model to continuously adapt to system changes through an online learning update mechanism, achieving adaptive dynamic optimization and greatly improving the efficiency and accuracy of energy complex flexible resource aggregation. BRIEF DESCRIPTION OF DRAWINGS

[0023] In order to more clearly illustrate the technical solutions of the embodiments of the present disclosure, the drawings needed to be used in the embodiments will be briefly introduced below. The drawings herein are incorporated into the specification and form a part of the specification, which show the embodiments consistent with the present disclosure, and are used to illustrate the technical solutions of the present disclosure together with the specification. It should be understood that the following drawings only show some embodiments of the present disclosure, and therefore should not be regarded as a limitation on the scope, and other related drawings can also be obtained by those skilled in the art without creative labor.

[0024] Figure 1 The flowchart of the energy complex flexible dynamic aggregation method based on reinforcement learning of the present application; Figure 2 The structure diagram of the energy complex flexible dynamic aggregation system based on reinforcement learning of the present application. DETAILED DESCRIPTION

[0025] In order to make the purpose, technical solutions and advantages of the embodiments of the present disclosure clearer, the technical solutions in the embodiments of the present disclosure will be described clearly and completely below in conjunction with the drawings in the embodiments of the present disclosure. Obviously, the described embodiments are only some of the embodiments of the present disclosure, but not all the embodiments. The components of the embodiments of the present disclosure described and shown in the drawings herein can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present disclosure provided in the drawings is not intended to limit the scope of the claimed present disclosure, but only represents selected embodiments of the present disclosure. Based on the embodiments of the present disclosure, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present disclosure.

[0026] It should be noted that: similar reference numerals and letters represent similar items in the following drawings, and therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.

[0027] The term "and / or", merely describes an associated relationship, which means that there can be three relationships, for example, A and / or B, which can represent three cases: A exists alone, A and B exist together, and B exists alone. In addition, the term "at least one" herein means any one of the plurality or any combination of at least two of the plurality, for example, including at least one of A, B and C, which means including any one or more elements selected from the set consisting of A, B and C.

[0028] The application will be further described in detail below in combination with the drawings and examples.

[0029] As shown in the drawings, Figure 1 The application provides a flexible and dynamic aggregation method for an energy complex based on reinforcement learning, comprising the following steps: S1: Obtain the operation data of flexible resources in the energy complex, collect real-time operation parameters of electric, gas, heat and multi-energy coupling devices, and use neural differential equation process modeling technology to build a continuous time dynamic system model; An energy complex refers to a complex system integrating multiple energy form conversion, storage, transmission and utilization functions, which realizes the coordinated and optimized operation of electric energy, thermal energy, gas and other multiple energies through a unified energy management platform. Such a comprehensive system usually contains diversified energy production equipment, storage facilities, conversion devices and terminal energy utilization equipment, forming an organic whole that works together in a coordinated manner. The core feature of the energy complex is its multi-energy coupling characteristic, that is, different forms of energy can be converted into each other through various conversion devices, thereby improving the overall energy utilization efficiency and enhancing the flexibility and reliability of system operation.

[0030] Within the energy complex, energy storage stations serve as important energy buffer facilities, bearing multiple functions such as peak load shifting, frequency regulation and voltage support. These energy storage facilities usually use large-capacity lithium ion battery systems or other advanced energy storage technologies, with the ability of fast response and precise control. Distributed photovoltaic power stations utilize scattered resources such as building rooftops and open spaces to convert solar energy into electricity and connect to the grid nearby, with their output characteristics significantly affected by weather conditions, showing obvious intermittency and volatility. Electric vehicle charging stations are not only energy consumption facilities, but also mobile energy storage resources that can participate in grid regulation after being equipped with bidirectional charging and discharging technology, and their aggregation effect is increasingly important in urban energy systems.

[0031] Flexible resources refer to various devices and systems that can quickly adjust their energy output, consumption, or storage status according to system demand. These resources have adjustable and controllable characteristics and can respond to system dispatch instructions at different time scales. Among electrical flexible resources, battery energy storage systems are the most important adjustment resources due to their millisecond response speed and high-precision power control capabilities. Although supercapacitors have relatively small energy storage capacity, their ultra-fast charging and discharging characteristics make them particularly suitable for handling transient power fluctuations and frequency regulation tasks. Interruptible loads, including industrial production lines and large refrigeration equipment, can temporarily stop running in emergency situations, providing emergency regulation capabilities for the system.

[0032] Gas flexible resources mainly include various gas storage and regulation facilities. Underground gas storage takes advantage of geological structures to provide large-scale, long-period natural gas storage capacity, playing an important role in supply and demand balance and emergency protection. Although the capacity of ground gas storage facilities such as high-pressure spherical tanks is relatively small, their regulation response is more flexible and can quickly respond to short-term gas regulation needs. The regulation characteristics of these gas resources are mainly reflected in pressure regulation, flow control, and storage management.

[0033] Thermal flexible resources participate in system regulation through the storage and release of thermal energy. Thermal storage tanks store thermal energy using the sensible heat characteristics of water or other media, playing a role in peak regulation and stable operation in heating systems. Thermal storage materials, especially phase change materials, can absorb or release a large amount of latent heat during the phase change process, achieving high-density thermal energy storage. The response characteristics of these thermal resources are influenced by heat conduction and convective heat transfer processes, and usually have certain thermal inertia, with relatively long response times but strong persistence.

[0034] Multi-energy coupling devices are key hubs that connect different energy forms and achieve mutual conversion between energy forms. Combined heat and power (CHP) units simultaneously generate electric and thermal energy, and by adjusting their operation mode, they can be optimally configured between the electricity market and the heat market. Electric boilers convert electric energy into thermal energy, increasing power consumption during periods of renewable energy abundance while providing a heat source for heating systems. The operating parameters of these coupling devices include conversion efficiency, response time, operating constraints, and other dimensions.

[0035] The operational data for various flexible resources encompasses a wealth of technical parameters and status information. Real-time power reflects the instantaneous operating status of equipment and serves as the direct basis for system scheduling. Capacity parameters, including rated capacity, available capacity, and remaining capacity, determine the resource's adjustment potential. Response time describes the time required for equipment to reach its target state from receiving a command, directly affecting its applicability in different adjustment services. Charging and discharging efficiency reflects the loss characteristics during energy conversion, impacting the overall economic efficiency of the system. Operational constraints include technical boundary conditions such as maximum and minimum output limits, ramp rate limits, and start-stop time constraints. These constraints ensure that equipment operates within a safe and reliable range and are also crucial factors that must be considered when formulating scheduling strategies.

[0036] In this step, the flexible dynamic aggregation system of the energy complex (hereinafter referred to as the system) collects operational data of various flexible resources within the energy complex in real time through a high-speed data channel established by the industrial Ethernet protocol. For electrical resources, the system collects electrical parameters including voltage, current, active power, reactive power, and power factor, as well as characteristic parameters of battery energy storage devices such as state of charge, charge-discharge cycle count, and health status. The energy storage devices employ a high-frequency sampling rate of 10Hz to capture rapid dynamic response characteristics. For thermal resources, the system collects thermal parameters such as temperature, pressure, flow rate, and thermal power, as well as characteristic parameters of thermal energy storage devices such as heat storage capacity, heat loss rate, and heat transfer coefficient, using a standard sampling rate of 1Hz. For gaseous resources, the system collects gas parameters such as pressure, flow rate, composition, and calorific value, as well as characteristic parameters of gas energy storage devices such as gas storage capacity, compression ratio, and gas release rate. For multi-energy coupling devices, such as combined heat and power units, electric refrigeration / electric heating equipment, and hydrogen power generation equipment, the system simultaneously collects input and output parameters and conversion efficiencies of multiple energy forms.

[0037] The collected real-time operating parameters (raw data) undergo preprocessing, starting with outlier detection and removal. The system employs the three sigma (3σ) rule to calculate the mean μ and standard deviation σ for each data category, marking data points falling outside the range [μ-3σ, μ+3σ] as outliers and removing them. Next, data standardization is performed using the Z-score standardization method, converting data of different dimensions and scales into a standard distribution with a mean of 0 and a standard deviation of 1, eliminating the impact of dimensional differences on subsequent modeling. Finally, time series alignment is performed to address inconsistencies in sampling frequencies across different devices and parameters, establishing a unified time baseline.

[0038] Based on the preprocessed high-quality data, the system uses the neural network ordinary differential equation process modeling technique to construct a continuous-time dynamic system model. This technique combines traditional ordinary differential equations with deep neural networks, enabling it to accurately describe the dynamic characteristics of complex nonlinear systems.

[0039] Firstly, the dynamic process of the energy device is represented as a differential equation form: where x(t) represents the system state vector, u(t) represents the control input vector, and θ represents the parameters to be learned. Then a deep neural network with a residual network architecture is used to approximate the function f, and the gradient is calculated by the adjoint sensitivity method to achieve efficient optimization of the parameters.

[0040] Representing the dynamic process of the energy device as a differential equation is a basic step in establishing an accurate mathematical model. This representation method is derived from the continuous-time characteristics of physical systems and can accurately describe the evolution of system state over time. In this differential equation framework, the system state vector contains all key variables that describe the current operating state of the energy device, such as the state of charge, temperature, internal resistance, and other internal state parameters of the energy storage device. The control input vector represents the external control signals applied to the system, such as the charge and discharge power set value, valve opening adjustment, and other controllable input quantities. The set of parameters to be learned includes all parameters in the system model that need to be determined through data learning, which may include device physical property parameters, conversion efficiency coefficients, time constants, and other parameters.

[0041] Traditional physical modeling methods usually construct differential equations based on first principles and empirical formulas, but this method often fails to accurately capture all nonlinear characteristics and coupling relationships when faced with complex energy systems. In particular, when considering factors such as device aging, environmental changes, and operating condition diversity, the limitations of pure physical models become more apparent. Therefore, the invention innovatively uses a deep neural network to approximate the state evolution function in the differential equation, which combines the interpretability of physical modeling with the powerful fitting ability of machine learning.

[0042] The introduction of the residual network architecture is a key innovation point of this method. The residual network uses a jump connection mechanism, allowing the network to learn the state change quantity rather than the absolute state value, which is highly consistent with the idea of describing the state derivative in the differential equation. Specifically, each basic block of the residual network contains a shortcut connection with an identity mapping, making the network output the result of the input plus a residual function. This design not only alleviates the gradient vanishing problem in deep network training, but more importantly, it matches the essence of the differential equation in a physical sense, because the differential equation itself is describing the incremental change of the system state.

[0043] The specific architecture of the neural network fully considers the characteristics of the energy system. The input layer of the network receives both the current state vector and the control input vector, forming a complete system input. The hidden layer uses a multi-layer fully connected structure, with each layer equipped with appropriate activation functions and regularization mechanisms. The choice of activation function is particularly important, as it needs to express nonlinear characteristics while maintaining good gradient propagation characteristics. The output layer of the network produces a state derivative vector, which has the same dimension as the state vector and directly corresponds to the right-hand function of the differential equation.

[0044] The adjoint sensitivity method is the core technology for training such neural ordinary differential equation models. This method is derived from the adjoint equation theory in optimal control theory. By constructing an adjoint equation that is dual to the original differential equation, the gradient of the objective function with respect to the model parameters can be efficiently calculated. Specifically, the adjoint sensitivity method first solves the original differential equation forward to obtain the state trajectory, and then solves the adjoint equation backward to calculate the gradient information. The biggest advantage of this method is that its computational complexity is independent of the time step number. No matter how long the integration time is, the complexity of gradient calculation remains constant, which is extremely important for processing long-time sequence data of energy systems.

[0045] The parameter optimization process is carried out iteratively, with each iteration including two stages: forward propagation and backward propagation. In the forward propagation stage, the system uses the current parameters to solve the differential equation and obtains the predicted state trajectory. In the backward propagation stage, the adjoint sensitivity method is used to calculate the gradient of the loss function with respect to all model parameters, and then a gradient descent algorithm is used to update the parameters. This process continues until the model performance reaches a satisfactory level or the convergence condition is met. The entire optimization process maintains the continuity and interpretability of the physical model, while fully utilizing the advantages of deep learning in handling complex nonlinear relationships, providing a new technical approach for energy system modeling.

[0046] In addition, in S1, for the battery energy storage system, the first model considers the state of charge variation, the variation of internal resistance with temperature and aging degree, the nonlinear characteristics of power conversion efficiency, and the capacity attenuation effect. For the gas turbine, the second model includes start-stop characteristics, operating inertia, output adjustment rate limit, and the variation of thermal-electric conversion efficiency with load. For the thermal energy storage system, the third model integrates the dynamic characteristics of heat transfer process, heat loss mechanism, and nonlinear characteristics of phase change materials. These sub-models (first, second, and third models) are trained and verified by a variable step ordinary differential equation solver, and finally integrated into a continuous-time system model that can accurately describe the dynamic characteristics of various resources, providing a solid foundation for subsequent aggregation problem identification and strategy optimization.

[0047] S2: Based on the continuous-time dynamic system model, a static spiking neural network is used to initially identify the energy aggregation problem. After identifying the potential aggregation demand, an adaptive dynamic spiking neural network classifier is activated to accurately classify the specific energy aggregation problem. Step S2 implements the energy aggregation problem identification and classification using a two-layer spiking neural network. This step innovatively employs biomimetic neural computing principles, simulating the information encoding and processing mechanisms of biological neurons to construct an efficient and robust framework for identifying energy aggregation problems. First, in S2, the state data generated from the continuous-time dynamic system model in step S1 is normalized to ensure that all input feature dimensions are unified within the [0,1] interval, thereby improving the network's learning efficiency and stability.

[0048] A static spiking neural network serves as the first layer of the recognizer, responsible for rapidly detecting potential aggregation needs in the system. This network employs a biologically sound Leaky Integrate-and-Fire (LIF) neuron model, utilizing the membrane potential dynamic equation... This simulates the working mechanism of a biological neuron. When a neuron receives an input stimulus, the membrane potential v gradually accumulates; when it reaches the threshold voltage Vth... th At this time, the neuron generates a pulse and resets the membrane potential to the resting potential V. rest This mechanism enables the network to efficiently process time-varying signals and extract temporal information features. The input layer contains 64 LIF neurons that receive system state data; the hidden layer contains 128 fully connected LIF neurons, with connection weights initialized using the Xavier method; the output layer has 4 neurons, corresponding to four basic aggregation needs: load peak-valley regulation, frequency regulation, voltage support, and emergency response. The network converts continuous state values ​​into pulse sequences through a time-coding mechanism and uses a winner-take-all mechanism to determine the final recognition result.

[0049] Once the static spiking neural network identifies a potential aggregation requirement, an adaptive dynamic spiking neural network classifier (the second-layer classifier) ​​is activated for accurate problem classification. This classifier employs a Growing When Required (GWR) mechanism, which dynamically adjusts the network structure based on the novelty of the input data. Specifically, it calculates the Euclidean distance between the input feature vector and the existing neuron prototype vectors, and when the minimum distance exceeds a preset novelty threshold θ... sim (Usually set to 0.8) and network activity a t Below the activity threshold θ actWhen the value is set to 0.1, the current input is determined to be a new pattern, and new neurons are automatically added to the network to expand its capacity. The initial weights of the new neurons are determined by linear interpolation, and their prototype vector is the average of the weights of the best-matching neuron and the second-best-matching neuron.

[0050] Network connection weights are updated using an Adaptive Spike-Timing-Dependent Plasticity (Ad-STDP) rule. Inspired by biological synaptic plasticity mechanisms, this rule adjusts synaptic strength based on the time difference between preceding and following neuronal pulses. The current neuronal pulse occurs before the following neuronal pulse. When the value is greater than 0, the connection weight is increased; conversely, it is decreased. The update formula is as follows: (Activity), among which Let A be the learning rate. The function f is time-dependent, and f(activity) is an adaptive function that dynamically adjusts the learning intensity based on the historical activity level of neurons, effectively preventing overlearning and catastrophic forgetting. The classifier output layer contains 12 neurons, corresponding to more specific aggregate sub-problem types, including short-term peak shaving, medium-term peak shaving, long-term energy storage, primary frequency adjustment, and secondary frequency adjustment, providing precise guidance for subsequent strategy formulation.

[0051] This dual-layer spiking neural network architecture fully combines the rapid response capability of static networks with the adaptive learning capability of dynamic networks. It can efficiently identify complex aggregated demands in energy systems and continuously improve classification performance as system operation experience is accumulated, forming a deep understanding of the operation mode of energy complexes.

[0052] S3: Obtain the different aggregation demand types in the accurate classification results of the specific energy aggregation problem, use the multilateral stochastic flow matching method to perform alignment analysis on the energy data measured at non-equidistant time points, and use the measurement value spline enhancement technology and fraction matching to quantitatively evaluate the response characteristics of each flexibility resource. The response characteristics of each flexibility resource include the randomness of state transitions and the sequentiality of decision-making in the energy aggregation process. In this embodiment, step S3 achieves resource response characteristic assessment based on multilateral random flow matching. This step aims to solve the problem of high-dimensional data processing with non-equidistant time points commonly found in energy systems, and to accurately quantify the dynamic response characteristics of various flexible resources through advanced statistical learning methods.

[0053] Specifically, according to the specific aggregation problem type identified in step S2, a corresponding data processing strategy and time window parameter are determined. For fast response type requirements (such as frequency adjustment), a short time window (seconds to minutes) and a high sampling rate are used; for medium and long-term adjustment requirements (such as load peak valley adjustment), a longer time window (hours) and a moderate sampling rate are used.

[0054] Then, a plurality of marginal distributions are constructed, each distribution corresponding to a data feature of a specific resource type or time period. For electrical energy storage resources, the marginal distribution includes key dimensions such as power output, response time, charging and discharging efficiency, and cycle life decay coefficient. For thermal energy storage resources, the marginal distribution covers parameters such as temperature change gradient, heat storage capacity, heat release rate, and heat loss coefficient. For gas energy storage resources, the marginal distribution includes indicators such as pressure change rate, volume flow, and compression efficiency. The system uses kernel density estimation method for non-parametric modeling, without making prior assumptions about the form of data distribution, and uses Gaussian kernel function The original data points are smoothed. The bandwidth parameter h is determined by k-fold cross-validation optimization, which avoids overfitting while ensuring fitting accuracy.

[0055] To address the challenge of non-equidistant time series data processing, the measurement value spline technique is used to significantly enhance the processing capability of irregular snapshot times. In actual energy systems, due to communication delays, equipment failures, or network congestion, data collection time points are often irregularly distributed, resulting in poor results with traditional interpolation methods. First, determine an appropriate spline node sequence, set the node density according to the time interval distribution of the observed data, and increase the number of nodes in data-intensive areas. Then construct a cubic B-spline basis function , which is a cubic polynomial in four consecutive node intervals and zero in other intervals, with C2 continuity to ensure the smoothness of the interpolation results. The observed values at irregular time points are expressed as a linear combination of the spline function: The coefficients c i are determined by minimizing the reconstruction error, while introducing a smoothing term control function to control the smoothness and prevent overfitting of noise.

[0056] In high-dimensional data space, the fractional matching algorithm is used to avoid the problems of dimension disaster and numerical instability. This algorithm minimizes the difference between the real data distribution and the fractional function (gradient of the logarithmic probability density function) of the distribution model, effectively avoiding the difficulty of calculating the normalization constant in traditional maximum likelihood estimation. The system first calculates the logarithmic density gradient of the real data distribution as the real fractional function, and the logarithmic density gradient of the distribution model as the fractional function of the distribution model. Then construct the fractional matching objective function:

[0057] The least square method is used to calculate the square error between the two score functions and integrate the sum. Finally, the random gradient descent optimization parameter θ is realized by the Adam optimizer, and the accurate matching of each marginal distribution is realized through multiple iterations.

[0058] The system finally outputs the quantified resource response characteristics, including response time constant (reflecting resource adjustment speed), adjustment accuracy (reflecting resource control accuracy) and stability index (reflecting resource output stability). These characteristic indexes form a weight matrix , where T is the time step, N is the number of resource types, M is the number of characteristic dimensions, and is the subsequent reinforcement learning strategy optimization. The accurate resource capability description makes the aggregated decision fully consider the actual characteristics and coordination potential of each resource.

[0059] S4: Model the energy aggregation process as a Markov decision process and train and optimize it using the twin-delayed deep deterministic policy gradient algorithm to generate a reinforcement learning aggregation strategy optimization model; Step S4 realizes the construction of the reinforcement learning aggregation strategy optimization model. This step formalizes the dynamic aggregation problem of the energy complex into a Markov decision process (MDP) framework and learns the optimal decision strategy through advanced reinforcement learning algorithms. The system first defines the five-tuple of MDP, where S is the state space, A is the action space, P is the state transition probability function, R is the reward function, and γ is the discount factor (usually set to 0.95).

[0060] The construction of the state space S fully considers the multi-dimensional characteristics and time-varying nature of the energy system, and includes five main dimensions: load demand dimension, renewable energy dimension, energy storage state dimension, power quality dimension and economic signal dimension. The load demand dimension includes the total load demand (MW) at the current time, the load prediction uncertainty (standard deviation) and the load change rate (MW / min). The renewable energy dimension includes wind power output prediction (MW), photovoltaic output prediction (MW), prediction confidence interval and weather state coding. The energy storage state dimension includes the state of charge of electrical energy storage (%), the temperature distribution of thermal energy storage (℃), the pressure state of gas energy storage (bar) and the available capacity of each energy storage device. The power quality dimension includes system frequency (Hz), key node voltage amplitude (p.u.), power factor and harmonic distortion rate. The economic signal dimension includes real-time electricity price (yuan / MWh), auxiliary service price and carbon emission factor. The total dimension of the state vector is 17, which fully captures the system operating state information.

[0061] The action space A adopts a hybrid continuous-discrete design. The continuous action part is defined as the participation degree allocation vector of each type of resource , where represents the contribution proportion of the i-th type of resource in the aggregation, satisfying the normalization constraint . The discrete action part denotes the aggregation mode selection: 0 - economic optimization mode, 1 - fast response mode, 2 - stable operation mode, 3 - emergency support mode. The complete action vector is , which can flexibly express complex resource scheduling strategies.

[0062] The reward function adopts a weighted multi-objective form to comprehensively evaluate the multi-aspect performance of the aggregation strategy. The response speed reward item evaluates the timeliness of the system's response to the aggregated demand by comparing the deviation of the actual response time from the target response time and imposing an additional penalty on the overtime response. The regulation accuracy reward item evaluates the matching degree of the aggregation output and the target value, considering both the overall accuracy and the individual deviation of each resource. The operation cost reward item calculates the economic cost of strategy execution, including energy cost, equipment wear cost, and start-stop cost. The system stability reward item evaluates the impact of the strategy on system stability, including frequency stability, voltage stability, and power oscillation suppression effect. The environmental benefit reward item encourages renewable energy consumption and punishes fossil fuel use, promoting low-carbon operation. The total reward function combines these five indicators by weighting: w1=0.25, w2=0.30, w3=0.20, w4=0.15, w5=0.10, reflecting the relative importance of different performance targets.

[0063] The system adopts the Enhanced Twin Delayed Deep Deterministic Policy Gradient (ETD3) algorithm for training and optimization. This algorithm has undergone multiple optimization designs tailored to the characteristics of the energy aggregation problem. The Actor network adopts a multi-layer structure, with the input layer receiving a 17-dimensional state vector and performing normalization through Layer norm alization; the hidden layer uses the Swish activation function and introduces residual connections to alleviate the gradient vanishing problem; the output layer is divided into continuous action part (uses Softmax activation to ensure constraints) and discrete action part (uses Gumbel-Softmax to realize differentiable sampling). The Critic network adopts a dual-path design, with the state path and action path processed separately and merged in the fusion layer. The two Critic networks have the same structure but independent parameter initialization, improving the robustness of Q value estimation.

[0064] During training, the system adopts a priority experience replay mechanism, assigns sampling probability according to TD error, and improves the utilization rate of effective samples. The target network update adopts an adaptive soft update mechanism to achieve the balance between rapid learning in the early stage of training and stable convergence in the later stage. The strategy update introduces a delay mechanism and noise regularization, updates the Actor network every 3 steps, adds truncated Gaussian noise to the target action, and improves the exploration ability and robustness of the strategy. The learning rate is dynamically adjusted using the cosine annealing strategy, which changes periodically during training, helping the model to jump out of the local optimum. When the average reward of 100 consecutive episodes changes by less than 0.1% and the policy loss stabilizes at 10 -4 After the above, it is considered that the training is converged, and a high-performance reinforcement learning aggregated strategy optimization model is finally obtained.

[0065] S5: Real-time receiving state data and outputting the optimal aggregated strategy by using the reinforcement learning aggregated strategy optimization model, and using the optimal aggregated strategy to send adjustment instructions to each flexibility resource through a distributed control system, and collecting actual response data for online learning and updating to realize adaptive dynamic optimization.

[0066] This step deploys the trained reinforcement learning model to the actual running environment to realize real-time optimization and scheduling of the flexibility resources of the energy complex. At the beginning, the reinforcement learning aggregated strategy optimization model is deployed to the edge computing node of the energy management system, and a stable and low-delay communication connection is established with each resource controller through high-speed Ethernet. The edge computing architecture reduces data transmission delay and improves system response speed, and is particularly suitable for fast response scenarios that require millisecond-level decision-making.

[0067] During runtime, real-time state data is collected and processed at a period of 1 second, including device operating parameters obtained through the SCADA system, real-time weather data obtained from the weather station, short-term load prediction results obtained from the load prediction system, and real-time price signals obtained from the energy market. After preprocessing and feature extraction, these multi-source heterogeneous data form a 17-dimensional state vector input into the Actor network. The network quickly calculates and outputs a continuous action vector representing the optimal participation ratio of each type of resource, as well as a discrete aggregated mode selection result.

[0068] Before action execution, the system verifies the feasibility of the strategy through a constraint checking module. Constraint checking includes two levels: device physical constraints and system operation constraints. Device physical constraints include power limits ( ), capacity limits ( ), ramp rate limits ( ), and minimum operation / shutdown time of devices, etc. System operation constraints include power balance constraints ( ), line flow constraints ( ), and node voltage constraints ) and system reserve capacity constraints, etc. If the policy violates any constraints, the system projects the action vector into the feasible region using a quadratic programming projection algorithm, with the projection formula , ensuring that the finally executed policy satisfies both physical constraints and is as close as possible to the original optimized policy.

[0069] The adjusted policy is sent to each resource controller through a standard communication protocol (e.g., Modbus, IEC 61850, or custom API). The control instruction structure contains three key parts: target output value (e.g., power setpoint), response time requirement, and priority flag. After receiving the instruction, the resource controller confirms the feasibility of execution based on its own state and capability boundaries, and returns confirmation information. If a resource cannot execute the instruction, the policy is quickly adjusted and the task is redistributed to other available resources. During execution, the system continuously monitors the actual output and response characteristics of each resource, recording key performance indicators such as response delay, adjustment deviation, energy conversion efficiency, and operating cost.

[0070] To ensure that the model can continuously adapt to system changes, an online learning update mechanism is used. At the end of each control cycle, the system calculates the true reward value based on the actual execution effect, and stores the state transition sample (s t , a t , r t , s t+1 ) into the experience replay buffer. When the buffer accumulates enough new samples (usually 1000), the system starts the incremental learning process. Online learning uses a smaller learning rate (10 -5 ) and batch size (32) to avoid excessive disturbance to the trained effective knowledge. To prevent the problem of catastrophic forgetting, the system uses the Elastic Weight Consolidation (EWC) technology, adding a regularization term to the loss function, where is the importance weight, measuring the degree of influence of parameter on model performance, is the key parameter value. The importance weight is estimated by the Fisher information matrix, ensuring the stability of important parameters.

[0071] The system updates model parameters every hour, adapting to new operating modes and changes in device characteristics while maintaining historical valid knowledge. As the system runs for a longer period of time, the model performance continues to improve, significantly improving resource coordination efficiency and accuracy. Practical applications show that compared to traditional static optimization methods, this adaptive dynamic optimization system can improve renewable energy consumption by 20%-35%, reduce operating costs by 15%-25%, and shorten adjustment response time by 40%-60%, providing strong technical support for efficient, economic, and environmentally friendly operation of energy complexes.

[0072] An industrial park energy complex managed by a certain municipal power supply company is a typical multi-energy coupling system. The park covers an area of about 15 square kilometers and includes 120 manufacturing enterprises, three large commercial complexes, and 20,000 residential users. The park has 50 MW of distributed photovoltaic power stations, 20 MW / 40 MWh of lithium battery energy storage stations, 15 MW of gas turbine combined heat and power units, thermal energy storage systems, and electric vehicle charging station groups, among other multi-type energy facilities. As renewable energy penetration continues to increase and energy demand becomes increasingly diverse, the park faces challenges such as large load peak-valley differences, difficulty in renewable energy consumption, and fluctuations in power quality, necessitating the establishment of an intelligent energy aggregation management system.

[0073] The power supply company launched an energy complex dynamic aggregation project based on reinforcement learning in 2024, aiming to achieve coordinated and optimized operation of various flexible resources in the park through advanced artificial intelligence technology. The project team deployed a comprehensive data collection network within the park, established a centralized energy management platform, and developed an intelligent dispatching system with self-learning capabilities. After six months of system construction and three months of trial operation, the project has achieved significant results in improving renewable energy consumption, reducing operating costs, and improving power quality.

[0074] The first phase of the project focuses on establishing a comprehensive data collection system and dynamic system model. In terms of data collection, the power supply company installed over 500 intelligent sensors and monitoring devices within the park, building a high-speed and reliable communication network through industrial Ethernet. For the 20 MW lithium battery system of the energy storage station, a high-frequency sampling rate of 10 Hz is used to monitor key parameters such as voltage, current, temperature, and state of charge of each battery cluster in real time, as well as health status indicators such as charge and discharge cycle count, internal resistance change, and capacity decay. The monitoring system of the distributed photovoltaic power station collects output power, voltage, current, and environmental parameters such as module temperature and irradiance from each inverter at a frequency of 1 Hz.

[0075] The data acquisition of the gas turbine cogeneration unit covers multiple parameters such as the operation state, fuel consumption, power output, thermal output, and exhaust gas temperature. The thermal energy storage system mainly monitors the temperature distribution of the storage tank, hot water flow, and changes in the amount of stored heat. The electric vehicle charging station group collects real-time data such as the usage state of each charging pile, charging power, and vehicle information. These raw data are preprocessed by the preprocessing module to detect outliers and standardize the data. The three-sigma rule is used to remove about 2% of the abnormal data points, and the Z-score standardization is used to eliminate the dimensional differences between different parameters.

[0076] During the dynamic modeling phase, the project team used neural differential equation techniques to establish accurate dynamic models for each type of major equipment. Taking the lithium battery energy storage system as an example, the model unified the complex dynamic processes such as state of charge change, power response characteristics, temperature influence, and aging effects into continuous-time differential equation form. By collecting three months of operation data from the energy storage station, including response characteristics under different charging and discharging modes, performance under various environmental temperature conditions, and voltage response characteristics within different state of charge ranges, the trained model achieved an root mean square error of less than 1.5% in state of charge prediction, a voltage prediction error of less than 0.03V, and a power response time prediction error of less than 30ms.

[0077] For the gas turbine unit, the project team paid special attention to the dynamic characteristic modeling of its start-stop process and load regulation process. By recording the operation data of the unit under different seasons and load levels, including cold start, hot start, and detailed parameter changes during different load ramping processes, the dynamic model can accurately predict the start-up time, load response speed, and fuel consumption rate of the unit. Model verification results show that the power prediction accuracy reaches 2.1% of the rated power, and the thermal efficiency prediction error is controlled within 0.8 percentage points.

[0078] Based on dynamic modeling, the project team developed an intelligent identification system for aggregation problems based on a double-layer pulse neural network. The first layer of the static pulse neural network is responsible for quickly identifying potential aggregation needs in the system. In actual operation, when the park experiences a sharp increase in load during the summer peak electricity consumption period, the system can identify the load peak-valley regulation demand within 50 milliseconds. When a nearby 110kV substation fails, causing the system frequency to deviate from the rated value by 0.3Hz, the network immediately identifies the frequency regulation demand. When the voltage of an important load center in the park drops to 0.92 per unit, the system quickly identifies the voltage support demand.

[0079] The recognition accuracy of the static pulse neural network reached 91.5% in actual operation, basically meeting the demand of rapid preliminary screening. For the identified potential aggregation demand, the system automatically activates the second layer of adaptive dynamic pulse neural network for accurate classification. The dynamic classifier uses the grow-when-needed mechanism, which can automatically expand the network structure according to new emerging operation modes. In the first two months of the project operation, the network gradually expanded from the initial 36 neurons to 68 neurons, and the newly added neurons were mainly used to represent new aggregation scenarios such as concentrated charging of electric vehicles in the park and coordinated shutdown of industrial loads.

[0080] The dynamic classifier further subdivides the basic four types of aggregation demand into 12 specific subcategories. For example, load peak-valley regulation is subdivided into short-term peak shaving, medium-term peak shaving, and long-term energy storage; frequency regulation is subdivided into primary frequency regulation and secondary frequency regulation; voltage support is subdivided into voltage and reactive power support, and harmonic control, etc. This refined classification provides more accurate guidance for subsequent resource allocation and strategy formulation. In actual application, when the system identifies short-term peak shaving demand, it will preferentially call fast-responding energy storage resources; for long-term energy storage demand, it will coordinate thermal energy storage systems and interruptible loads, etc.

[0081] The third phase of the project focuses on the time alignment of multi-source heterogeneous data and the quantitative evaluation of resource characteristics. Due to the inconsistent data sampling frequencies of different devices in the park, the power monitoring system uses 50Hz high-frequency sampling, the thermal system uses 1Hz medium-frequency sampling, and the gas system only uses 0.2Hz low-frequency sampling. This inconsistency in time brings challenges to unified analysis. The project team uses the multi-marginal stochastic flow matching method to construct marginal distribution models for electric energy storage, thermal energy storage, and gas energy storage resources.

[0082] Taking the 20MW lithium battery energy storage system in the park as an example, the project team constructed the marginal distribution of key parameters such as power output, response time, charge and discharge efficiency, and cycle life attenuation coefficient. Through three months of continuous monitoring data analysis, it was found that the power output distribution of the energy storage system showed obvious bimodal characteristics, corresponding to the charging and discharging two main operation modes. The response time distribution was mainly concentrated in the range of 200-800 milliseconds, with an average response time of 450 milliseconds. The charge and discharge efficiency showed nonlinear characteristics at different power levels, with the highest efficiency of 94.5% at 50% rated power and an efficiency of about 91.2% at rated power.

[0083] For the thermal energy storage system, the project team focused on analyzing the distribution characteristics of the temperature gradient, heat storage capacity, and heat release rate of the thermal storage tank. The monitoring data showed that the temperature gradient of the thermal storage tank was typically in the range of 0.8-3.2°C / min, and the heat release rate reached a peak near the phase change material storage area. The response characteristic analysis of the gas turbine system showed that the cold start time was about 25 minutes, the hot start time was about 8 minutes, and the load ramping rate was 5% of the rated power per minute.

[0084] Through the measurement spline technique, the project team successfully converted these non-equidistant observation data into continuous time series representation. Using the fractional matching algorithm, the response characteristics of various resources were accurately quantified, and a complete set of characteristic parameters including response time constant, regulation accuracy, and stability index were obtained. These parameters provide accurate resource capability description for subsequent reinforcement learning model training.

[0085] Based on the aforementioned analysis results, the project team modeled the park energy aggregation problem as a Markov decision process, constructing a 17-dimensional state space and a mixed continuous-discrete action space. The state space includes the current load demand of the park, renewable energy output prediction, state of each energy storage device, power quality indicators, and real-time electricity price. In actual operation, these state quantities are updated every second, providing real-time system state perception for intelligent decision-making.

[0086] The action space is designed as a participation allocation vector for various flexible resources plus an aggregation mode selection. The continuous action part includes the output allocation proportion of resources such as energy storage systems, combined heat and power units, and interruptible loads, while the discrete action part includes economic optimization, fast response, stable operation, and emergency support four aggregation modes. The project team designed a multi-objective reward function that considers response speed, regulation accuracy, operation cost, system stability, and environmental benefits, and determined the weight coefficients of each objective through the analytic hierarchy process.

[0087] The reinforcement learning model uses an improved twin-delayed deep deterministic policy gradient algorithm for training. The actor network receives 17-dimensional state input and outputs resource allocation strategy through a three-layer hidden network, while the critic network uses a double-path design to evaluate action value. The model training uses the historical operation data of the park for the first six months, including typical operation scenarios in spring, summer, autumn, and winter, as well as handling cases of various abnormalities and emergencies. The training process uses priority experience replay and adaptive soft update mechanism, and the model converges after 1500 iterations, with an average reward of -42.3 on the validation set and an aggregation accuracy of 96.8%.

[0088] The trained reinforcement learning model was officially deployed in the park energy management system in May 2025, starting to guide the coordinated operation of various flexible resources in real time. The system receives real-time state data every 1 second, sends adjustment instructions to resource controllers after verifying the feasibility of the strategy through the constraint checking module. In actual operation, when the system detects that the photovoltaic output suddenly drops by 8 MW due to cloud cover, the intelligent scheduling system formulates a response strategy within 1.2 seconds: the energy storage system immediately outputs 5 MW of power, the gas turbine increases its output by 2 MW, and 1 MW of interruptible load is reduced, successfully maintaining system power balance.

[0089] In a typical summer peak load regulation scenario, when the park load reaches 52 MW at 19:30 in the evening, exceeding the power supply capacity by 3 MW, the system automatically activates the multi-resource coordinated regulation strategy. The energy storage system outputs 20 MW of rated power for 30 minutes, the combined heat and power unit increases its output from 13 MW to 15 MW, and 2 MW of industrial user peak load shifting is coordinated. The entire regulation process is smooth and orderly, successfully resolving the power supply tension and avoiding the occurrence of power cut.

[0090] The adaptive learning capability of the system is fully demonstrated in actual operation. When a batch of electric vehicle charging piles are put into operation in the park, the system automatically adapts to the new load mode through online learning mechanism. In the first two weeks of operation, the system continuously optimizes model parameters by collecting actual response data of electric vehicle charging, enabling the newly added resources to effectively participate in the coordinated scheduling of the park. Similarly, when the response characteristics of a gas turbine change due to equipment aging, the system can also automatically adjust the corresponding control strategy through continuous learning.

[0091] After three months of trial operation, the project has achieved significant results in many aspects. In terms of renewable energy consumption, the utilization rate of photovoltaic power generation in the park has increased from 82.4% to 94.7%, with an average monthly increase of clean energy consumption of about 180 MWh. In terms of operation cost control, through intelligent resource scheduling and peak valley difference arbitrage, the park saves about 150,000 yuan in monthly electricity charges, while reducing the number of invalid starts and stops of equipment, prolonging the service life of equipment.

[0092] In terms of power quality improvement, the voltage qualification rate of the park after system regulation has increased from 96.8% to 99.2%, the standard deviation of frequency deviation has decreased from 0.08 Hz to 0.03 Hz, and the power factor has remained above 0.95. In terms of emergency response capability, the response time of the system to various disturbances has been shortened by an average of 45%, and the regulation accuracy has been improved by more than 30%. In terms of environmental benefits, through optimization of energy allocation and improvement of renewable energy utilization rate, the park reduces carbon dioxide emissions by about 850 tons per month.

[0093] The successful implementation of the project provides more stable and reliable power supply for enterprises in the park, improves the business environment, and also accumulates valuable experience for the construction of smart grid for the power supply company. With the extension of system operation time and the increase of data accumulation, the system performance is expected to be further improved, providing strong technical support for the construction of new power system.

[0094] In step S1, the operation data of flexible resources in the energy complex are obtained, real-time operation parameters of electric, gas, heat and multi-energy coupling devices are collected through industrial Ethernet protocol, and a continuous time dynamic system model is constructed by using neural differential equation process modeling technology, which specifically includes: S1.1: Obtain the operation data of the flexible resources, remove the data points outside the normal range by using the three sigma rule, and eliminate the dimension difference by Z-score standardization to generate the preprocessed operation data; In step S1.1, real-time monitoring of various resources in the energy complex is realized through industrial Ethernet protocol, real-time power, capacity, charge and discharge efficiency and other parameters of electric resources, pressure, flow, reserve and other states of gas resources, temperature, heat capacity, heat transfer coefficient and other characteristics of heat resources, and operation mode and conversion efficiency of multi-energy coupling devices are collected. The data sampling frequency is determined according to the dynamic response characteristics of the resources, the electric energy storage device adopts 10Hz high frequency sampling, and the thermal energy storage device adopts 1Hz standard sampling. The collected data is subjected to outlier detection, and the data points outside the normal range are removed by using the 3σ rule, and the dimension difference is eliminated by Z-score standardization.

[0095] Specifically, step S1.1 mainly completes the acquisition and preprocessing of the operation data of various flexible resources in the energy complex. A high-performance, low-delay data acquisition network is established through industrial Ethernet protocol to realize comprehensive monitoring of electric, gas, heat and multi-energy coupling devices. Different sampling strategies are adopted for different types of energy equipment during data acquisition. For electric resources with fast response speed and significant dynamic characteristics, such as battery energy storage systems, super capacitors and power electronic devices, the system adopts a high-frequency sampling rate of 10Hz to ensure that millisecond-level transient response characteristics can be captured; for heat resources with relatively slow response, such as heat storage devices and phase change material heat storage systems, a standard sampling rate of 1Hz is adopted; for gas resources and multi-energy coupling devices, the sampling rate is set to 0.5-5Hz according to their dynamic characteristics. This differentiated sampling strategy not only ensures the accuracy of key data acquisition, but also avoids the storage and processing pressure caused by data redundancy.

[0096] The collected raw data types are diverse, including voltage, current, power, frequency, temperature, pressure, flow, energy storage state, conversion efficiency, and hundreds of other parameters. These raw data inevitably have problems such as outliers, noise, and missing values, which need to be systematically preprocessed. First, the three-sigma (3σ) rule is applied for outlier detection and removal. Specifically, for each type of data, its mean μ and standard deviation σ are calculated, and then the data points deviating from the mean by more than three standard deviations (i.e., falling outside the range [μ-3σ, μ+3σ]) are marked as outliers and removed. This statistical distribution-based outlier detection method assumes that the data approximately follows a normal distribution, which is suitable for most engineering measurement data. For specific parameters known not to conform to the normal distribution, the system uses an improved interquartile range (IQR) method or a density-based local outlier factor (LOF) algorithm for anomaly detection to improve the accuracy and applicability of the detection.

[0097] After removing outliers, the remaining data is standardized to eliminate the dimensional differences between different parameters. Standardization uses the Z-score method to convert the original data x to standardized data z: , so that the converted data has a standard normal distribution characteristic with a mean of 0 and a standard deviation of 1. This processing allows parameters of different dimensions and orders of magnitude to be compared and calculated on the same numerical scale, providing a unified data basis for subsequent modeling. For some parameters that need to maintain physical meaning, the system uses the Min-Max normalization method to linearly map the data to the [0, 1] interval: , preserving the relative relationships between parameters.

[0098] In addition to outlier processing and standardization, missing data in the data also needs to be handled. For short-term data missing (less than 3 times the sampling period), the system uses linear interpolation to fill in; for medium-length missing (3-10 sampling periods), spline interpolation is used to ensure data smoothness; for long-term missing (more than 10 sampling periods), the system constructs a temporary prediction model based on historical data or similar working condition data to estimate and fill in, and marks these filled data, giving them lower confidence weights in subsequent analysis. Through this series of processing, the original heterogeneous, incomplete, and inconsistent data is transformed into high-quality, standardized preprocessed data sets, laying a solid foundation for neural ordinary differential equation modeling.

[0099] S1.2: Based on the preprocessed operating data, the energy system state evolution process is represented as a differential equation, a neural network with a residual network architecture is used to approximate the differential equation function and calculate the gradient through the adjoint sensitivity method, and a battery energy storage model, a gas turbine model, and a thermal energy storage model are established; In step S1.2, a continuous-time dynamic system model is constructed based on the theory of neural ordinary differential equations, and the state evolution process of the energy system is represented as a differential equation form, where x(t) is a system state vector, u(t) is a control input vector, and θ is a neural network parameter. The neural network approximation function f using the ResNet architecture is used to calculate the gradient by the adjoint sensitivity method to achieve efficient optimization of the parameters. An independent sub-model is established for each type of resource, the battery energy storage model considers the dynamic change of the state of charge and the capacity attenuation effect, the gas turbine model includes the start-stop characteristics and the influence of thermal inertia, and the thermal energy storage model integrates the heat transfer process and the phase change characteristics.

[0100] Specifically, step S1.2 uses preprocessed operation data to represent the state evolution process of the energy system as a differential equation form, and approximates these differential equations by a neural network. The core idea of this step is to integrate traditional physical modeling methods with modern deep learning techniques, maintaining the interpretability of the model while improving the modeling accuracy of complex nonlinear systems. The system first represents the dynamic process of the energy equipment as a differential equation form: where x(t) represents the system state vector (such as the state of charge of energy storage, equipment temperature, pressure, etc.), u(t) ∈ R m represents the control input vector (such as charging and discharging power, valve opening, etc.), θ represents the parameter set to be learned, and f represents the state derivative function, which describes how the system state evolves over time.

[0101] Traditional methods usually construct the function f based on simplified physical models or linearized approximations, making it difficult to accurately describe the nonlinear and time-varying characteristics of complex energy systems. The invention innovatively uses a deep neural network with a residual network (ResNet) architecture to approximate the function f, fully utilizing the advantages of deep learning in nonlinear function fitting. The key feature of the residual network is the introduction of skip connections, which enables the network to learn the state change rather than the absolute state value, which is highly consistent with the idea of describing the state derivative in differential equations. The neural network structure uses a multi-layer perceptron design, with a typical configuration of: input layer (n+m neurons, receiving state and control input) → first hidden layer (256 neurons) → second hidden layer (128 neurons) → third hidden layer (64 neurons) → output layer (n neurons, output state derivative). The hidden layer uses the Swish activation function f(x) = x·sigmoid(βx), where the β parameter is learnable. This function combines the sparse activation characteristics of ReLU and the smoothness of the hyperbolic tangent function, making it suitable for differential equation approximation tasks.

[0102] The training of neural network parameters θ employs an adjoint sensitivity method to compute the gradient, which is derived from optimal control theory and is capable of efficiently computing the sensitivity of a differential equation system to parameters. Specifically, a loss function is constructed: , which measures the deviation of the model predicted trajectory from the measured data. Then the gradient of the loss function to the parameters is computed by solving the adjoint equation

[0103] , where a is the adjoint variable. This method is more suitable for handling differential equation systems than the traditional backpropagation, with high computational efficiency and small memory footprint, and is capable of handling long time series and high-dimensional state space.

[0104] Based on this framework, dynamic models of three typical energy devices are established respectively. The battery energy storage model considers the state of charge (SOC) dynamic variation, internal resistance variation with temperature and cycle number, nonlinear characteristics of charging and discharging efficiency, and capacity degradation effect. The core differential equation of the model is , where η is the efficiency factor (related to SOC, temperature T and power P), is the rated capacity, and the constant 3600 is used for unit conversion. The gas turbine model contains start-stop characteristics, operating inertia, output regulation rate limit, and variation characteristics of thermal-electric conversion efficiency with load, and the core equations include the power dynamic equation and the temperature dynamic equation , where is the power time constant, is the thermal capacity. The thermal energy storage model integrates the dynamic characteristics of heat transfer process, heat loss mechanism and nonlinear characteristics of phase change material, and the core equation is For phase change materials, special treatment is needed for the nonlinear behavior near the phase change temperature, becomes a function of temperature , which presents a sharp peak near the phase change point.

[0105] Each device model is trained by a dedicated neural network, while introducing physical knowledge to guide network learning, such as adding physical constraint terms in the loss function to ensure the satisfaction of basic physical laws such as energy conservation, the second law of thermodynamics, etc. This "physics-informed neural network" method significantly improves the generalization ability and physical reasonableness of the model, making it still maintain reasonable prediction outside the coverage of the training data.

[0106] S1.3: A variable step-size ordinary differential equation solver is used to train the battery energy storage model, the gas turbine model and the thermal energy storage model respectively, to obtain the continuous-time dynamic system model. ​

[0107] In step S1.3, the model training employs a variable-step ordinary differential equation solver to ensure numerical stability and computational efficiency. The model performance is evaluated by comparing the mean squared error between the model predicted values and the actual measured values, and the early stopping method is used to prevent overfitting. The final continuous-time dynamic system model can accurately describe the dynamic characteristics of various energy resources, providing a foundation for subsequent aggregation problem identification and strategy optimization.

[0108] where step S1.3 employs a variable-step ordinary differential equation solver to train and validate the battery energy storage model, gas turbine model, and thermal energy storage model, and finally integrates them into a complete continuous-time dynamic system model. The training of neural ordinary differential equation models is significantly different from traditional deep learning models, requiring special numerical solution techniques.

[0109] Neural ordinary differential equation models are a new modeling paradigm that combines traditional ordinary differential equation theory with modern deep learning techniques, specifically designed for modeling and prediction problems of continuous-time dynamic systems. The core idea of this model is to consider the solution of a differential equation as the continuous limit of a neural network, parameterizing the vector field function of the differential equation through a neural network, thereby achieving accurate modeling of complex dynamic systems. Unlike traditional discrete-time neural networks, neural ordinary differential equations can handle predictions at any time point, not limited by fixed time steps, making them particularly suitable for handling irregular time sampling and multiple time scales in energy systems.

[0110] The basic principles of neural ordinary differential equations are based on the numerical solution theory of initial value problems of ordinary differential equations. At the mathematical level, any continuous-time dynamic system can be described by a system of first-order ordinary differential equations, and neural ordinary differential equations approximate the vector field of this differential equation system by introducing trainable neural networks. This approximation is not a simple function fitting, but rather learning the internal laws of system state evolution while maintaining the continuity of system dynamics. The role of the neural network here is similar to that of empirical formulas or physical laws in traditional modeling, but its expressive power is stronger, capable of capturing complex nonlinear relationships and high-order coupling effects.

[0111] In the context of energy system modeling, neural ordinary differential equations are particularly suitable for describing the dynamic response characteristics of energy devices. Energy devices such as batteries, capacitors, heat exchangers, etc. have obvious dynamic characteristics, and their state changes follow specific physical laws but are influenced by various complex factors. Traditional modeling methods often require a large number of simplifying assumptions, while neural ordinary differential equations can learn these complex relationships through data-driven methods while maintaining physical meaning. For example, the state of charge of a battery not only depends on the charging and discharging current, but also is affected by temperature, aging degree, cycle number, and other factors. These complex relationships are difficult to accurately describe with analytical formulas, but neural ordinary differential equations can establish accurate dynamic models by learning patterns in historical data.

[0112] The training process of neural ordinary differential equations is significantly different from traditional neural networks. Traditional neural network training is based on discrete input-output pairs, while neural ordinary differential equation training needs to consider the continuity of the entire time trajectory. Training data is usually in the form of time series, containing state observations of the system at different times. The training goal is to make the predicted state trajectory as close as possible to the observed trajectory, which requires solving the differential equation to get the predicted trajectory, and then calculating the difference between the prediction and observation. Due to the involvement of differential equation solving, gradient calculation becomes complex, which requires the use of the aforementioned adjoint sensitivity method to efficiently calculate the gradient.

[0113] Numerical solution of the model is a key technical link in the application of neural ordinary differential equations. Due to the highly nonlinear nature of the vector field parameterized by the neural network, direct analytical solution is almost impossible, and numerical integration methods must be used. Common numerical solvers include Euler method, Runge-Kutta method, multi-step method, etc., each with its own range of applications and accuracy characteristics. In practical applications, the choice of solver needs to balance between calculation accuracy and computational efficiency. For real-time application scenarios such as energy systems, computational efficiency is particularly important, so adaptive step solvers are often used to dynamically adjust the integration step size based on the smoothness of the solution.

[0114] Another important feature of neural ordinary differential equation models is their memory efficiency. Traditional recurrent neural networks require storing all intermediate states when processing long sequences, with memory requirements growing linearly with sequence length. However, neural ordinary differential equations only need to store the final state and the adjoint state when calculating the gradient using the adjoint sensitivity method, making the memory requirement independent of sequence length. This allows the model to handle very long time series data. For long time series of energy system monitoring data, this feature has important practical value.

[0115] The neural differential equation model also has good generalization ability in the application of energy complex. Because the model learns the internal dynamics of the system rather than a simple input-output mapping, the trained model can be extended beyond the training time range to predict the evolution of the system for a longer time. This time extrapolation ability is of great significance for long-term planning and predictive control of energy systems. At the same time, the continuous time characteristics of the model enable it to naturally process data with different sampling frequencies, providing an effective way for the fusion of multi-source heterogeneous data.

[0116] An advanced variable step solver is used to automatically adjust the integration step size according to the severity of the state change, ensuring numerical stability and improving computational efficiency. For rigid differential equation systems (such as thermal systems with multi-time scale characteristics), implicit solvers (such as backward Euler method or backward difference formula method) are used; for non-rigid systems, explicit solvers (such as Runge-Kutta method or Adams method) with higher computational efficiency are used.

[0117] The training process uses a batch training strategy, with each batch containing 10-20 time series samples, each with a length of 100-500 time steps. The Adam optimizer is used for optimization, with an initial learning rate of 0.001 and a learning rate decay mechanism introduced, reducing the learning rate to 0.8 times every 50 epochs. To prevent overfitting, the early stopping method is used, stopping training when the loss function on the validation set does not improve for 10 consecutive epochs. At the same time, the L2 regularization term (weight decay coefficient is 0.0001) is introduced to suppress the excessive growth of model parameters. In addition to the main state prediction error, the loss function of the training process also includes a derivative matching term , prompting the model to learn the correct dynamic derivative relationship.

[0118] The training of the battery energy storage model pays special attention to the dynamic response characteristics under different working conditions. The training data includes various charging and discharging modes: constant power mode, impulse response mode, random fluctuation mode, and high-frequency fluctuation mode under actual grid frequency modulation conditions. The model needs to accurately capture the nonlinear efficiency characteristics under different charging and discharging powers, performance changes under different temperature conditions, and voltage response characteristics within different SOC ranges. The verification criteria include an SOC prediction root mean square error (RMSE) of less than 2%, a voltage prediction RMSE of less than 0.05V, and a power response time error of less than 50ms.

[0119] Gas turbine model training focuses on dynamic characteristics during start-up and load change processes. Training data covers cold start, hot start, normal shutdown, emergency shutdown, and steady-state operation and load ramping processes at different load levels such as 25%, 50%, 75%, and 100%. The model needs to accurately predict start-up time, time to reach rated load, fuel consumption rate, exhaust temperature, and thermal efficiency changes at different loads. Verification criteria include power prediction RMSE less than 3% of rated power, thermal efficiency prediction RMSE less than 1 percentage point, and temperature prediction RMSE less than 5°C.

[0120] Thermal energy storage model training focuses on temperature distribution dynamics and heat loss characteristics during charging and discharging processes. For phase change material thermal storage systems, emphasis is placed on capturing nonlinear behavior near the phase change temperature; for stratified hot water thermal storage systems, emphasis is placed on modeling temperature stratification effects and mixing processes. Verification criteria include temperature prediction RMSE less than 2°C, heat metering error less than 5%, and thermal response time error less than 2 minutes.

[0121] After training each sub-model, a unified continuous-time dynamic system model is constructed through ensemble learning. The ensemble uses a multi-model switching strategy to automatically select the most suitable sub-model based on the current system state and operating mode. At the same time, a soft switching mechanism is used to ensure the smoothness of the model switching process and avoid jumps in the prediction results. The final dynamic system model is tested on a validation dataset to ensure accurate state prediction and dynamic response characteristic description under various operating conditions. The model can predict the state evolution trajectory of the system over future time periods, providing a reliable simulation environment and prediction basis for subsequent aggregation problem identification and optimal control.

[0122] In step S2, based on the continuous-time dynamic system model, a static pulse neural network is used to preliminarily identify the energy aggregation problem, and after identifying potential aggregation needs, an adaptive dynamic pulse neural network classifier is activated to accurately classify specific energy aggregation problems, specifically including: S2.1: Obtain normalized state data, use the leaky integrate-and-fire neuron model in the static pulse neural network for pulse coding, and identify four types of aggregation needs: load peak-valley regulation, frequency regulation, voltage support, and emergency response; In step S2.1, the static pulse neural network serves as the first layer identifier, responsible for quickly detecting potential aggregation needs in the system. The network uses the LIF (Leaky Integrate-and-Fire) neuron model, which simulates the integrate-and-fire characteristics of biological neurons, with the membrane potential dynamic equation being where τ is the membrane time constant (usually set to 20 ms), R is the membrane resistance (value 1 MΩ), and I is the input current. When the membrane potential v exceeds the threshold voltage Vth (-50 mV), the neuron generates a spike and resets the membrane potential to the resting potential V rest (-70 mV). The input layer contains 64 neurons that receive normalized system state data, including key indicators such as: load prediction error (relative error range ±20%), renewable energy output fluctuation (15-minute sliding variance), energy storage state of charge (SOC percentage), grid frequency deviation (±0.5 Hz range), and node voltage amplitude (0.95-1.05 normalized). The hidden layer uses a fully connected structure with 128 LIF neurons, and the connection weights are initialized using the Xavier method with a random distribution of synaptic delays set to 1-5 ms. The output layer is set to four neurons corresponding to load peak-valley adjustment, frequency adjustment, voltage support, and emergency response, and uses the Winner-Take-All mechanism to determine the final recognition result.

[0123] The core of step S2.1 is to replace the traditional artificial neural network with a spiking neural network (SNN) that is closer to the working mechanism of the biological nervous system, and to process the dynamic information of the energy system using a time coding mechanism. First, normalized state data from the output of the continuous-time dynamic system model constructed in step S1 is obtained. After preprocessing, these data form a unified normalized representation in the [0, 1] interval, including load prediction error (relative error range ±20%, normalized to [0, 1]), renewable energy output fluctuation (15-minute sliding variance, normalized to [0, 1]), energy storage state of charge (SOC percentage, originally in the [0, 1] interval), grid frequency deviation (±0.5 Hz range, normalized to [0, 1]), and node voltage amplitude (0.95-1.05 normalized value, normalized to [0, 1]) key indicators.

[0124] The static spiking neural network uses a biologically reasonable Leaky Integrate-and-Fire (LIF) neuron model that processes information by simulating the membrane potential dynamics of biological neurons. The membrane potential dynamic equation of the LIF neuron is where τ is the membrane time constant (usually set to 20 ms), representing the time characteristic of the natural decay of the membrane potential; R is the membrane resistance (value 1 MΩ), representing the proportionality coefficient for converting input current to membrane potential; I is the input current, representing the input signal from the previous layer of neurons; v is the membrane potential, representing the activation state of the neuron. When the membrane potential v exceeds the threshold voltage V th (-50 mV), the neuron generates a spike and resets the membrane potential to the resting potential V rest(-70 mV), and then enters a refractory period (about 2 ms) during which it does not respond to any input. This mechanism enables SNNs to naturally handle timing information, making them particularly suitable for dynamic aggregation problem identification in energy systems.

[0125] The network structure of the static spiking neural network adopts a three-layer architecture: input layer, hidden layer, and output layer. The input layer contains 64 LIF neurons, each responsible for receiving a normalized system state parameter. The input data is converted into a spike sequence through rate coding, with larger values resulting in higher spike frequencies. The typical encoding formula is where x is the normalized input value, f min and f max are the minimum and maximum spike frequencies (usually set to 10 Hz and 100 Hz, respectively). The hidden layer uses a fully connected structure and contains 128 LIF neurons. The connection weights are initialized using the Xavier method to ensure stable signal propagation in the network. The synaptic transmission characteristics include synaptic delay and synaptic weight. The synaptic delay is set to a random distribution between 1-5 ms to simulate the signal transmission delay in biological neural systems, enhancing the network's time sensitivity.

[0126] The output layer is set to 4 LIF neurons, each corresponding to one of the four basic aggregation needs: load peak-valley adjustment (for demand-side management), frequency regulation (for grid stability), voltage support (for local grid quality), and emergency response (for emergencies). The network uses the Winner-Take-All (WTA) mechanism to determine the final identification result, i.e., the class represented by the output neuron with the most spikes within the simulation time window (usually 100 ms) is determined as the main aggregation demand of the current system. In practical applications, if the number of spikes of two or more output neurons is close (with a difference of less than 10%), the system will mark multiple possible aggregation demands for further judgment by a subsequent precise classifier.

[0127] The network training uses a supervised learning method combined with Temporal Backpropagation and approximate gradient descent techniques. Since the spike function itself is not differentiable, the system uses a surrogate gradient method with a smooth sigmoid function The gradient of the approximated impulse function enables end-to-end error backpropagation. The training dataset contains 10,000 labeled samples, covering various typical aggregation scenarios, such as morning and evening peak load periods, periods of significant renewable energy fluctuations, and power grid fault recovery periods. The training objective is to maximize the pulse frequency of the target class neurons while minimizing the pulse frequency of non-target class neurons. The final static spiking neural network can quickly identify potential aggregation demands in the system with an accuracy of over 90%, providing preliminary screening results for subsequent accurate classification.

[0128] S2.2: Based on the four types of aggregation demands, activate the adaptive dynamic spiking neural network classifier, and dynamically adjust the network structure according to the novelty of the input data through the adaptive dynamic spiking neural network classifier, automatically increasing neurons to expand network capacity; Step S2.2 activates the adaptive dynamic spiking neural network classifier based on the four basic aggregation demands identified by the static spiking neural network for accurate classification. Unlike static SNN, dynamic SNN uses a growing when required (GWR) mechanism that can dynamically adjust the network structure according to the novelty of the input data. This feature makes it particularly suitable for handling evolving new aggregation demands in energy systems. Dynamic SNN receives the output results of static SNN and combines more detailed system state data to form a high-dimensional feature vector for accurate classification.

[0129] The core of adaptive dynamic SNN is the growing when required mechanism, which is based on competitive learning theory and gradually forms a topological mapping of the input space by dynamically adding neurons and adjusting connection weights. The network initially contains only a small number of neurons (usually 2-3 per class), which gradually expands during the learning process. Each neuron has an associated prototype vector w, representing its position in the input feature space. When a new feature vector x is input into the network, the system calculates the Euclidean distance d(x, w) between x and all existing neuron prototype vectors, and finds the minimum distance corresponding neuron, called the best matching unit (BMU).

[0130] The growing when required decision depends on two key indicators: the minimum distance d bmu and the network activity a t . The network activity is defined as a t = exp(-d bmu ), reflecting the matching degree of the current input and the existing knowledge of the network. When d bmu is greater than the preset novelty threshold θ sim (typically set to 0.8) and the network activity a t is less than the activity threshold θ actWhen the distance between the current input and the prototype vector of the BMU is less than a second threshold (usually set to 0.1), the system determines that the current input represents a new pattern and requires the addition of a new neuron to expand the network capacity. This dual-threshold mechanism can identify true new patterns while avoiding excessive responses to noise and outliers.

[0131] The addition of a new neuron follows specific rules. First, the system needs to determine the location of the new neuron in the network, i.e., its prototype vector. An effective strategy is to set the prototype vector of the new neuron w new as a linear combination of the current input x and the prototype vector of the BMU w bmu : w new = 0.5x + 0.5w bmu . This allows the new neuron to better represent the input region that is currently not well covered. In addition, the system also needs to find the Second Best Matching Unit (SBMU) and establish a lateral connection between the BMU and the SBMU, forming the topology of the input space. Each neuron also maintains an age counter, recording the number of training samples since the last activation, for subsequent invalid neuron cleaning.

[0132] In addition to structural adjustment, the network updates the prototype vectors of existing neurons through the competitive Hebbian learning rule. For the BMU, its prototype vector is updated according to the following rule: where is the learning rate (usually initially set to 0.1 and gradually reduced to 0.01 as training progresses). For the direct topological neighbors of the BMU, a similar but smaller learning rate update rule is applied: where , ensuring smooth changes in the topology.

[0133] As training progresses, there may be neurons in the network that are not active for a long time, these neurons occupy computing resources but no longer provide useful information. The system regularly checks the age counters of all neurons, when a neuron's age exceeds a pre-set threshold (usually 100 training samples) and its topological connections are less than 2, it is removed from the network, releasing resources for new effective neurons. This dynamic structural adjustment mechanism enables the network to efficiently represent the input space while adapting to changes in data distribution.

[0134] The final adaptive dynamic SNN not only can identify known aggregate demand types, but also can detect and learn new emerging patterns, forming continuous learning and adaptation to the operation characteristics of the energy complex. When a new aggregate demand pattern is identified and stable learning, the system will assign it a specific semantic label, and provide an explanation to the administrator through the human-computer interaction interface, realizing the explainability and transparency of the artificial intelligence system.

[0135] In step S2.2, based on the four types of aggregated demand, an adaptive dynamic spiking neural network classifier is activated, and the network structure is dynamically adjusted according to the novelty of the input data by the adaptive dynamic spiking neural network classifier, and the number of neurons is automatically increased to expand the network capacity, specifically including: S2.2.1: Obtain the feature vector of the four types of aggregated demand, calculate the Euclidean distance between the input data and the existing neuron prototype vector, and determine a new mode when the minimum distance exceeds the preset novelty threshold, triggering the network structure adjustment mechanism; In step S2.2.1, the dynamic classifier adopts a Growing When Required (GWR) mechanism based on the self-organizing mapping theory to dynamically adjust the network structure according to the novelty of the input data. Novelty evaluation is performed by calculating the Euclidean distance d bmu between the input vector and the best matching neuron. When d bmu > θ sim (the threshold is set to 0.8) and the network activity a t < θ act (the threshold is set to 0.1), a new neuron is automatically added to expand the network capacity.

[0136] In the embodiments of the present application, in this stage, first, the feature vectors of the four types of basic aggregated demand output from the static spiking neural network are obtained, which contain rich multi-dimensional information. For load peak-valley regulation demand, the features include load prediction curve shape, peak-valley difference ratio, duration estimate, and distribution of load types that can participate in regulation, etc.; for frequency regulation demand, the features include frequency deviation size, change rate, oscillation characteristics, system inertia state, etc.; for voltage support demand, the features include voltage deviation amplitude, reactive power distribution, sensitive node position, fluctuation period, etc.; for emergency response demand, the features include event type identification, severity assessment, impact range prediction, and necessary response time, etc. After standardization processing, these features form high-dimensional feature vectors in a unified format, and the typical dimension is between 20-30.

[0137] The core of novelty detection is to calculate the distance between the input feature vector and the prototype vector of the existing neuron in the network, and to evaluate the similarity between the current input and the learned knowledge. The system uses Euclidean distance as the similarity measure, and the calculation formula is where x is the input feature vector, and w is the neuron prototype vector. To improve computational efficiency, especially in high-dimensional (more than three-dimensional) feature space, a nearest neighbor search algorithm based on KD tree is implemented, reducing the search complexity from O(n) to O(log n), where n is the number of neurons in the network. In practical applications, the system also introduces an early stopping strategy, which immediately stops searching when the distance of the nearest neighbor found is less than the novelty threshold, further optimizing the calculation performance.

[0138] Two conditions need to be satisfied simultaneously to determine whether an input represents a new pattern: the minimum distance exceeds the preset novelty threshold, and the network activity is lower than the activity threshold. The novelty threshold θsim is a key parameter that directly affects the growth rate and generalization ability of the network. A too low threshold will lead to excessive growth of the network, and it is too sensitive to noise and variation; a too high threshold may ignore important new patterns. In this system, θsim adopts an adaptive setting strategy, with an initial value of 0.8, which gradually increases as the network size grows, and the calculation formula is:

[0139] where n is the current number of neurons, and n0 is the reference value (set to 50). This adaptive strategy allows the network to grow rapidly in the early stage to cover the input space, and then adds new neurons more cautiously in the later stage to avoid redundancy and overfitting.

[0140] The network activity at is the second indicator to evaluate the matching degree of the current input to the network, defined as at = exp(-dbmu), with a value range of (0, 1]. When the best matching distance dbmu is smaller, the activity is closer to 1, indicating that the network has stronger expression ability for the current input. The activity threshold θact is usually set to 0.1, indicating that when the network's expression activity for the input is lower than this value, a new neuron is considered to be added. In practical applications, the activity threshold also adopts a dynamic adjustment strategy, setting a lower value (0.05) in the early stage of training to promote the network to quickly cover the input space, and gradually increasing to a stable value (0.15) as the training progresses, so that the network pays more attention to subtle differences.

[0141] When it is determined that the current input represents a new pattern, the network structure adjustment mechanism will be triggered to prepare to add a new neuron at the corresponding position. This data-driven dynamic structure adjustment allows the network to adaptively represent the distribution characteristics of the input space, providing a flexible learning framework for capturing new types of aggregated demand evolving in energy systems. Compared with traditional fixed structure networks, this dynamic growth mechanism greatly improves the model's expression ability and adaptability, especially suitable for processing energy system data with non-stationary characteristics. At the same time, through the double threshold judgment mechanism and adaptive parameter strategy, the system effectively controls the network complexity while maintaining the learning ability, achieving a good balance between expression ability and computational efficiency.

[0142] S2.2.2: Based on the network structure adjustment mechanism, create a new hidden layer neuron in the corresponding output category, use the feature vector of the current input data as the prototype vector of the new neuron, and initialize the connection weights between the new neuron and the input layer and the output layer; In step S2.2.2, the initial weights of the new neuron are determined by the linear interpolation method: wnew = 0.5(w bmu +w smu ), where w bmu and w smu are the weight vectors of the best and the second best matching neurons, respectively.

[0143] When it is determined that a new neuron needs to be added, the first step is to determine which output class this new neuron should be associated with. In the adaptive dynamic spiking neural network of the present invention, each hidden layer neuron is associated with a specific output class, forming a class-specific representation structure. There are two strategies to determine the class association: supervised and semi-supervised. In the supervised training phase, the system directly uses the label information of the training samples to associate the new neuron to the corresponding output class; in the deployment phase or in the case of unlabeled samples, the system adopts a semi-supervised strategy to infer the most likely class based on the class association of the best matching neuron and the current network activation pattern.

[0144] After determining the class association, the system creates a new hidden layer neuron in the corresponding output class. This process includes three key steps: prototype vector initialization, connection weight setting, and topology structure updating. The prototype vector initialization adopts a hybrid strategy, combining the feature vector of the current input data with the prototype vector of the best matching neuron to form the initial representation of the new neuron. The specific formula is wnew = λx + (1-λ)wbmu, where λ is the mixing factor, initially set to 0.5, and gradually increased as the network develops, so that the neurons added later pay more attention to the features of new samples. This hybrid strategy not only preserves the general features of the current class, but also expresses the uniqueness of new samples, ensuring smooth expansion of the network.

[0145] Connection weight setting is a key step in building the functional role of the new neuron in the network. The new neuron needs to establish two types of connections: forward connection with the input layer and horizontal connection with other hidden layer neurons. The forward connection weight is initialized through supervised learning, using a similar hybrid strategy to the prototype vector, combining the optimal weight of the current sample (calculated through temporary backpropagation) and the weight of the best matching neuron, to ensure that the new neuron can produce appropriate activation for similar inputs. The horizontal connection is established through a competitive learning mechanism, the new neuron establishes excitatory connections (positive weights) with spatially adjacent neurons of the same class and inhibitory connections (negative weights) with adjacent neurons of different classes, forming a local competition-global cooperation activation pattern.

[0146] Topology update is a core feature of self-organizing networks, which dynamically adjusts the connection between neurons to form a topological preserving mapping of the input space. After adding a new neuron, the system needs to update the network topology, mainly including two aspects: first, establish a direct connection between the new neuron and the best matching neuron, and set the initial age to 0; second, disconnect the connection between the best matching neuron and its neighbor farthest away, and connect these neurons with the new neuron. This connection reorganization mechanism ensures that the network structure accurately reflects the topological relationship of the input space, helping to form a semantically coherent representation space.

[0147] In addition to the basic mechanism described above, the present application also introduces several innovative designs to enhance the stability and adaptability of the network. First is the probation period mechanism of the new neuron, the newly added neuron enters a probation period of 30 training samples, during which its prototype vector can be quickly adjusted (learning rate increased by 50%), but it does not participate in further network structure reorganization. This mechanism allows the new neuron to fully adapt to its representation area, while avoiding network structure fluctuations caused by immature neurons. Second is the neuron integration mechanism, the system periodically checks the similarity of adjacent neurons, and when the prototype vectors of two neurons are highly similar (Euclidean distance less than 0.2) and belong to the same category, they are merged into one neuron, reducing redundant representation. Finally, the elimination mechanism of inactive neurons, neurons that have not been activated for a long time (more than 200 training samples) will be removed, releasing computing resources.

[0148] Through this dynamic and balanced structure adjustment mechanism, the adaptive dynamic pulse neural network can efficiently represent and adapt to the complex and variable aggregated demand patterns in the energy system, while maintaining computational efficiency and continuously improving representation ability and classification accuracy. This self-organizing learning method is particularly suitable for handling emerging patterns and concept drift problems in the energy system, enabling the system to continuously adapt to changing energy environments and user demands.

[0149] S2.2.3: updating the activation strength of the new neuron in the network through a competitive learning mechanism, determining the optimal response neuron, and finally determining the expanded network capacity.

[0150] Step S2.2.3 realizes the updating of the activation strength of the new neuron in the network through a competitive learning mechanism, determines the optimal response neuron, and thus realizes the effective expansion of the network capacity. Competitive learning is an unsupervised learning method inspired by biological neural systems, whose core idea is that neurons compete with each other to respond to a specific input pattern, and only the "winning" neuron and its neighbors can update their weights, thus forming a self-organizing mapping of the input space. In the adaptive dynamic pulse neural network, the competitive learning mechanism plays a key role in ensuring that the network can efficiently represent the input space and optimize the allocation of computing resources.

[0151] The competitive learning process begins with the calculation of neuron activation strength. When an input feature vector x enters the network, the system calculates the activation strength of each hidden layer neuron, with the basic formula where is the prototype vector of neuron i, and σ is the width parameter of the Gaussian function (affects the receptive field size of the neuron). To enhance the dynamic representation capability of the network, the present invention adopts a context-sensitive activation mechanism, which modifies the activation strength as:

[0152] where is the context modulation coefficient (dynamically adjusted according to the historical activation pattern of the neuron), is the temporal modulation factor (considers the time structure of the input sequence). This enhanced activation calculation enables the network to better capture the temporal dependencies and contextual information in energy systems, improving the accuracy and stability of classification.

[0153] Determining the optimal response neuron is the core of competitive learning. In the traditional Winner-Take-All strategy, only the neuron with the highest activation strength is allowed to respond and update the weights. However, this extreme competitive mechanism easily leads to uneven resource utilization and over-specialization problems. The present invention adopts an improved k-Winners-Take-All strategy, allowing the top k most strongly activated neurons to participate in response and learning, where k is dynamically adjusted according to the network size, usually set to 5%-10% of the total number of hidden layer neurons. The learning strength of each winning neuron is weighted according to its activation strength ranking, ensuring that the dominant neuron gets the most learning opportunities while allowing suboptimal neurons to learn moderately. This soft competition mechanism significantly improves the network's expression efficiency and generalization ability.

[0154] After determining the winning neuron, the system updates its prototype vector and connection weights through the competitive Hebbian learning rule. For the update of the prototype vector, the basic formula is where is the learning rate specific to the neuron. Unlike the traditional fixed learning rate, the present invention implements a multi-level adaptive learning rate mechanism. At the global level, the basic learning rate gradually decreases as training progresses, with an initial value of 0.1, decaying to 90% of the original value every 1000 training samples, until a minimum value of 0.01. At the neuron level, each neuron adjusts the overall learning rate where is the age modulation function (new neurons have a higher learning rate), frequency modulation function (frequent update of neuron learning rate properly reduced to prevent over-specialization). This multi-level adaptive learning mechanism ensures that the network can balance the need for rapid learning of new patterns and maintaining learned knowledge, significantly improving learning efficiency and stability.

[0155] Connection weight update adopts similar competitive rules, but emphasizes more on the optimization of topological structure. For connections between the winning neuron and its direct topological neighbors, the system resets their age to 0, indicating that these connections are "refreshed" by recent activity. For connections that do not participate in activity, their age increases by 1. When the connection age exceeds a pre-set threshold (usually 50-100), the connection is removed, indicating that the corresponding neuron no longer represents the adjacent input region. At the same time, the system also dynamically creates new connections when two non-adjacent neurons simultaneously produce strong responses to a certain input, establishing new connections between them, reflecting the topological structure of the input space. This dynamic connection management mechanism enables the network to continuously optimize its topological structure, accurately reflecting the inherent relationships of the input space.

[0156] Through the above competitive learning mechanism, the adaptive dynamic spiking neural network can efficiently utilize computational resources to form an optimized representation of the input space. The network capacity dynamically expands during the learning process, enabling it to represent both frequently occurring common patterns and rare but important special cases. As training progresses, the network gradually forms a stable representation structure, with different categories of aggregation needs forming clear boundaries in the feature space, providing a solid foundation for subsequent accurate classification. This self-organizing mapping capability makes the system particularly suitable for handling complex pattern recognition problems in energy complexes, enabling automatic discovery and adaptation to potential structures in data without the need for extensive manual feature engineering and prior knowledge.

[0157] S2.3: Update the connection weights of the dynamic spiking neural network classifier through the adaptive spike-timing-dependent plasticity rule, and perform accurate classification of a specific energy aggregation problem through the dynamic spiking neural network classifier.

[0158] In step S2.3, the network connection weights are updated through the adaptive spike-timing-dependent plasticity (Ad-STDP) rule, which adjusts the synaptic strength based on the time difference Δt between pre- and post-spikes, with the update formula being (activity), where is the learning rate (initial value 0.01, decaying with training), is a time-dependent function, taking a double exponential form or with time constant f(activity) is an adaptive function defined as where is an adaptive coefficient (value 0.1). The Ad-STDP rule can adaptively adjust the learning strength according to the historical activity level of neurons, avoiding overlearning and catastrophic forgetting problems by maintaining the activity history record (sliding window length 1000 time steps). The classifier output layer contains 12 neurons, corresponding to the subdivided aggregate problem types: short-term peak regulation (<15 minutes), medium-term peak regulation (15 minutes-4 hours), long-term energy storage (>4 hours), primary frequency regulation, secondary frequency regulation, voltage reactive power support, harmonic governance, island operation support, emergency load reduction, multi-energy coordinated optimization, demand response management, auxiliary service provision, etc., providing accurate guidance for subsequent strategy making.

[0159] Specifically, step S2.3 updates the connection weights of the dynamic pulse neural network by the adaptive spike-timing-dependent plasticity (Ad-STDP) rule, completing accurate classification of specific aggregate problems. STDP is a learning rule inspired by the biological nervous system, which adjusts synaptic strength based on the relative relationship between the pulse occurrence times of pre- and post-neurons, and is an important learning mechanism discovered in neuroscience. In traditional STDP, if the pulse of the pre-neuron causes the post-neuron to generate a pulse (i.e., the pre-pulse occurs before the post-pulse), the synaptic connection is strengthened; otherwise, if the post-neuron pulse occurs before the pre-neuron, the connection is weakened.

[0160] The improved Ad-STDP rule adopted by the application not only considers the pulse timing relationship, but also incorporates neuron activity history information, realizing a more flexible learning process. Specifically, the synaptic weight update formula is

[0161] (activity degree), where is the learning rate (initial value set to 0.01, gradually reduced during training), is a time-dependent function describing the pulse time difference on the weight change, and f(activity degree) is an adaptive function that adjusts the learning strength according to the historical activity level of neurons.

[0162] Time-dependent function The classic double exponential form is adopted: when >0 (pre-pulse earlier than post-pulse), ; when <0 (post-pulse earlier than pre-pulse), . Where A + and A -The intensities of long-term potentiation (LTP) and long-term depression (LTD) are controlled respectively, usually set to 0.5 and 0.55, so that learning is slightly biased towards inhibition, enhancing the selectivity of the network; The time constants, which control the width of the time windows, are usually set to 20 ms, consistent with the typical values of the biological nervous system.

[0163] The adaptive function f(activity) is an important innovation of the invention, defined as where α is the adaptive coefficient (value 0.1), is the historical activity level of the neuron, calculated by exponential moving average: and β is the smoothing factor (value 0.95). This mechanism reduces the learning rate of neurons with high-frequency activity and increases the learning rate of neurons with low-frequency activity, avoiding the dominance of a few neurons in network behavior, promoting more balanced feature expression, while preventing the problem of catastrophic forgetting.

[0164] To further improve classification accuracy, the system adds a supervised learning layer based on the dynamic SNN. This layer contains 12 output neurons, corresponding to more subdivided aggregate problem types: short-term peak regulation (< 15 minutes), medium-term peak regulation (15 minutes-4 hours), long-term energy storage (> 4 hours), primary frequency regulation, secondary frequency regulation, voltage reactive power support, harmonic governance, island operation support, emergency load reduction, multi-energy coordination optimization, demand response management, and auxiliary service provision. The supervised learning adopts the teacher forcing strategy, providing additional stimulation to target category neurons during the training phase to accelerate the formation of correct mapping.

[0165] The training data is processed in an online learning manner, and the system continuously collects new samples from the actual operating environment and updates the network parameters regularly. To evaluate the classification performance, the system uses confusion matrix analysis to calculate the precision, recall, and F1 score of each aggregate demand category. The target accuracy is set to no less than 95%, ensuring that subsequent strategy optimization is based on accurate problem classification. For borderline cases that are difficult to classify, the system uses a soft classification method to output multiple possible categories and their probability distributions for the decision-making module to reference.

[0166] The final accurate classification result not only contains the type label of the aggregate demand, but also includes key attribute information such as response urgency (divided into immediate response, short-term response, and regular response), expected duration (from minute-level to day-level), resource preference index (indicating the suitability of different types of resources), etc. These rich semantic information provides accurate guidance for subsequent resource response characteristic evaluation and strategy optimization, ensuring the selection of the most suitable resource combination and control strategy for specific aggregate demands.

[0167] Through this precise classification, the system can transform raw system state data into aggregated demand descriptions with clear semantics, establishing a bridge between the physical world of energy and the decision-making and control system, laying the foundation for subsequent optimization decisions. As system operating experience accumulates, the classifier's performance will continue to improve, gradually forming a deep understanding of the operating modes of the energy complex, realizing the transformation from data to knowledge.

[0168] In step S3, the multilateral stochastic flow matching method is used to perform alignment analysis on energy data measured at non-equidistant time points. Through measured value spline augmentation and fractional matching, the response characteristics of each flexibility resource are quantitatively evaluated, specifically including: S3.1: Obtain the data distribution characteristics of different types of energy equipment, construct multiple marginal distributions of electric energy storage resources, thermal energy storage resources and gas energy storage resources, perform nonparametric modeling through kernel density estimation method, and generate a smoothed distribution model; In step S3.1, the construction and modeling of multiple marginal distributions first involves constructing multiple marginal distributions, each corresponding to the data distribution characteristics of a specific resource type or time period. Specifically, for electrical energy storage resources (such as lithium batteries and supercapacitors), the marginal distribution includes key dimensions such as power output (kW), response time (ms level), charge / discharge efficiency (85%-95%), and cycle life decay coefficient; for thermal energy storage resources (such as molten salt thermal storage and phase change materials), the marginal distribution covers parameters such as temperature change gradient (°C / min), thermal storage capacity (MWh), heat release rate (MW), and heat loss coefficient; for gas energy storage resources (such as compressed air energy storage), the marginal distribution includes pressure change rate (bar / min), volumetric flow rate (… Indicators such as compression efficiency, etc., are used. Each marginal distribution is modeled nonparametrically using kernel density estimation, avoiding prior assumptions about the data distribution. A Gaussian kernel function is employed.

[0169] The original data points are smoothed, where h is the bandwidth parameter, which directly affects the smoothness and accuracy of the estimation. The bandwidth parameter is determined through k-fold cross-validation (k=5), with the goal of minimizing the mean squared error. This avoids both overfitting and underfitting.

[0170] Specifically, this step first acquires the data distribution characteristics of different types of energy equipment, constructs multiple marginal distributions, and performs nonparametric modeling using the kernel density estimation method to generate a smoothed distribution model. The multi-marginal distribution framework is a key technology for processing high-dimensional heterogeneous energy data. It allows the system to model the characteristics of different energy dimensions separately, and then combine them through a complex dependency structure, effectively avoiding the curse of dimensionality problem faced by directly modeling high-dimensional joint distributions.

[0171] For electrical energy storage resources, the system constructs marginal distributions of multiple key characteristics. The power output distribution describes the output capability of electrical energy storage devices under different load conditions, typically exhibiting a bimodal distribution reflecting the two main operating modes of charging and discharging. The response time distribution characterizes the time characteristics of the device from receiving instructions to reaching the target output, with significant differences between different technical routes: supercapacitors and flywheel energy storage have response times concentrated in the millisecond range (5-50 ms), lithium-ion batteries in the second range (0.5-5 s), and lead-acid batteries and sodium-sulfur batteries in the ten-second range (10-60 s). The charge-discharge efficiency distribution reflects the energy conversion process loss, which is not only related to the type of technology, but also closely related to the operating state: the efficiency is usually lower under high-power working conditions, and the efficiency is higher under shallow charge-discharge cycles. The cycle life decay coefficient distribution describes the degradation law of the capacity of energy storage devices with the use of cycles, showing obvious technical correlation and dependence on use patterns.

[0172] For thermal energy storage resources, the system constructs marginal distributions of key parameters such as temperature change gradient, heat storage capacity, heat release rate, and heat loss coefficient. The temperature change gradient distribution usually exhibits skewness, with a higher density in the lower gradient region (0.5-2℃ / min), reflecting the thermal inertia characteristics of the thermal system. The heat storage capacity distribution often exhibits a multimodal structure corresponding to different scales and technical types of device clusters. The heat release rate distribution is highly dependent on device type and operating mode: water storage heat storage systems usually have uniform distribution, while phase change material heat storage systems exhibit temperature-dependent non-uniform distribution, with the heat release rate reaching a peak near the phase change temperature. The heat loss coefficient distribution reflects the energy retention capability under different insulation technologies and environmental conditions, typically exhibiting a lognormal distribution characteristic.

[0173] The marginal distribution of gas energy storage resources includes key parameters such as pressure change rate, volume flow, and compression efficiency. The pressure change rate distribution reflects the charging and discharging dynamics of the gas storage system, typically exhibiting conditional distribution, i.e., the distribution form depends on the current pressure level and gas storage capacity. The volume flow distribution is directly related to the power output capability of the system, with large-scale compressed air energy storage systems typically operating in the range. The compression efficiency distribution reflects the energy conversion process loss, which is closely related to the compression ratio, temperature, and device characteristics, with a typical distribution range of 65%-85%.

[0174] After obtaining these raw distribution data, the system uses the Kernel Density Estimation (KDE) method for non-parametric modeling. Compared with parametric distribution (such as normal distribution, Weibull distribution, etc.), KDE does not make prior assumptions about the form of data distribution, and can more flexibly adapt to the complex characteristics of various energy equipment. The basic idea of KDE is to regard each data point as the center of the kernel function, and form a smooth probability density estimate by superimposing these kernel functions.

[0175] The system uses a Gaussian kernel function , which has good mathematical properties and computational efficiency. The bandwidth parameter h directly affects the smoothness of the estimate, and is the most critical hyperparameter in KDE. Too small bandwidth leads to an estimate that is too "rough" and overfits noise; too large bandwidth will be over-smoothed, masking the essential structure of the data. The present invention uses an improved Silverman rule to adaptively determine the initial bandwidth , where σ is the sample standard deviation, IQR is the interquartile range, and n is the sample size. On this basis, the bandwidth parameter is further optimized through k-fold cross-validation (k=5), with the goal of minimizing the log-likelihood difference between the estimated distribution and the actual data.

[0176] To handle sparse areas and extreme values in the marginal distribution, the system introduces an adaptive kernel bandwidth technique. A smaller bandwidth is used in data-intensive areas to preserve the detailed structure, and a larger bandwidth is used in data-sparse areas to improve the stability of the estimate. This adaptive bandwidth strategy significantly improves the estimation quality of the edge region of the distribution, and is particularly important for capturing the limit performance and abnormal behavior of energy equipment.

[0177] Through the above method, the system finally obtains a set of high-quality smooth marginal distribution models, which accurately describe the characteristics of various energy equipment in different dimensions. These distribution models not only provide a statistical description of the performance of the equipment, but also reflect the randomness and variability of the response characteristics of the equipment, providing a solid statistical foundation for subsequent precise matching and response characteristic evaluation.

[0178] S3.2: Based on the distribution model after smoothing, the system uses the measured value spline technique to construct a cubic B-spline basis function to process the observation values at irregular time points, generating a continuous time series; Step S3.2 is based on the distribution model after smoothing, and a cubic B-spline basis function is constructed using the measurement value spline technology to process the observation values at irregular time points and generate a continuous time sequence. This step solves the non-equidistant sampling problem commonly encountered in energy system monitoring, realizes reliable conversion from discrete observation points to continuous time functions, and provides a unified time reference framework for subsequent analysis. Measurement splines are an advanced interpolation technology specially designed for irregular sampling data, which takes into account data fitting accuracy and function smoothness, and is particularly suitable for processing multi-source heterogeneous time series data in energy systems.

[0179] Initially, the system needs to obtain the observation values at irregular time points from the distribution model after smoothing and determine the appropriate spline node sequence. The monitoring data of an energy complex usually comes from different systems and different manufacturers' equipment, with inconsistent sampling frequency and time stamp: the power monitoring system may use high-frequency sampling of 50Hz or 100Hz, the heat system usually uses medium-frequency sampling of 1Hz or 0.2Hz, and the gas system may use a sampling rate of 0.1Hz or lower. This time inconsistency makes it difficult to apply traditional uniform grid interpolation methods, requiring special irregular node processing technology.

[0180] The determination of the node sequence directly affects the quality and computational efficiency of the spline interpolation. The present invention uses an adaptive node distribution strategy to set the node density according to the time interval distribution of the observation data, increasing the number of nodes in data-intensive areas and reducing the number of nodes in data-sparse areas. Specifically, the system first calculates the time interval distribution of adjacent observation points to obtain its percentiles (usually 10%, 25%, 50%, 75%, and 90%). Then, according to these statistical characteristics, the node distribution strategy is designed: in the area where the time interval is less than the 25% percentile, the node density is set to one node per original observation point; in the area where the time interval is between the 25% and 75% percentiles, the node density is one node per two observation points; in the area where the time interval is greater than the 75% percentile, the node density is dynamically adjusted according to the interval size, usually not more than 3 times the original observation point interval. This adaptive node strategy ensures the expression accuracy in data-intensive areas and avoids overfitting and waste of computational resources in sparse areas.

[0181] After determining the node sequence, the system constructs a cubic B-spline basis function, which is the core component of spline interpolation. B-spline is a group of polynomial functions with compact support, each B-spline basis function is non-zero in a limited interval and zero in other intervals. This local support characteristic makes B-spline have significant advantages in computational efficiency and numerical stability. The cubic B-spline basis function is expressed as a cubic polynomial in four consecutive node intervals. These basis functions have Continuity, i.e. the continuity of function values and their first and second derivatives at the nodes, guarantees the high smoothness of the interpolation results, which is particularly suitable for representing the continuous change process of physical quantities in energy systems.

[0182] For a given sequence of nodes , the number of cubic B-spline basis functions is , each spanning four adjacent nodes and being zero outside this interval. The computation of the basis functions is done using a recursive definition, starting with the 0thorder B-spline, i.e. the piecewise constant function. Higher order B-splines are then computed using a recursive formula. This recursive definition gives the B-spline basis functions good numerical properties, avoiding the Runge's phenomenon commonly seen in polynomial interpolation of high order.

[0183] After the basis functions are constructed, the system determines the coefficients of the spline function by least squares fitting, so that the spline curve passes through or approximates all observed points. Let the observed data be (tj, yj), j = 1, 2,..., n, where tj is the time point and yj is the corresponding observed value.

[0184] To avoid overfitting and control the smoothness of the curve, the system introduces a regularization term, forming a least squares problem with a penalty. The system uses the Generalized Cross-Validation (GCV) method to automatically select the optimal λ value, where Sλ is the smoothing matrix, and tr(Sλ) is its trace, representing the "equivalent degrees of freedom".

[0185] Through the above measurement value spline technology, the system successfully converts the discrete observed values at irregular time points into a continuous time function representation, solving the key challenge of inconsistent time of energy data. The generated continuous time series has good interpolation accuracy and smoothness, and can accurately reflect the dynamic response process of various energy equipment, providing a unified time reference framework and high-quality data basis for subsequent fractional matching and performance evaluation. This continuous representation is particularly advantageous in capturing dynamic characteristics at different time scales, from millisecond-level power transient response to hour-level thermal system changes, which can be accurately described.

[0186] In step S3.2, based on the smoothed distribution model, the measurement value spline technology is used to construct cubic B-spline basis functions to process the observed values at irregular time points, and a continuous time series is generated, which specifically includes: S3.2.1: Obtain the observed values at irregular time points in the smoothed distribution model, determine the sequence of spline nodes, set the node density according to the time interval distribution of the observed data, and increase the number of nodes in data-intensive areas; Step S3.2.1 details the feature extraction and analysis process for non-uniformly spaced time series data. In the actual operation environment of the energy complex, the data acquisition equipment is often affected by factors such as communication delay, equipment failure, or network congestion, resulting in observation sequences with irregular time intervals. This non-uniform interval characteristic poses a challenge to traditional time series analysis methods, as most standard methods assume that data points are uniformly distributed in time. To address this issue, the system first performs comprehensive feature extraction on the collected raw time series, capturing the basic statistical characteristics, time structure characteristics, and frequency domain characteristics of the data.

[0187] Basic statistical feature extraction is the first step in understanding data distribution. The system calculates the central tendency indicators (mean, median, mode), dispersion indicators (standard deviation, interquartile range, range), distribution shape indicators (skewness, kurtosis), and extreme value characteristics (maximum and minimum values and their occurrence times) for each time series. Unlike uniformly spaced sequences, the calculation of these statistics needs to take into account the non-uniformity of time intervals. For example, the weighted mean calculation formula is where the weight is inversely proportional to the time interval between adjacent observation points, ensuring that long time intervals do not excessively affect the overall statistical characteristics. This time weighting method enables the statistics to more accurately reflect the actual distribution of physical quantities, rather than being affected by the sampling pattern.

[0188] Time structure feature extraction focuses on the dynamic evolution pattern of the data. The system uses the Lomb-Scargle periodogram to analyze the periodic components in non-uniformly spaced data, which does not require equal-interval sampling and is particularly suitable for handling irregular observation data in energy systems. The Lomb-Scargle method calculates the normalized power spectrum Through this method, the system can identify the main periodic components in the data, even in the case of irregular sampling intervals, and obtain reliable results.

[0189] In addition to periodic analysis, the system also detects trend components and change point features of time series. Trend detection uses the Mann-Kendall non-parametric test, which does not require equal-interval distribution of data, only focuses on the relative size relationship of data points, and calculates the statistic S. Frequency domain feature extraction converts time series into the frequency domain for analysis, revealing the contribution of different frequency components. Since traditional fast Fourier transform (FFT) requires equal-interval sampling, the system uses non-uniform FFT (NUFFT) technology to process non-equal-interval data. NUFFT maps non-equal-interval data to a uniform grid through interpolation and convolution methods, then applies the standard FFT algorithm, and finally performs the corresponding compensation correction. This method retains the computational efficiency of FFT while accommodating the characteristics of non-equal-interval sampling. From the frequency spectrum, the system extracts features such as main frequency components, frequency band energy distribution, and power spectral density patterns, which are of great value for identifying the oscillation characteristics, control response characteristics, and noise characteristics of energy equipment.

[0190] For different types of energy equipment, the system also extracts specific domain features. For electrical energy storage equipment, features such as charge and discharge state transition frequency, maximum power change rate, and response delay distribution are extracted; for thermal energy storage equipment, features such as thermal response curve shape, temperature gradient distribution, and thermal inertia characteristics are extracted; for gas energy storage equipment, features such as pressure fluctuation characteristics, flow variation patterns, and compression response curves are extracted. These domain-specific features fully consider the physical characteristics and technical constraints of different energy forms, providing a rich feature basis for subsequent resource response characteristic evaluation.

[0191] Through the above comprehensive feature extraction process, the system converts non-equal-interval time series data into structured feature vectors, which capture various aspects of energy equipment response characteristics, providing rich and accurate information input for subsequent measurement value spline processing and fractional matching algorithms, effectively overcoming the analysis challenges brought by the time irregularity of raw data.

[0192] S3.2.2: Based on the node sequence, construct cubic B-spline basis functions, each of which is a cubic polynomial within four consecutive node intervals and zero in other intervals; The detailed description of step S3.2.2 is as follows: The quality of spline interpolation is highly dependent on the selection of node sequences, and improper node distribution can lead to problems of over-smoothing (too few nodes) or over-fitting (too many nodes). This step uses a multi-level adaptive strategy to intelligently determine the optimal node sequence based on the time distribution characteristics, value range variation characteristics, and physical model constraints of the data, ensuring the best balance between computational efficiency and fitting accuracy.

[0193] The system first performs preliminary node distribution analysis based on the characteristics of the original non-equidistant time series. The core idea is to increase node density in areas with dramatic data changes and reduce node density in relatively stable areas, thereby achieving the best fitting effect with the least number of nodes. In specific implementation, the system calculates the time interval and the value change between adjacent data points, and then calculates the change rate . Based on these basic quantities, the system constructs an initial node distribution strategy: for areas with a change rate exceeding the 90th percentile (dramatic change area), a node is set for each original data point; for areas with a change rate between the 50th and 90th percentiles (moderate change area), a node is set for every 2-3 data points; for areas with a change rate below the 50th percentile (stable area), a node is set for every 5-10 data points. This data-driven adaptive strategy ensures a high match between node distribution and data change characteristics.

[0194] In addition to the basic strategy based on change rate, the system also considers the distribution characteristics of time interval. For areas with a particularly large time interval (greater than the 95th percentile), which may represent data interruption or device offline period, the system increases special processing logic: if the data changes significantly on both sides of the large interval (more than 20% of the data range), a transition node is added on each side of the interval; if the change is not significant, a small number of nodes are uniformly inserted in the large interval area to avoid "information vacuum area". This special processing for abnormal time intervals significantly improves the system's adaptability to data interruption.

[0195] For energy devices with obvious physical model constraints, the system introduces a model-guided node optimization strategy. For example, for a battery energy storage system, the state of charge (SOC) change usually follows a specific charge-discharge curve, with a significant change in slope when SOC approaches 0% or 100%. The system uses this prior knowledge to automatically increase node density in the SOC critical region (0-10% and 90-100%), ensuring the physical reasonableness of the fitted curve even if the original data is insufficiently sampled in these areas. Similarly, for a thermal energy storage system, node density is increased near the phase change temperature (±2°C range) to accurately capture the nonlinear thermal characteristics during phase change; for a compressed air energy storage system, node density is increased in the region where pressure approaches the maximum rated value to reflect the nonlinear compression characteristics in the high-pressure state.

[0196] After the initial node sequence is determined, the system further applies an adaptive node optimization algorithm to optimize the node distribution through an iterative process. The core idea is to identify areas where the current node distribution is insufficient through local fitting error analysis, and to increase or adjust the node positions accordingly. Specifically, the system first uses the initial node sequence for spline fitting, then calculates the fitting error of each original data point where s(t) is the spline function. Based on the error distribution, error hotspots (error greater than twice the average error) are identified, and nodes are added in these areas. The number of added nodes is proportional to the local error size, but an upper limit (usually 5 new nodes) is set to prevent overfitting. The precise location of the new nodes is determined by minimizing the local error function, ensuring the maximum reduction of fitting error.

[0197] To prevent unlimited growth of the number of nodes, the system also introduces a node redundancy detection and merging mechanism. When the spline function between two adjacent nodes is almost linear (the second derivative is close to zero) and the fitting error is small, one of the nodes can be removed without significantly affecting the fitting quality. The system calculates a node importance indicator where is a distance weighting function that evaluates the contribution of each node. Nodes with importance below a threshold are marked as redundant and can be removed while maintaining overall fitting accuracy. This dynamic balancing mechanism ensures that the number of nodes is maintained within a reasonable range, avoiding underfitting while preventing overfitting and wasting computational resources.

[0198] Through the above multi-level adaptive strategy, the system ultimately obtains a highly optimized spline node sequence that matches the time structure, value range changes and physical characteristics of the data in distribution, laying a solid foundation for the construction of cubic B-spline basis functions in the next step. Compared with fixed interval or simple heuristic methods, this intelligent node determination strategy significantly improves the accuracy and computational efficiency of spline interpolation, especially suitable for handling complex and variable non-equidistant time series data in energy systems.

[0199] S3.2.3: Determine the coefficients of the cubic B-spline basis function by the least squares fitting method, so that the spline curve passes through or approximates all observation points, obtaining the continuous time series with time continuity and smoothness.

[0200] Step S3.2.3 completes the construction of cubic B-spline basis functions based on the optimized node sequence, and applies the regularization least squares method to determine the spline function coefficients, finally generating high-quality continuous time series representation. B-spline is a class of piecewise polynomial functions with local support, which has the advantages of high computational efficiency, good numerical stability and strong flexibility compared with traditional interpolation polynomials, especially suitable for handling large-scale complex energy system data.

[0201] The system first determines the optimized node sequence based on step S3.2.2 A set of cubic B-spline basis functions is constructed. Cubic B-spline refers to B-spline of order k = 3, which is a cubic polynomial in each segment, with C 2 continuity at the nodes (i.e. both the function value and its first and second derivatives are continuous). This smoothing property makes cubic B-spline particularly suitable for representing continuous state changes in physical systems, accurately capturing the dynamic response characteristics of energy equipment while avoiding unnatural oscillations or jumps.

[0202] S3.3: Calculate the difference between the minimized real data distribution of the continuous time sequence and the score function of the distribution model using the score matching algorithm to obtain the response characteristics of each flexibility resource.

[0203] Step S3.3 calculates the difference between the real data distribution of the continuous time sequence and the score function of the distribution model using the score matching algorithm to quantitatively evaluate the response characteristics of each flexibility resource. Score matching is an advanced distribution estimation technique that is particularly suitable for handling density estimation problems in high-dimensional data space. By comparing the score function (gradient of the log probability density function) of the distribution instead of directly comparing the probability density, it skillfully avoids the difficulty of calculating the normalization constant in traditional methods, greatly improving the computational efficiency and numerical stability.

[0204] First, the system obtains the real data samples of the continuous time sequence generated in the previous step, which have been converted to continuous representation under a unified time reference after measurement value spline processing. Based on these samples, the system calculates the log density gradient of the real data distribution pdata(x) as the real score function. Specifically, for each data point x, its score function is calculated by the gradient of kernel density estimation: where K is the kernel function, and its derivative. This non-parametric estimation method avoids prior assumptions about the form of the data distribution, and can adapt to the complex characteristic distribution of various energy equipment.

[0205] At the same time, the system establishes a parameterized model distribution and calculates its log density gradient , as the score function of the distribution model. The choice of the model distribution depends on the characteristics of the problem and the computational complexity. For simpler unimodal distributions, the system adopts exponential family distributions (e.g., normal distribution, gamma distribution, beta distribution, etc.) and their mixture models; for complex multimodal distributions or distributions with special constraints, the system adopts more flexible representations such as normalizing Flows or Energy-Based Models. Regardless of the choice of model, its score function can usually be calculated directly through analytical derivatives or automatic differentiation techniques, avoiding the difficulty of directly calculating the probability density and its normalization constant.

[0206] Based on the true score function and the score function of the distribution model, the system constructs a score matching objective function, calculates the squared error between the two score functions using the least squares method, and integrates the calculation results. Mathematically represented as , where Epdata represents the expectation of the true data distribution. This objective function can be further simplified by integral division, avoiding the direct estimation of , where represents the trace of the log density function Hessian matrix, i.e., the sum of the second-order derivatives of each dimension. This conversion makes score matching only need to calculate the derivative of the model distribution, without estimating the derivative of the true distribution, further improving the computational efficiency and stability.

[0207] The optimization of the score matching objective function adopts gradient descent type algorithms, among which the Adam optimizer is widely used due to its ability to adaptively adjust the learning rate for different parameters. The parameter settings of the optimization process are: the initial learning rate is 0.001, the momentum parameters β1=0.9, β2=0.999, and the gradient clipping threshold is 10.0 to prevent gradient explosion. To prevent overfitting, the system uses early stopping method (early stopping), which stops training when the loss function on the validation set does not improve for 30 consecutive epochs. The optimization process usually requires 500-1000 iterations, and the convergence criterion is that the relative change of the loss function is less than 10 -6 .

[0208] Through score matching optimization, the system obtains the optimal model parameters θ*, which makes the score function of the distribution model most approximate to the true score function. Based on these optimized parameters, the system quantitatively evaluates the response characteristics of each flexibility resource, mainly including three key indicators: 1. Response time constant: Characterizes the time behavior of a resource from receiving an instruction to reaching the target output. For electrical energy storage devices, the response time constant is usually in the range of milliseconds to seconds, reflecting its fast regulation capability; for thermal energy storage devices, the response time constant is in the range of minutes to hours, reflecting the thermal inertia characteristics of the thermal system; for gas energy storage devices, the response time constant is in the range of seconds to minutes, reflecting the dynamic characteristics of the compression and release process.

[0209] 2. Regulation accuracy: Characterizes the degree of matching between the actual output of the resource and the target instruction. Regulation accuracy is usually expressed in terms of relative error or root mean square error, reflecting the accuracy and stability of the control system. High-precision regulation (error <1%) is suitable for fine grid frequency modulation services, medium-precision regulation (error 1%-5%) is suitable for conventional load tracking, and low-precision regulation (error >5%) is only suitable for extensive peak shaving.

[0210] 3. Stability index: Characterizes the volatility and sustainability of the resource output. The stability index includes multiple dimensions such as output variance, maximum duration, decay rate, etc., comprehensively reflecting the reliability and durability of the resource at different time scales. High-stability resources are suitable for providing long-term basic support services, medium-stability resources are suitable for medium- and short-term regulation services, and low-stability resources are only suitable for short-term pulse response services.

[0211] These response characteristic indicators constitute the capability profile of each flexible resource, forming a response characteristic matrix that accurately describes the state transition randomness and decision sequence of each resource in the energy aggregation process. These indicators not only reflect the physical characteristics and technical limitations of the resource, but also include the comprehensive effects of factors such as the reaction speed of the control system, communication delay, measurement accuracy, etc., providing a comprehensive and accurate description of the resource characteristics for subsequent reinforcement learning strategy optimization, ensuring that the optimization decision can fully consider the actual capabilities and coordination potential of each resource, and achieve the optimization of the overall performance of the system.

[0212] S3.3.1: Obtain the real data distribution of the continuous time sequence, calculate the logarithmic density gradient of the real data distribution and take it as the real score function, and calculate the logarithmic density gradient of the distribution model and take it as the score function of the distribution model; In step S3.3.1, the score matching algorithm effectively solves the problem of dimension disaster and numerical instability in high-dimensional data space by skillfully avoiding the calculation of the normalization constant. The core idea of this algorithm is to minimize the difference between the score function (gradient of the logarithmic probability density function) of the real data distribution and the distribution model.

[0213] In the first key step of the fractional matching algorithm, the system needs to obtain real data samples from continuous time series and calculate the log-density gradient of the real data distribution based on these samples. This process first extracts representative sample points from the high-quality continuous time series generated in the previous step, which not only contains the original observation data but also the estimated values of the intermediate time points generated by the spline interpolation technique. The sample selection strategy considers the uniformity of the time distribution and the representativeness of the data changes, ensuring that both typical running states and various abnormal and transitional states are covered.

[0214] The calculation of the real fractional function adopts a non-parametric method based on kernel density estimation. The system first constructs a Gaussian kernel function for each data point, and then estimates the log-density gradient of the overall distribution by calculating the weighted gradient of all kernel functions. The advantage of this method is that it does not require any prior assumptions about the form of the data distribution, and can automatically adapt to various complex distribution patterns. In actual calculation, the system adopts an adaptive bandwidth selection strategy to dynamically adjust the width parameter of the kernel function according to the local density characteristics of the data, ensuring fine gradient estimation in data-intensive areas and stability of the estimation in data-sparse areas.

[0215] At the same time, the system establishes a parameterized model distribution to approximate the real data distribution. The selection of the model distribution is based on the physical characteristics of the energy system and the statistical characteristics of the data. For simple distributions that exhibit unimodal characteristics, the system adopts appropriate members of the exponential family distribution, such as normal distribution for symmetric distribution, gamma distribution for right-skewed distribution, beta distribution for bounded distribution, etc. For distributions that exhibit multi-modal characteristics or complex shapes, the system adopts mixed models or more flexible energy models for representation. The log-density gradient of the model distribution can usually be calculated directly through analytical methods or automatic differentiation techniques, which greatly simplifies the subsequent optimization process.

[0216] When dealing with high-dimensional data, the calculation of the fractional function faces the challenge of dimensionality curse. The system adopts multiple dimension reduction and simplification strategies to address this problem. First, the main variation direction of the data is identified through principal component analysis or other dimension reduction techniques, and the high-dimensional problem is projected into a low-dimensional subspace for processing. Secondly, the decomposition strategy is adopted to decompose the high-dimensional joint distribution into a combination of multiple low-dimensional marginal distributions and conditional distributions, and then the fractional function is estimated and synthesized after decomposition. This decomposition not only reduces the computational complexity, but also improves the stability and interpretability of the estimation.

[0217] S3.3.2: Based on the real fractional function and the fractional function of the distribution model, construct a fractional matching objective function, calculate the squared error between the real fractional function and the fractional function of the distribution model using the least squares method, and integrate and sum the calculation results; In step S3.3.2, the objective function is defined as where is the real data distribution, is the parametric model distribution, denotes the score function. In practical computation, the real score function is obtained by empirical estimation: where is the sample point.

[0218] Based on the real score function and the score function of the distribution model calculated in the previous step, the system enters the construction phase of the score matching objective function. The core idea of the design of this objective function is to minimize the squared error between the two score functions, and in this way the gradient characteristics of the model distribution are as close as possible to the gradient characteristics of the real data distribution. Compared with the traditional maximum likelihood estimation method, the score matching method skillfully avoids the calculation difficulty of the probability density function normalization constant, which has a significant computational advantage in dealing with complex high-dimensional distributions.

[0219] The construction of the objective function adopts the square loss of the expectation form, that is, the square expectation of the difference between the two score functions under the real data distribution. In actual calculation, this expectation is approximated by the average of a limited number of samples. In order to improve the accuracy of the estimation, the system uses advanced sampling techniques such as importance sampling and stratified sampling. Importance sampling reduces the estimation variance by adjusting the sample weight, which is particularly effective in dealing with tail distributions. Stratified sampling ensures uniform coverage of samples in the entire data space, avoiding insufficient samples in some important areas.

[0220] In the specific implementation of the objective function, the system also considers the balance between computational efficiency and numerical stability. For large-scale data sets, direct calculation of the score difference of all sample pairs is computationally infeasible. The system uses the mini-batch stochastic gradient method, which only calculates the objective function value of a subset of samples each time, and approximates the global optimal solution through multiple iterations. The size of the mini-batch needs to be carefully selected, as too small a batch size will result in excessive gradient estimation noise, and too large a batch size will reduce computational efficiency and memory usage efficiency.

[0221] In order to enhance the stability of the optimization process, the system introduces regularization terms in the objective function. These regularization terms include the L2 norm penalty of the parameters, which is used to prevent model overfitting; the smoothness constraint of the gradient, which is used to ensure that the learned score function has good mathematical properties; and the penalty term based on physical constraints, which is used to ensure the reasonableness of the model distribution in physics. The selection of regularization parameters is automatically determined by methods such as cross-validation, achieving the best balance between fitting accuracy and generalization ability.

[0222] The objective function also needs to handle the problem of scale difference between different dimensions and different physical quantities. Different parameters in the energy system often have different orders of magnitude and units, and direct error calculation may lead to some dimensions dominating the optimization process. The system adopts an adaptive weight adjustment strategy, dynamically adjusts the weight coefficient according to the variance and importance of each dimension. At the same time, a multi-scale objective function is introduced, and sub-objectives are defined at different time and spatial scales. Through multi-objective optimization technology, the global optimal solution is found.

[0223] S3.3.3: Optimize the model parameters of the score matching objective function through the gradient descent algorithm, so that the score function of the distribution model approximates the true score function, and obtain the response characteristics of each flexible resource containing response time constant, adjustment precision and stability index.

[0224] In step S3.3.3, the Adam optimizer is used to realize the random gradient descent method to optimize the parameters θ, and the learning rate is set to 0.001, and the momentum parameters β1=0.9, β2=0.999. Through multiple iterations (usually 500-1000 rounds), the accurate matching of each marginal distribution is realized, and the convergence criterion is that the loss function changes less than 10 -6 .

[0225] Specifically, after the objective function is constructed, the system enters the iterative optimization phase of the model parameters. This phase adopts a gradient-based optimization algorithm, among which the Adam optimizer is the first choice due to its ability to adaptively adjust the learning rate of different parameters. The Adam optimizer combines the advantages of the momentum method and the RMSprop algorithm, and can automatically adjust the learning rate of each parameter during training. It is particularly effective for problems with sparse or noisy gradients. In such a complex optimization problem as the energy system, the gradient information of different parameters often has great differences, and the adaptive characteristics of Adam can ensure that all parameters are properly updated.

[0226] The hyperparameter settings of the optimization process need to be tuned according to the specific problem characteristics. The initial value of the learning rate is set to 0.001, which is a widely validated default value that neither leads to divergence of the optimization process nor makes the convergence speed too slow. The momentum parameter is set to 0.9, which controls the influence of historical gradient information, and a higher value helps to accelerate convergence and reduce oscillation. The decay rate of the second moment estimate is set to 0.999, which is used to estimate the variance of the gradient, and a higher value helps to maintain stable parameter updates in the later training period. To prevent the problem of gradient explosion, the system sets the gradient clipping threshold to 10.0, and when the norm of the gradient exceeds this threshold, it will be scaled proportionally.

[0227] The convergence judgment of the training process adopts multiple criteria. First, the trend of the loss function, when the relative change of the loss function of the last ten iterations is less than one ten-thousandth, it is considered to have approached the convergence state. Second, the amplitude of parameter change, when the update amplitude of all parameters is less than one thousandth of their current values, it indicates that the parameters have stabilized. Finally, the change of the validation set performance, by monitoring the performance of the model on the independent validation set to prevent overfitting. The early stopping mechanism is set to stop training when the validation set loss has not improved for thirty consecutive periods, which helps to select the model with the best generalization performance.

[0228] The optimization process usually requires hundreds to thousands of iterations to reach a satisfactory convergence state. During the training process, the system monitors the trend of various indicators in real time, including training loss, validation loss, gradient norm, parameter norm, etc. By visualizing the change curve of these indicators, problems in the training process can be found in time, such as gradient vanishing, gradient explosion, overfitting, etc., and appropriate adjustment measures can be taken. At the same time, the system also records key intermediate results for subsequent analysis and debugging.

[0229] After the model training is completed, the system extracts the response characteristic indicators of various flexibility resources from the optimal parameters obtained by optimization. The extraction of response time constant is based on the time decay characteristics of the model distribution, and the change rate of the fractional function in the time dimension is analyzed to determine it. For electric energy storage devices, the response time constant is usually in the range of milliseconds to seconds, reflecting its rapid regulation capability; for thermal energy storage devices, the response time constant is in the range of minutes to hours, embodying the thermal inertia characteristics of the thermal system; for gas energy storage devices, the response time constant is in the range of seconds to minutes, reflecting the dynamic characteristics of the compression and release process.

[0230] The regulation accuracy indicator is calculated by analyzing the deviation distribution between the model predicted value and the actual observed value. The system calculates various accuracy indicators such as mean absolute error, root mean square error, maximum error, etc., and determines the corresponding accuracy level according to the requirements of different application scenarios. High-precision regulation capability makes the resource suitable for participating in fine grid auxiliary services, medium-precision regulation capability is suitable for conventional load tracking tasks, and low-precision regulation capability is mainly used for extensive peak shaving and valley filling applications.

[0231] The calculation of stability indicators covers multiple time scales and multiple performance dimensions. Short-term stability is evaluated by analyzing the high-frequency fluctuation characteristics of resource output, including output variance, fluctuation frequency, decay rate, etc. Medium-term stability focuses on the performance retention ability of the resource during continuous operation, including maximum continuous operation time, performance decay rate, etc. Long-term stability evaluates the reliability and life characteristics of the resource in the long-term use process. These stability indicators provide an important basis for resource scheduling strategy formulation and maintenance plan arrangement.

[0232] In step S4, a twin-delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregation policy optimization model, specifically including: S4.1: model the energy aggregation process as a Markov decision process, construct a state space including load demand, renewable energy output and real-time state of resources, and an action space including resource combination and output allocation ratio, design a reward function including response speed, adjustment accuracy and operation cost, based on the Markov decision process framework, state space, action space and reward function, use a twin-delayed deep deterministic policy gradient algorithm for training and optimization to generate a reinforcement learning aggregation policy optimization model.

[0233] In step S4.1, based on the Markov decision process framework, state space, action space and reward function, a twin-delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregation policy optimization model, specifically including: S4.1.1: obtain multi-dimensional characteristic data, construct a state space including load demand dimension, renewable energy dimension, energy storage state dimension, power quality dimension and economic signal dimension, and an action space including resource participation allocation vector and aggregation mode selection, and establish a Markov decision process framework; In step S4.1.1, the energy complex dynamic aggregation problem is formalized as a Markov decision process (MDP), defined as a five-tuple, where S is the state space, A is the action space, P is the state transition probability, R is the reward function, and γ is the discount factor. The construction of the state space S fully considers the multi-dimensional characteristics and time-varying nature of the energy system. Specifically, it includes: load demand dimension (current total load demand, load prediction uncertainty, load change rate), renewable energy dimension (wind power output prediction, photovoltaic output prediction, prediction confidence interval, weather state coding), energy storage state dimension (electric energy storage state of charge, thermal energy storage temperature distribution, gas energy storage pressure state, available capacity of each energy storage device), power quality dimension (system frequency, key node voltage amplitude, power factor, harmonic distortion rate), and economic signal dimension (real-time electricity price, auxiliary service price, carbon emission factor). The action space A is designed as a mixed continuous-discrete action space, with the continuous action part defined as the participation allocation vector of various resources, and the discrete action part representing the aggregation mode selection.

[0234] S4.1.2: in the Markov decision process framework, design a multi-objective reward function including response speed reward item, adjustment accuracy reward item, operation cost reward item, system stability reward item and environmental benefit reward item, and generate a comprehensive performance evaluation index; In step S4.1.2, the reward function is in the form of weighted multi-objective, which comprehensively evaluates the performance of the aggregation policy: response speed reward item:

[0235] wherein measuring response time deviation, heavy penalty for timeout.

[0236] Adjusting precision reward term:

[0237] The first term evaluates overall precision, and the second term considers individual deviation of each resource.

[0238] Operating cost reward term:

[0239] Including energy cost, wear and tear cost and start-stop cost.

[0240] System stability reward term:

[0241] Evaluate the effect of frequency stability, voltage stability and power oscillation suppression.

[0242] Environmental benefit reward term:

[0243] Encourage renewable energy consumption and punish fossil energy use.

[0244] Total reward function:

[0245] Weight parameters are determined by AHP: .

[0246] S4.1.3: Based on the comprehensive performance evaluation index, an improved twin delayed deep deterministic policy gradient algorithm is used to construct actor network and critic network to obtain the reinforcement learning aggregated policy optimization model.

[0247] In step S4.1.3, an improved twin delayed deep deterministic policy gradient (Enhanced TD3, ETD3) algorithm is used to optimize the design according to the characteristics of the energy aggregation problem. The Actor network structure: the input layer receives a 17-dimensional state vector, which is normalized through Layer normalization The first hidden layer has 256 neurons, and uses Swish activation function , β = 1.702. The second hidden layer has 128 neurons, and a residual connection is introduced to alleviate the vanishing gradient problem. The third hidden layer has 64 neurons, and a Dropout layer (p = 0.1) is added to prevent overfitting. The output layer is divided into two parts: the continuous action output adopts a Softmax activation to ensure the discrete action output adopts a Gumbel-Softmax to realize differentiable sampling. The Critic network structure: adopts a double-path design, and the state path and the action path are processed respectively and then fused. The state path: 17-dimensional input → 128-dimensional hidden layer → 64-dimensional representation. The action path: n+1-dimensional input → 64-dimensional hidden layer → 64-dimensional representation. The fusion layer: 128-dimensional concatenation → 64-dimensional → 32-dimensional → 1-dimensional Q value output. Two Critic networks Q1 and Q2 have the same structure but are independently initialized.

[0248] The experience replay buffer adopts a priority sampling mechanism, and the sampling probability is allocated according to the TD error . , where α = 0.6 controls the priority intensity. The buffer capacity is set to 2 × 10 6 , to ensure sufficient experience diversity. The target network update adopts an adaptive soft update mechanism: , where τ base = 0.001, τ max = 0.01, and λ = 0.001, to achieve a balance between fast learning in the early stage of training and stable convergence in the later stage. The policy update introduces a delay mechanism and noise regularization: the Actor network is updated every d = 3 steps, and the target action adds truncated Gaussian noise , -c, c), where σ = 0.1 and c = 0.2, to improve the robustness of the policy. The loss function optimization: the Critic loss adopts Huber loss or , to reduce the influence of outliers. The Actor loss adds an entropy regularization term H(π) = -Σπ(a|s)logπ(a|s), to encourage policy exploration. The learning rate scheduling: adopts a cosine annealing strategy , where , , and T is the total number of training steps, to realize periodic adjustment of the learning rate. When the average reward of consecutive 100 episodes changes by less than 0.1% and the policy loss stabilizes at 10 -4 or less, it is considered that the training converges. The average reward of the final model on the validation set is not less than -50, and the aggregate accuracy is more than 95%.

[0249] In step S5, the state data is received in real time by using the reinforcement learning aggregation policy optimization model, and the optimal aggregation policy is output, and the optimal aggregation policy is used to send adjustment instructions to the flexible resources through the distributed control system, and actual response data is collected for online learning and updating, to realize adaptive dynamic optimization, specifically including: S5.1: Obtain real-time state data by using the reinforcement learning aggregation policy optimization model, verify the feasibility of the policy through the constraint checking module, if the policy violates the constraint, project the action vector into the feasible region by using the projection algorithm, and generate the optimal aggregation policy; In step S5.1, the trained reinforcement learning model is deployed to the edge computing node of the energy management system, and a communication connection is established with each resource controller through high-speed Ethernet. When the system is running, the model receives real-time state data at a period of 1 second, including device operating parameters collected by the SCADA system, weather prediction data, load prediction results, etc. After data preprocessing and feature extraction, the Actor network is input, and the network outputs a continuous action vector representing the optimal participation ratio of each type of resource. The feasibility of the policy is verified by the constraint checking module before the action is executed.

[0250] After the reinforcement learning model is trained and deployed to the actual running environment, the system enters the real-time policy generation and execution phase. The first key link of this phase is to obtain real-time state data by using the trained reinforcement learning aggregation policy optimization model and generate a preliminary control policy. The system collects the latest state information from each monitoring point at a period of one second, which covers the full range of operation state of the energy complex. Real-time data acquisition is realized through high-speed Ethernet, ensuring low delay and high reliability of data transmission.

[0251] The preprocessing of state data is an important step to ensure the quality of model input. The original monitoring data may contain measurement noise, communication interference or abnormal values caused by device failure, which will seriously affect the decision quality if not handled in time. The system adopts a multi-level data cleaning and verification mechanism, first filters out obvious abnormal values through rationality check, and then detects and corrects potential data errors through time series analysis. For short-time data missing, the system uses interpolation method for completion; for longer time data interruption, the system will enable the backup sensor or use historical data for estimation.

[0252] The preprocessed state data is input into the Actor network for policy inference. The Actor network quickly calculates the optimal participation ratio and aggregation mode selection of each type of resource based on the current state. This calculation process usually takes tens of milliseconds to complete, meeting the delay requirements of real-time control. The policy output by the network includes a continuous resource allocation vector and a discrete mode selection result, forming a complete control instruction framework.

[0253] Constraint checking of the policy is a key mechanism to ensure the safe and reliable operation of the system. The constraint checking module needs to verify whether the generated policy meets various physical constraints and operational constraints. Physical constraints include the upper and lower limits of the power of various devices, the ramp rate limit, the minimum operating time, and the minimum downtime, and other technical boundary conditions. Operational constraints include power balance requirements of the system, line transmission capacity limits, node voltage range constraints, and system stability requirements. The constraint checking uses a fast numerical calculation method, which can complete the verification of all constraint conditions within milliseconds.

[0254] When it is detected that the policy violates certain constraint conditions, the system automatically starts the projection algorithm to modify the infeasible action vector to the feasible region. The projection algorithm uses a quadratic programming method to make the minimum adjustment to the part that violates the constraints while maintaining the basic intention of the policy. This adjustment strategy ensures that the modified control instructions not only meet all the constraint conditions, but also as close as possible to the original optimization policy. The projection process also considers the priority of the constraints, giving the highest priority to safety-related hard constraints, and allowing moderate relaxation for economic-related soft constraints.

[0255] After the constraint checking and projection correction are completed, the system generates the final aggregated policy. This policy not only contains the target output values of each resource, but also contains the time window, priority identification, and backup solution of the execution. The representation of the policy uses a standardized data format, which is convenient for subsequent transmission and execution. At the same time, the system also generates the confidence assessment of the policy, which assesses the reliability of the policy execution effect based on the uncertainty of the current state and the prediction accuracy of the model.

[0256] S5.2: Based on the optimal aggregated policy, send adjustment instructions to the flexible resources through a communication protocol and collect actual execution effect data to obtain actual response data; In step S5.2, the adjusted policy is sent to each resource controller through the Modbus protocol to achieve distributed coordinated control. The control instruction contains two parts of information: target output power and response priority. After receiving the instruction, the resource controller confirms the execution feasibility according to its own state and capability boundary, and returns the feedback information to the central control system. During the execution process, the system continuously monitors the actual output and response characteristics of each resource, records the response delay, adjustment deviation, and energy conversion efficiency, and other key indicators.

[0257] After the strategy generation, the system enters the instruction transmission and execution monitoring phase. Based on the generated optimal aggregation strategy, the system needs to convey specific control instructions to the controllers of each flexibility resource. Instruction transmission adopts a variety of standardized industrial communication protocols, including Modbus, IEC61850, DNP3, etc., ensuring reliable communication with devices of different manufacturers and different technical routes. The selection of communication protocols is optimized according to the characteristics of device types and application scenarios. For energy storage devices that require high real-time performance, high-speed Ethernet protocols are used, while for thermal devices with relatively low real-time requirements, lower-cost serial communication protocols can be used.

[0258] The design of control instructions adopts a structured information format, each instruction contains multiple key components. The target output value is the core content of the instruction, which clearly specifies the power output level or other control targets that the device needs to achieve. The response time requirement specifies how long the device must complete the adjustment action, which is dynamically adjusted according to the current system demand and device capability. The priority identifier helps the device controller determine the execution order when facing multiple instructions, ensuring that the most important adjustment action can be executed first. The backup instruction provides an alternative when the main instruction cannot be executed, enhancing the fault tolerance of the system.

[0259] After receiving the adjustment instruction, the resource controller first performs a local feasibility check. This check process considers factors such as the current state, health status, maintenance status of the device, and external environmental conditions. If the controller determines that the instruction can be safely executed, it will send confirmation information to the central system and begin executing the adjustment action. If the controller finds that the instruction has safety risks or is technically unfeasible, it will feedback specific problem information to the central system and provide possible alternative solutions. This distributed safety check mechanism establishes an effective balance between central decision-making and local safety.

[0260] During the execution of the instruction, the system continuously monitors the actual response of each resource. The monitoring content includes multiple dimensions such as the actual output power of the device, response delay, adjustment accuracy, energy conversion efficiency, and device state. The monitoring of actual output power helps the system understand whether the adjustment effect meets expectations, the monitoring of response delay helps evaluate the dynamic performance of the device, the monitoring of adjustment accuracy reflects the accuracy of the control system, and the monitoring of conversion efficiency relates to the economy of the system.

[0261] The system also establishes a detection and handling mechanism for abnormal responses. When the actual response of a certain resource deviates significantly from the expected response, the system automatically starts the abnormal handling program. Minor deviations may only trigger minor adjustments to the adjustment parameters, while serious deviations may result in the resource being temporarily excluded from the scheduling range and other resources being activated for compensation. The abnormality detection uses a statistical learning-based method to distinguish between normal performance fluctuations and real abnormal conditions, avoiding false positives and false negatives.

[0262] To optimize communication efficiency and reduce network burden, the system adopts an intelligent data transmission strategy. For frequently changing key parameters, high-frequency transmission is adopted, and for relatively stable state information, lower transmission frequency is adopted. At the same time, the system also realizes data compression and incremental transmission technology, only transmitting the changed data part, greatly reducing the occupation of communication bandwidth. In the case of poor network conditions, the system can automatically adjust the data transmission strategy to ensure that the most critical information can be transmitted in a timely and reliable manner.

[0263] S5.3: Calculate the actual reward value for the actual response data and store it in the experience replay buffer, start online learning update to realize the adaptive dynamic optimization of the system.

[0264] In step S5.3, the online learning update mechanism ensures that the model continuously adapts to system changes. After each control cycle, actual execution effect data is collected, including the actual output power, response time, adjustment accuracy, etc. of each resource, the actual reward value is calculated and stored in the experience replay buffer. When the buffer accumulates enough new experience, online learning update is started. Online learning uses a small learning rate (10 -5 ) and batch size (32) to avoid excessive disturbance to trained parameters. To prevent catastrophic forgetting, the Elastic Weight Consolidation (EWC) technique is used to add a regularization term to the loss function, where is the importance weight, is the key parameter value. Model parameters are updated every hour to ensure that historical knowledge is maintained while adapting to new operating modes.

[0265] Based on instruction execution and response monitoring, the system enters the continuous learning and adaptive optimization phase. The core goal of this phase is to enable the reinforcement learning model to continuously improve decision-making performance based on actual operating experience, adapt to changes in system characteristics and new operating modes. First, the system needs to collect and organize actual response data, which includes not only the output performance of each resource, but also the overall system operating effect, user satisfaction, economic benefit, and other comprehensive indicators.

[0266] The calculation of actual reward values is a key step in online learning. Unlike the theoretical reward function used in the training phase, actual reward values are calculated based on real running results, more accurately reflecting the actual effect of policy execution. Response speed reward is calculated based on measured response time data, adjustment accuracy reward is calculated based on the deviation of actual adjustment effect from the target value, running cost reward is calculated based on actual fees incurred, system stability reward is calculated based on power quality monitoring data, and environmental benefit reward is calculated based on actual energy consumption and emission data. This reward calculation based on actual running results makes the learning process more realistic and improves the practicality of the model.

[0267] The management of the experience replay buffer adopts advanced storage and retrieval strategies. Due to the huge amount of data generated by online operation, the system needs to intelligently select and retain the most valuable experience samples. The priority sampling mechanism assigns storage priorities according to the learning value of experience, and those containing rare conditions, important decisions or significant learning effects are given higher retention priorities. At the same time, the system also implements the timeliness management of experience, and the too old experience will be gradually eliminated to make room for new experience. This dynamic experience management ensures that the buffer always contains the most relevant and valuable learning samples.

[0268] The parameter update strategy of online learning needs to balance between learning new knowledge and preserving historical knowledge. The learning rate is set much lower than the initial training phase, usually one-tenth to one percent of the training phase, which avoids the impact of new experience on learned knowledge. The batch size is also adjusted to a smaller value, usually between 16 and 64, ensuring that each update is a gradual improvement rather than a dramatic change. The update frequency is set to once an hour, maintaining learning sensitivity while avoiding excessive frequent parameter adjustment.

[0269] To prevent the problem of catastrophic forgetting, the system implements the flexible weight merging technique. This technique calculates the importance of each model parameter to the performance of historical tasks, and applies stronger constraints to important parameters when learning new tasks. In specific implementation, the system maintains an importance weight matrix, which records the sensitivity of each parameter to key performance indicators. During parameter update, the change range of important parameters is limited to a small range, while less important parameters are allowed to adjust more. This selective constraint mechanism effectively protects key knowledge from being destroyed by new learning processes.

[0270] Adaptive optimization is also reflected in the dynamic adjustment of model structure. When the system detects new operating modes or resource types, it may need to expand the expressive power of the model. The system implements a gradual network structure expansion mechanism, which can add new neurons or connections without affecting existing functions. The newly added network components quickly adapt to new task requirements through transfer learning, while the original network structure remains relatively stable. This gradual expansion ensures that the system can continuously adapt to changing application environments.

[0271] The system also establishes performance monitoring and evaluation mechanisms to regularly assess the effectiveness of online learning and the trend of model performance changes. Key performance indicators include decision accuracy, response speed, economic benefits, system stability, and other dimensions. By comparing the performance changes before and after online learning, the system can objectively evaluate the learning effect and adjust the learning strategy or roll back to the previous model version if necessary. This continuous performance monitoring ensures that the system always remains in the best operating state.

[0272] As shown in Figure 2 Corresponding to the above method, the application also provides a dynamic aggregation system for energy complex based on reinforcement learning, comprising an acquisition module, an identification module, an analysis module, a construction module and an output module.

[0273] The acquisition module is used to acquire the operation data of flexible resources in the energy complex, collect real-time operation parameters of electric, gas, heat and multi-energy coupling devices through an industrial Ethernet protocol, and use neural differential equation process modeling technology to construct a continuous-time dynamic system model; The identification module is used to preliminarily identify the energy aggregation problem based on the continuous-time dynamic system model, activate an adaptive dynamic pulse neural network classifier after identifying potential aggregation needs, and accurately classify specific energy aggregation problems using a static pulse neural network; The analysis module is used to determine corresponding data processing strategies and time window parameters using different aggregation demand types in the accurate classification results of the specific energy aggregation problem, align and analyze energy data measured at non-equidistant time points using a multi-marginal random flow matching method, and quantitatively evaluate the response characteristics of each flexible resource through measurement value spline enhancement technology and fractional matching, wherein the response characteristics of each flexible resource include the randomness of state transition and the sequence of decision-making in the energy aggregation process; The construction module is used to model the energy aggregation process as a Markov decision process, train and optimize using a twin-delayed deep deterministic policy gradient algorithm, and generate a reinforcement learning aggregation strategy optimization model; The output module is configured to receive state data in real time by using the reinforcement learning aggregated policy optimization model, output an optimal aggregated policy, send adjustment instructions to the flexible resources by using the optimal aggregated policy through a distributed control system, collect actual response data for online learning and updating, and realize self-adaptive dynamic optimization.

[0274] The application combines advanced methods such as neural ordinary differential equation modeling technology, double-layer pulse neural network classification, multi-marginal random flow matching and reinforcement learning policy optimization, and constructs a complete energy complex flexible dynamic aggregation method and system. The method can accurately model the dynamic characteristics of a complex multi-element heterogeneous energy system, effectively process high-dimensional, non-equidistant time series data, realize dynamic aggregation and optimal scheduling of flexible resources in an energy complex, and has a broad application prospect.

[0275] The above description is only preferred embodiments of the present application and is not intended to limit the present application. The present application can have various changes and modifications for those skilled in the art. Any modification, equivalent replacement, improvement, etc. within the spirit and principles of the present application shall be included in the protection scope of the present application.

[0276] The present disclosure also provides a computer readable storage medium having a computer program stored thereon, wherein the computer program is run by a processor to perform the steps of the energy complex flexible dynamic aggregation system method based on reinforcement learning provided in the above method embodiments. The storage medium can be a volatile or non-volatile computer readable storage medium.

[0277] In addition, the present disclosure also provides a computer program product having a computer program stored thereon, wherein the computer program is run by a processor to perform the steps of the energy complex flexible dynamic aggregation system method based on reinforcement learning provided in any of the above embodiments. For details, please refer to the above method embodiments, which will not be repeated here.

[0278] The computer program product can be specifically implemented by hardware, software or a combination thereof. In an optional embodiment, the computer program product is specifically embodied as a computer storage medium, which can be a volatile or non-volatile computer readable storage medium. In another optional embodiment, the computer program product is specifically embodied as a software product, such as a software development kit (Software Development Kit, SDK) and the like.

[0279] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the device and the apparatus described above can refer to the corresponding process in the foregoing method embodiment, and will not be repeated here. In several embodiments provided in the present disclosure, it should be understood that the disclosed device, apparatus and method can be implemented in other ways. The apparatus embodiments described above are only schematic, for example, the division of the units is only a logical function division, and there can be another division in actual implementation, for example, a plurality of units or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the units shown or discussed can be indirect coupling or communication connection through some communication interfaces, devices or units, and can be electrical, mechanical or other forms.

[0280] The units described as separate components can or can not be physically separate, and the components shown as units can or can not be physical units, i.e., they can be located in one place or distributed on multiple network units. Some or all of the units can be selected according to actual needs to achieve the purpose of the embodiment.

[0281] In addition, the functional units in each embodiment of the present disclosure can be integrated in one processing unit, or each unit can be physically present separately, or two or more units can be integrated in one unit.

[0282] If the functions are realized in the form of software function units and sold or used as independent products, they can be stored in a non-volatile computer readable storage medium executable by a processor. Based on this understanding, the technical solutions of the present disclosure or the part of the present disclosure that essentially contributes to the prior art or the part of the technical solutions can be embodied in the form of a software product, which is stored in a storage medium and includes a plurality of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present disclosure. The foregoing storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and various program code storage media.

[0283] Finally, it should be noted that the above-described embodiments are merely specific embodiments of the present disclosure, used to illustrate the technical solutions of the present disclosure, and are not intended to limit the present disclosure. The protection scope of the present disclosure is not limited thereto. Although the present disclosure has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can make modifications or easy changes to the technical solutions described in the foregoing embodiments, or easily think of changes or equivalent replacements for some of the technical features; and these modifications, changes or replacements do not cause the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present disclosure, and should be covered within the protection scope of the present disclosure. Therefore, the protection scope of the present disclosure should be subject to the protection scope of the claims.

Claims

1. A reinforcement learning based energy complex dynamic optimization method, characterized in that, The method comprises the following steps: acquiring operation data of flexible resources in an energy complex, collecting real-time operation parameters of electric, gas, heat and multi-energy coupling equipment, and using neural differential equation process modeling technology to build a continuous time dynamic system model; based on the continuous time dynamic system model, using a static pulse neural network to preliminarily identify an energy aggregation problem, activating an adaptive dynamic pulse neural network classifier after identifying potential aggregation demand, and accurately classifying a specific energy aggregation problem; acquiring different aggregation demand types in the accurate classification result of the specific energy aggregation problem, using a multi-marginal random flow matching method to align and analyze energy data measured at non-equidistant time points, quantitatively evaluating the response characteristics of each flexible resource through measurement value spline enhancement technology and fractional matching, wherein the response characteristics of each flexible resource include the randomness of state transition and the sequence of decision-making in the energy aggregation process; modeling the energy aggregation process as a Markov decision process, and using a twin-delayed deep deterministic policy gradient algorithm for training and optimization to generate a reinforcement learning aggregation strategy optimization model; using the reinforcement learning aggregation strategy optimization model to receive state data in real time and output an optimal aggregation strategy, using the optimal aggregation strategy to send adjustment instructions to each flexible resource through a distributed control system, and collecting actual response data for online learning and updating to realize adaptive dynamic optimization.

2. The method of claim 1, wherein, acquiring operation data of flexible resources in an energy complex, collecting real-time operation parameters of electric, gas, heat and multi-energy coupling equipment, and using neural differential equation process modeling technology to build a continuous time dynamic system model, comprising: acquiring the operation data of each flexible resource, using a three-sigma rule to remove data points outside the normal range and eliminating dimension differences through Z-score standardization to generate preprocessed operation data; based on the preprocessed operation data, representing the energy system state evolution process as a differential equation, using a neural network with a residual network architecture to approximate the differential equation function and calculating the gradient through the adjoint sensitivity method to establish a battery energy storage model, a gas turbine model and a thermal energy storage model; using a variable step ordinary differential equation solver to train the battery energy storage model, the gas turbine model and the thermal energy storage model respectively to obtain the continuous time dynamic system model.

3. The method of claim 1, wherein, based on the continuous time dynamic system model, using a static pulse neural network to preliminarily identify an energy aggregation problem, activating an adaptive dynamic pulse neural network classifier after identifying potential aggregation demand, and accurately classifying a specific energy aggregation problem, comprising: acquiring normalized state data, using a leaky integrate-and-fire neuron model to perform pulse coding through the static pulse neural network, identifying four types of aggregation demand, including load peak-valley adjustment, frequency adjustment, voltage support and emergency response; based on the four types of aggregation demand, activating an adaptive dynamic pulse neural network classifier, and automatically increasing neurons to expand network capacity according to the novelty of input data through the adaptive dynamic pulse neural network classifier to dynamically adjust the network structure; The connection weights of the dynamic spiking neural network classifier are updated by an adaptive spike-timing-dependent plasticity rule, and a specific energy aggregation problem is accurately classified by the dynamic spiking neural network classifier.

4. The method of claim 3, wherein, Based on the four types of aggregation demands, an adaptive dynamic spiking neural network classifier is activated, and the network structure is dynamically adjusted according to the novelty of the input data, the network capacity is automatically expanded by the adaptive dynamic spiking neural network classifier, including: The feature vectors of the four types of aggregation demands are obtained, the Euclidean distance between the input data and the existing neuron prototype vector is calculated, and when the minimum distance exceeds the preset novelty threshold, it is determined that it is a new mode, and the network structure adjustment mechanism is triggered; Based on the network structure adjustment mechanism, new hidden layer neurons are created in the corresponding output category, the feature vector of the current input data is taken as the prototype vector of the new neuron, and the connection weights of the new neuron with the input layer and the output layer are initialized; The activation strength of the new neuron in the network is updated through the competitive learning mechanism, the optimal response neuron is determined, and the expanded network capacity is finally determined.

5. The method of claim 1, wherein, The multi-marginal stochastic flow matching method is used to align and analyze the energy data measured at non-equidistant time points, the measured value spline enhancement technology and fractional matching are used to quantitatively evaluate the response characteristics of each flexible resource, including: The data distribution characteristics of different types of energy equipment are obtained, and multiple marginal distributions of electric energy storage resources, thermal energy storage resources and gas energy storage resources are constructed, and non-parametric modeling is performed through kernel density estimation method to generate a smoothed distribution model; Based on the smoothed distribution model, a cubic B-spline basis function is used to process the observation values at irregular time points through the measured value spline technology to generate a continuous time series; The fractional function difference between the minimum real data distribution of the continuous time series and the distribution model is calculated by using the fractional matching algorithm to obtain the response characteristics of each flexible resource.

6. The method of claim 5, wherein, The fractional function difference between the minimum real data distribution of the continuous time series and the distribution model is calculated by using the fractional matching algorithm to obtain the response characteristics of each flexible resource, including: The real data distribution of the continuous time series is obtained, the logarithmic density gradient of the real data distribution is calculated and taken as the real fractional function, and the logarithmic density gradient of the distribution model is calculated and taken as the fractional function of the distribution model; Based on the real fractional function and the fractional function of the distribution model, a fractional matching objective function is constructed, the least squares method is used to calculate the square error between the real fractional function and the fractional function of the distribution model, and the calculation result is integrated and summed; The model parameters of the fractional matching objective function are optimized by the gradient descent algorithm to make the fractional function of the distribution model approach the real fractional function, and the response characteristics of each flexible resource including the response time constant, the regulation accuracy and the stability index are obtained.

7. The method of claim 5, wherein, Based on the smoothed distribution model, a measurement value spline technique is used to construct a cubic B-spline basis function to process the observation values at irregular time points, and a continuous time sequence is generated, including: Obtain the observation values at irregular time points in the smoothed distribution model, determine the spline node sequence, set the node density according to the time interval distribution of the observation data, and increase the number of nodes in the data-intensive area; Based on the node sequence, the cubic B-spline basis function is constructed, and each cubic B-spline basis function is a cubic polynomial within four consecutive node intervals and zero in other intervals; The coefficients of the cubic B-spline basis function are determined by the least squares fitting method, so that the spline curve passes through or approximates all observation points, and the continuous time sequence with time continuity and smoothness is obtained.

8. The method of claim 1, wherein, The energy aggregation process is modeled as a Markov decision process, and a twin delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregation policy optimization model, including: Model the energy aggregation process as a Markov decision process, construct a state space and an action space, and design a reward function; Based on the Markov decision process framework, state space, action space and reward function, a twin delayed deep deterministic policy gradient algorithm is used for training and optimization to generate the reinforcement learning aggregation policy optimization model.

9. The method of claim 8, wherein, Based on the Markov decision process framework, state space, action space and reward function, a twin delayed deep deterministic policy gradient algorithm is used for training and optimization to generate a reinforcement learning aggregation policy optimization model, including: Obtain multi-dimensional characteristic data, construct a state space containing load demand dimension, renewable energy dimension, energy storage state dimension, power quality dimension and economic signal dimension, and an action space containing resource participation degree allocation vector and aggregation mode selection, and establish a Markov decision process framework; In the Markov decision process framework, a multi-objective reward function is designed, including response speed reward item, adjustment accuracy reward item, operation cost reward item, system stability reward item and environmental benefit reward item, and a comprehensive performance evaluation index is generated; Based on the comprehensive performance evaluation index, an improved twin delayed deep deterministic policy gradient algorithm is used to construct an actor network and a critic network to obtain the reinforcement learning aggregation policy optimization model.

10. The method of claim 1, wherein, The reinforcement learning aggregation policy optimization model is used to receive real-time state data and output an optimal aggregation policy, and the optimal aggregation policy is used to send adjustment instructions to each flexible resource through a distributed control system, and actual response data is collected for online learning and updating to realize adaptive dynamic optimization, including: The reinforcement learning aggregation policy optimization model is used to obtain real-time state data, and a constraint checking module is used to verify the feasibility of the reinforcement learning aggregation policy. If the reinforcement learning aggregation policy violates the constraints, a projection algorithm is used to project the action vector into the feasible region to generate the optimal aggregation policy; Based on the optimal aggregation policy, adjustment instructions are sent to each flexible resource through a communication protocol, and actual execution effect data is collected to obtain actual response data; Actual reward values of the actual response data are calculated and stored to an experience replay buffer, online learning update is started, and adaptive dynamic optimization is realized.

11. A reinforcement learning based energy complex dynamic optimization system, characterized in that, The method comprises the following steps: An acquisition module is configured to acquire operation data of flexible resources in an energy complex, collect real-time operation parameters of electric, gas, heat and multi-energy coupling devices through an industrial Ethernet protocol, and construct a continuous-time dynamic system model by using a neural differential equation process modeling technology; An identification module is configured to preliminarily identify an energy aggregation problem by using a static pulse neural network based on the continuous-time dynamic system model, activate an adaptive dynamic pulse neural network classifier after identifying a potential aggregation demand, and accurately classify a specific energy aggregation problem; An analysis module is configured to determine corresponding data processing strategies and time window parameters according to different aggregation demand types in the accurate classification result of the specific energy aggregation problem, align and analyze energy data measured at non-equidistant time points by using a multi-marginal random flow matching method, and quantitatively evaluate response characteristics of each flexible resource by using a measured value spline enhancement technology and fractional matching, wherein the response characteristics of each flexible resource include randomness of state transition and sequence of decision-making in the energy aggregation process; A construction module is configured to model the energy aggregation process as a Markov decision process, train and optimize the energy aggregation process by using a twin-delayed deep deterministic policy gradient algorithm, and generate a reinforcement learning aggregation strategy optimization model; An output module is configured to receive state data in real time and output an optimal aggregation strategy by using the reinforcement learning aggregation strategy optimization model, send adjustment instructions to each flexible resource through a distributed control system by using the optimal aggregation strategy, collect actual response data for online learning update, and realize adaptive dynamic optimization.

Images