Virtual power plant multi-energy collaborative scheduling and real-time optimization control method and system

Abstract

The application provides a kind of virtual power plant multi-energy collaborative scheduling and real-time optimization control method and system, it is related to energy management technical field, including: obtaining energy state data identifies coupling relationship, constructs space-time correlation diagram and calculates embedding representation, identifies strong space-time correlation node pair based on embedding representation and carries out sub-domain division, configures boundary interaction protocol, collects response data and calculates deviation, predicts future influence trajectory, when trajectory exceeds allowable fluctuation range, boundary adjustment is executed and modified instruction is obtained by rolling optimization.It can improve the accuracy of multi-energy collaborative scheduling and real-time control effect.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of energy management, and in particular to a virtual power plant multi-energy collaborative scheduling and real-time optimization control method and system. BACKGROUND

[0002] With the diversification of energy structure and large-scale access of renewable energy, centralized energy management system is facing great challenges. As a new energy management mode, virtual power plant can realize the collaborative optimization and scheduling of multi-energy system by integrating distributed energy resources. The core is to use information communication technology to aggregate dispersed energy resources into a schedulable whole, participate in power market transactions and system regulation, thereby improving flexibility and reliability.

[0003] With the development of energy internet and smart grid, the existing virtual power plant technology mainly uses centralized or hierarchical distributed methods for energy scheduling optimization. By establishing a mathematical model and considering the constraint conditions of various energy resources, the optimization problem is solved to obtain the optimal scheduling scheme. However, the existing technology still lacks accurate identification and characterization of complex coupling relationships between energy nodes, which affects the rationality and feasibility of the scheduling scheme. The static partition of the scheduling domain structure cannot adaptively adjust the scheduling strategy according to the dynamic changes of the spatio-temporal correlation characteristics between energy nodes, making it difficult to respond to real-time changes in system operating conditions and external environments. In addition, there is a lack of effective deviation identification, evolution prediction and correction mechanisms, which cannot timely adjust the scheduling strategy to cope with system disturbances and uncertainties. SUMMARY

[0004] The present application provides a virtual power plant multi-energy collaborative scheduling and real-time optimization control method and system, which can at least solve some of the problems existing in the prior art.

[0005] In a first aspect of the present application, a virtual power plant multi-energy collaborative scheduling and real-time optimization control method is provided, comprising:

[0006] Obtaining state data of energy resources and identifying coupling relationships between energy nodes to obtain a coupling relationship set;

[0007] Constructing a spatio-temporal correlation graph according to the coupling relationship set and calculating the embedding representation of each energy node, identifying strong spatio-temporal correlation node pairs based on the embedding representation, and dividing energy nodes containing strong spatio-temporal correlation node pairs into the same scheduling sub-domain through spectral clustering, configuring a boundary interaction protocol for each scheduling sub-domain, identifying the causal influence strength of interaction variables in the boundary interaction protocol on the scheduling sub-domain through causal inference method, and modifying the information transmission weight in the boundary interaction protocol according to the causal influence strength to obtain multiple enhanced scheduling sub-domains;

[0008] The actual response data of the energy nodes in each enhanced scheduling sub-domain is collected and compared with the expected response data to obtain a response deviation. The response deviation is sequentially matched with historical scheduling scenario data to obtain a time sequence evolution mode of the response deviation and predict a future influence trajectory corresponding to the response deviation. When the future influence trajectory exceeds a permitted fluctuation range corresponding to the coupling relationship set, a corresponding affected scheduling sub-domain is determined and boundary coupling relationship adjustment is performed. A forward iteration is solved by combining a rolling time domain optimization method to obtain a modified scheduling instruction and perform the same.

[0009] In an optional implementation,

[0010] The state data of the energy resources is obtained and the coupling relationship between the energy nodes is identified to obtain a coupling relationship set, including:

[0011] The power output data, energy conversion data, and power exchange data between resources in the state data are collected as state data.

[0012] Based on the power exchange data between resources in the state data, a power transfer matrix is constructed and singular value decomposition is performed to obtain singular values and corresponding singular vectors. A dominant path is identified from the singular vectors corresponding to the singular values exceeding a preset singular value threshold, and adjacent energy node pairs in the dominant path are taken as candidate node pairs.

[0013] For each node pair in the candidate node pair, the power transfer amount between the source node and the target node accounts for the power ratio of the total power output of the source node based on the power output data. If it exceeds a preset ratio threshold, the current node pair is marked as a coupling node pair.

[0014] The power transfer amount and transfer delay corresponding to each coupling node pair in the power transfer matrix are extracted, and the conversion efficiency boundary is extracted from the energy conversion data. The power transfer amount, conversion efficiency boundary, and transfer delay are taken as coupling relationship parameters and encapsulated with node pair identification as a coupling relationship set.

[0015] In an optional implementation,

[0016] A space-time correlation graph is constructed according to the coupling relationship set and the embedded representation of each energy node is calculated. Based on the embedded representation, a strong space-time correlation node pair is identified, and the energy nodes containing the strong space-time correlation node pair are divided into the same scheduling sub-domain by spectral clustering, including:

[0017] The coupling node pairs are extracted from the coupling relationship set as edge topology, the energy nodes are taken as vertex topology, and the pre-obtained coupling relationship parameters are taken as edge attributes to construct a space-time correlation graph. The power change gradient and phase characteristics corresponding to the state data are calculated.

[0018] For each vertex in the spatiotemporal correlation graph, an initial representation is constructed based on the power change gradient, phase feature and edge attribute, a Hadamard product of the initial representations of adjacent energy nodes and power transfer amount is performed to obtain a power representation, an element-wise reciprocal weighting is performed on the power representation and the pre-acquired transfer delay, a nonlinear mapping is applied to obtain a projected representation, a cosine distance between the projected representation and the initial representation is calculated for normalized aggregation and skip connection to obtain an embedded representation;

[0019] A bilinear measure between embedded representations is calculated, a bidirectional power flow proportion of the coupled node pair within a preset time window is extracted, when the bidirectional power flow proportion exceeds a preset power flow threshold and the bilinear measure exceeds a preset similarity threshold, the coupled node pair is marked as a strong spatiotemporal correlation node pair;

[0020] An affinity matrix is constructed based on the embedded representation and a diffusion kernel processing is applied to obtain a diffusion matrix, spectral decomposition is performed based on the diffusion matrix to obtain a spectral embedding, and spectral clustering is performed on the strong spatiotemporal correlation node pair, and the energy nodes containing the strong spatiotemporal correlation node pair are divided into the same scheduling sub-domain according to the clustering result.

[0021] In an optional implementation,

[0022] A boundary interaction protocol is configured for each scheduling sub-domain, and a causal inference method is used to identify the causal influence strength of the interaction variables in the boundary interaction protocol on the scheduling sub-domains, including:

[0023] For each scheduling sub-domain, the energy node coupling relationship parameters between the current scheduling sub-domain and the adjacent scheduling sub-domain are taken as boundary interaction variables and a boundary interaction protocol is configured, a multivariate time series of the boundary interaction variables and the state data in the scheduling sub-domain is established, a vector autoregression fitting is performed to extract a residual sequence, a time-varying precision matrix is calculated based on the residual sequence, an adaptive lasso regularization constraint is applied to the time-varying precision matrix, and a sparse precision matrix is obtained by solving through coordinate descent iteration, and a variable pair corresponding to a non-zero element in the sparse precision matrix is taken as a direct causal edge to construct a time-varying causal graph;

[0024] Based on the direct causal edges in the time-varying causal graph, an exogenous processing is performed on the historical observation data corresponding to the energy nodes to obtain a control sample set, a propensity score of the boundary interaction variables on the state data in the scheduling sub-domain is calculated, a weighted observation data is obtained by weighting the observation data of the control sample set based on the reciprocal of the propensity score, a counterfactual intervention is applied to the boundary interaction variables, and a time difference component is obtained by calculating the difference of the power change gradient corresponding to the state data in the scheduling sub-domain before and after the intervention based on the weighted observation data, a group difference component is obtained by twice differencing the control sample set based on the time difference component, and the group difference component is output as the causal influence strength.

[0025] In an alternative embodiment,

[0026] According to the causal influence intensity, the information transmission weight in the boundary interaction protocol is modified to obtain a plurality of enhanced scheduling sub-domains, including:

[0027] For each scheduling sub-domain, a sensitivity coefficient of the boundary interaction variable to the state data in the scheduling sub-domain is calculated based on the causal influence intensity, an adjustment amplitude is obtained by nonlinearly mapping the information transmission weight of each boundary interaction variable in the boundary interaction protocol through the sensitivity coefficient, the adjustment amplitude is weighted and fused with the information transmission weight to obtain a modified information transmission weight;

[0028] Based on the modified information transmission weight, the information transmission time delay between adjacent scheduling sub-domains is calculated, the inter-domain interaction features of the boundary interaction variables in the boundary interaction protocol are aggregated according to the modified information transmission weight, the inter-domain interaction features are time-aligned through the information transmission time delay and spliced with the state data in the scheduling sub-domain to obtain a plurality of enhanced scheduling sub-domains.

[0029] In an alternative embodiment,

[0030] The actual response data and the expected response data of the energy nodes in each enhanced scheduling sub-domain are collected and compared to obtain a response deviation, and the response deviation is sequence-matched with the historical scheduling scene data to obtain a time evolution pattern of the response deviation and predict a future influence trajectory corresponding to the response deviation, including:

[0031] The actual power output of the energy nodes in each enhanced scheduling sub-domain is collected as actual response data, and a preset target power output is taken as expected response data, the response deviation between the actual response data and the expected response data is calculated and singular spectrum analysis is performed to obtain a reconstruction component and a residual component, the Lyapunov index of the reconstruction component is calculated and spliced with the reconstruction component and the residual component to obtain a response deviation sequence;

[0032] extracting boundary interaction variables from pre-acquired historical scheduling scenario data to construct a probability distribution, calculating a maximum mean discrepancy between boundary interaction variables of a current scheduling scenario and the probability distribution, extracting historical scheduling scenarios with a maximum mean discrepancy less than a preset discrepancy threshold to construct a historical scenario subset, extracting a common latent space representation of the response bias sequence and historical response bias sequences in the historical scenario subset through kernel canonical correlation analysis, calculating Mahalanobis distance based on the common latent space representation to determine sequence similarity and taking historical response bias sequences with sequence similarity exceeding a preset similarity threshold as matching sequences, performing Bayesian structural time series decomposition on the matching sequences to extract trend evolution paths and periodic disturbance components and determine a time series evolution mode, applying the time series evolution mode to the response bias sequence and performing multi-step forward deduction through integrated Kalman filtering to obtain the future impact trajectory.

[0033] In an optional implementation,

[0034] When the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, a corresponding affected scheduling sub-domain is determined and boundary coupling relationship adjustment is performed, forward iteration is solved in combination with a rolling time domain optimization method, modified scheduling instructions are obtained and execution includes:

[0035] determining whether the deviation amplitude of each prediction time in the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, if it exceeds, extracting the corresponding exceeding time as a prediction time, obtaining the energy node identifier corresponding to the prediction time and determining the affected scheduling sub-domain based on the energy node identifier, extracting the boundary coupling relationship between the affected scheduling sub-domain and the adjacent scheduling sub-domain, constructing the Jacobian matrix corresponding to the power exchange constraint in the boundary coupling relationship and calculating the directional derivative to obtain the sensitivity coefficient, modifying the power exchange constraint based on the sensitivity coefficient to obtain an optimized boundary coupling relationship;

[0036] taking the optimized boundary coupling relationship as a consensus constraint, setting the prediction time domain and the control time domain of the rolling time domain optimization, constructing an objective function in the prediction time domain and solving to obtain a Pareto optimal control sequence in the control time domain through an alternating direction multiplier method, extracting the control quantity of the first period in the Pareto optimal control sequence as a current execution action, converting the current execution action into a scheduling instruction and issuing it to the energy nodes in the affected scheduling sub-domain, collecting actual execution data of the energy nodes to update the actual response data and reconstructing the prediction time domain for forward iteration to obtain the modified scheduling instructions and execute them.

[0037] In a second aspect of the embodiments of the application, a virtual power plant multi-energy collaborative scheduling and real-time optimization control system is provided, comprising:

[0038] A first unit is configured to acquire state data of energy resources and identify coupling relationships between energy nodes to obtain a coupling relationship set;

[0039] A second unit is configured to construct a space-time correlation graph according to the coupling relationship set and calculate an embedding representation of each energy node, identify a strong space-time correlation node pair based on the embedding representation, divide energy nodes containing the strong space-time correlation node pair into a same scheduling sub-domain through spectral clustering, configure a boundary interaction protocol for each scheduling sub-domain, identify a causal influence strength of an interaction variable on a scheduling sub-domain in the boundary interaction protocol through a causal inference method, and correct an information transmission weight in the boundary interaction protocol according to the causal influence strength to obtain a plurality of enhanced scheduling sub-domains.

[0040] A third unit is configured to compare actual response data and expected response data of energy nodes in each enhanced scheduling sub-domain to obtain a response deviation, perform sequence matching on the response deviation and historical scheduling scenario data to obtain a time sequence evolution mode of the response deviation and predict a future influence trajectory corresponding to the response deviation, determine an affected scheduling sub-domain when the future influence trajectory exceeds an allowable fluctuation range corresponding to the coupling relationship set, and perform boundary coupling relationship adjustment, and perform forward iteration solving through a rolling time domain optimization method to obtain a corrected scheduling instruction and execute the corrected scheduling instruction.

[0041] In a third aspect of the embodiment of the present application, an electronic device is provided, comprising:

[0042] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0043] In a fourth aspect of the embodiment of the present application, a computer readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0044] In the present application, according to the coupling relationship between energy nodes, a space-time correlation graph is constructed and a strong space-time correlation node pair is identified, intelligent division of scheduling sub-domains is realized, the overall coordination complexity of a multi-energy system is effectively reduced, the causal influence strength of an interaction variable between scheduling sub-domains is accurately identified through a boundary interaction protocol and a causal inference method, and the information transmission weight is corrected, the collaborative scheduling efficiency of the multi-energy system is improved, the stability in the process of energy conversion and transmission is ensured, the time sequence evolution mode analysis and future influence trajectory prediction of the response deviation are performed, boundary coupling relationship adjustment and rolling time domain optimization can be performed in time, real-time optimization control of the multi-energy system is realized, and the reliability and economy of the virtual power plant operation are effectively improved. BRIEF DESCRIPTION OF DRAWINGS

[0045] Figure 1A flowchart of the virtual power plant multi-energy collaborative scheduling and real-time optimization control method of the embodiment of the present application is shown in Fig.

[0046] Figure 2 A sensitivity optimization multi-energy coordination flowchart of the virtual power plant multi-energy collaborative scheduling and real-time optimization control method of the embodiment of the present application is shown in Fig. DETAILED DESCRIPTION

[0047] To make the purpose, technical solutions and advantages of the embodiment of the present application clearer, the technical solutions in the embodiment of the present application will be described clearly and completely below in combination with the drawings in the embodiment of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0048] The technical solutions of the present application will be described in detail below with specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes can not be described in some embodiments.

[0049] Figure 1 A flowchart of the virtual power plant multi-energy collaborative scheduling and real-time optimization control method of the embodiment of the present application is shown in Fig. Figure 1 As shown, the method comprises:

[0050] Obtain the state data of energy resources and identify the coupling relationship between energy nodes to obtain a coupling relationship set;

[0051] According to the coupling relationship set, a space-time correlation graph is constructed and the embedding representation of each energy node is calculated. Based on the embedding representation, strong space-time correlation node pairs are identified, and the energy nodes containing the strong space-time correlation node pairs are divided into the same scheduling sub-domain through spectral clustering. A boundary interaction protocol is configured for each scheduling sub-domain. The causal inference method is used to identify the causal influence strength of the interaction variables in the boundary interaction protocol on the scheduling sub-domains. According to the causal influence strength, the information transmission weight in the boundary interaction protocol is modified to obtain a plurality of enhanced scheduling sub-domains;

[0052] The actual response data of the energy nodes in each enhanced scheduling sub-domain is collected and compared with the expected response data to obtain a response deviation. The response deviation is sequentially matched with historical scheduling scenario data to obtain a time sequence evolution mode of the response deviation and predict a future influence trajectory corresponding to the response deviation. When the future influence trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, the corresponding affected scheduling sub-domain is determined and boundary coupling relationship adjustment is performed. A forward iteration is solved by combining a rolling time domain optimization method to obtain a modified scheduling instruction and execute it.

[0053] In an alternative embodiment,

[0054] The state data of the energy resources and the coupling relationship between the energy nodes are obtained to obtain a coupling relationship set, which comprises:

[0055] The power output data, energy conversion data and power exchange data between resources of each energy resource in a continuous time window are collected as state data.

[0056] A power transfer matrix is constructed based on the power exchange data between resources in the state data and singular value decomposition is performed to obtain singular values and corresponding singular vectors. The dominant path is identified from the singular vectors corresponding to the singular values exceeding the preset singular value threshold, and the adjacent energy node pairs in the dominant path are taken as candidate node pairs.

[0057] For each node pair in the candidate node pair, the power transfer amount between the source node and the target node accounts for the power ratio of the total power output of the source node based on the power output data. If it exceeds the preset ratio threshold, the current node pair is marked as a coupling node pair.

[0058] The power transfer amount and transfer delay corresponding to each coupling node pair in the power transfer matrix are extracted, and the conversion efficiency boundary is extracted from the energy conversion data. The power transfer amount, conversion efficiency boundary and transfer delay are taken as coupling relationship parameters and encapsulated with node pair identification as a coupling relationship set.

[0059] The power output data, energy conversion data and power exchange data between resources of each energy resource in a continuous time window are collected as state data, which can be realized by a data collection module of a distributed energy management system. The sampling period is usually set to 5 minutes, and the time window can be set to 24 hours. The collected state data includes the power output of each photovoltaic array, the power curve of the wind turbine generator, the charge and discharge power of the energy storage device, the electric and heat output of the combined heat and power unit, and the power flow data between each energy node. For example, for a virtual power plant containing 10 energy nodes, in a 24-hour sampling period, each node will generate 288 sampling points, forming a multi-dimensional state data matrix.

[0060] Based on the collected state data, a power transfer matrix is constructed and singular value decomposition is performed. The power transfer matrix is an n x n matrix, where n is the number of energy nodes, and the matrix element P_ij represents the power transfer amount from node i to node j. When singular value decomposition is performed, the power transfer matrix is decomposed into the product of three matrices: the left singular vector matrix, the singular value diagonal matrix, and the transpose of the right singular vector matrix. Among them, the singular values are arranged in descending order in the diagonal matrix. In actual implementation, a preset singular value threshold can be set as 80% of the sum of the total singular values, that is, the singular values and their corresponding singular vectors whose contribution accumulates more than 80% are selected. For each selected singular vector, the energy nodes corresponding to the elements with the largest absolute values are identified, and these nodes form the dominant path. For example, assuming that the element values of the photovoltaic node, the energy storage node and the load node in a singular vector are 0.85, 0.72 and 0.65 respectively, then these three nodes form a dominant path, and the adjacent node pairs (photovoltaic-energy storage, energy storage-load) are marked as candidate node pairs.

[0061] For the identified candidate node pairs, the power proportion needs to be further calculated to determine the true coupling node pairs. For each candidate node pair (i, j), the proportion of the power amount transferred from node i to node j to the total output power of node i is calculated. In specific calculation, the average power transfer amount of node i to node j in a 24-hour window is counted, and divided by the average total power output of node i. If the proportion exceeds the preset proportion threshold (usually set to 15%), the node pair is marked as a coupling node pair. For example, if the average transmission power of the photovoltaic node to the energy storage node is 300kW, and the total output power of the photovoltaic node is 1500kW, the power proportion is 20%, which exceeds the threshold of 15%, so (photovoltaic-energy storage) is determined as a coupling node pair.

[0062] For each determined coupling node pair, three types of key parameters need to be extracted to represent its coupling relationship: power transfer amount, transfer delay and conversion efficiency boundary. The power transfer amount is directly extracted from the power transfer matrix, which can be represented as the maximum value, minimum value and average value in 24 hours; the transfer delay is obtained by analyzing the time difference between the source node power change and the target node response, and the cross-correlation analysis method is usually used for calculation; the conversion efficiency boundary is extracted from the energy conversion data, including the highest efficiency and the lowest efficiency observed in the historical running process. These parameters together constitute the coupling relationship parameter set, which is packaged with the node pair identification as a coupling relationship set.

[0063] Taking a virtual power plant including photovoltaic power generation, wind power generation, energy storage system, combined heat and power and controllable load as an example, after the foregoing analysis, the following coupling relationship can be obtained: the coupling relationship parameters between photovoltaic power generation and energy storage system are {power transfer amount: [50kW, 450kW, 250kW], transfer time delay: 30 seconds, conversion efficiency boundary: [0.92, 0.98]}; the coupling relationship parameters between combined heat and power and regional heating load are {power transfer amount: [200kW, 800kW, 500kW], transfer time delay: 120 seconds, conversion efficiency boundary: [0.75, 0.85]}.

[0064] In this embodiment, by introducing singular value decomposition means to decompose the power transfer matrix, the dominant power transfer path can be extracted from the complex energy interaction network, the automatic identification of the key energy flow direction is realized, the coupling node pairs with high contribution degree are screened out through the power proportion judgment mechanism, the redundancy association is effectively reduced, the reliability and calculation efficiency of the coupling identification are improved, the power transfer amount, the conversion efficiency boundary and the transfer time delay are collectively taken as the coupling relationship parameters to encapsulate the coupling relationship set, so that the interaction between the energy resources is not only reflected at the topological level, but also quantitatively described at the energy efficiency and dynamic characteristic level, which provides more physically interpretable basis support for subsequent energy optimization scheduling, system stability analysis and coupling characteristic prediction.

[0065] In an alternative embodiment,

[0066] According to the coupling relationship set, a space-time correlation graph is constructed and an embedding representation of each energy node is calculated, and based on the embedding representation, a strong space-time correlation node pair is identified and energy nodes containing the strong space-time correlation node pair are divided into the same scheduling sub-domain through spectral clustering, comprising:

[0067] A coupling node pair is extracted from the coupling relationship set as an edge topology, an energy node is taken as a vertex topology, and a coupling relationship parameter obtained in advance is taken as an edge attribute to construct a space-time correlation graph, and a power change gradient and a phase feature corresponding to the state data are calculated;

[0068] For each vertex in the space-time correlation graph, an initial representation is constructed based on the power change gradient, the phase feature and the edge attribute, a Hadamard product of the initial representations of adjacent energy nodes and the power transfer amount is taken to obtain a power representation, an element-wise reciprocal weighting is performed on the power representation and the transfer time delay obtained in advance, a nonlinear mapping is applied to obtain a projection representation, and a cosine distance between the projection representation and the initial representation is normalized and aggregated and a skip connection is performed to obtain an embedding representation;

[0069] calculating a bilinear metric between the embedding representations, extracting bidirectional power flow proportions of the coupling node pair within a preset time window, and marking the coupling node pair as a strong spatiotemporal correlation node pair when the bidirectional power flow proportions both exceed a preset power flow threshold and the bilinear metric exceeds a preset similarity threshold;

[0070] constructing an affinity matrix based on the embedding representations and applying a diffusion kernel process to obtain a diffusion matrix, performing spectral decomposition based on the diffusion matrix to obtain spectral embedding and performing spectral clustering on the strong spatiotemporal correlation node pair, and dividing energy nodes containing the strong spatiotemporal correlation node pair into the same scheduling sub-domain according to the clustering result.

[0071] From the coupling relationship set, extract the coupling node pair as the edge topology, take the energy node as the vertex topology, and take the pre-acquired coupling relationship parameters as the edge attribute. For each energy node in the virtual power plant, such as a photovoltaic power generation node, a wind power generation node, and an energy storage node, it is taken as a vertex of a spatiotemporal correlation graph. The identified coupling node pair, such as a photovoltaic power generation-energy storage system and a combined heat and power supply-heat load, is taken as an edge of the graph. The edge attribute includes the aforementioned acquired power transfer amount, transfer time delay, and conversion efficiency boundary. In addition, the power change gradient and the phase feature corresponding to the state data also need to be calculated. The power change gradient is obtained by calculating the power difference between adjacent time points divided by the time interval, reflecting the power change rate. The phase feature is extracted by performing Hilbert transform on the power time series data, representing the periodic characteristics of power fluctuations. For example, for a photovoltaic node with a sampling interval of 5 minutes, its 24-hour power output sequence can calculate 288 gradient values and corresponding phase feature values.

[0072] For each vertex in the spatio-temporal correlation graph, the corresponding initial representation, projected representation and final embedding representation are constructed. The initial representation is constructed based on the power change gradient, phase feature and edge attribute. Specifically, the power gradient sequence of the energy node is spliced with the phase feature sequence to form an initial feature vector. For example, the initial representation of a certain photovoltaic node can be a vector of length 576, with the first 288 dimensions representing the power gradient and the last 288 dimensions representing the phase feature. For each pair of adjacent nodes, the initial representation of the source node is multiplied by the power transfer amount between the source node and the target node to obtain the power representation, which reflects the influence strength of the node feature in the power transmission process. Then, the power representation is element-wise weighted by the reciprocal of the transmission delay, i.e. each element of the power representation is divided by the corresponding delay value, so that the power transmission with shorter delay has higher weight. The weighted representation is subjected to a nonlinear mapping (such as a ReLU activation function) to obtain the projected representation. The cosine distance between the projected representation and the initial representation is calculated as the importance weight of the neighbor node information, which is normalized and aggregated. The aggregation result is concatenated with the initial representation to obtain the embedding representation. For example, for the energy storage node, its embedding representation integrates the information of the photovoltaic node and the load node connected to it, representing the role and behavior characteristics of the node in the energy network.

[0073] After obtaining the embedding representation, the strong spatio-temporal correlation node pairs are identified. The bilinear measure between the embedding representations of any two nodes is calculated, i.e. the inner product of the two vectors after transformation by a bilinear transformation matrix is calculated to obtain the similarity score between the node pairs. The bidirectional power flow proportion of the coupled node pairs within a preset time window is extracted. The bidirectional power flow proportion is defined as the percentage of the bidirectional power transmission time between the nodes in the total time. For example, for a photovoltaic-energy storage node pair, if there are 4 hours of power transmission from photovoltaic to energy storage and 2 hours of power transmission from energy storage to photovoltaic in a 24-hour window, the forward power flow proportion is 16.7% and the reverse power flow proportion is 8.3%. When the bidirectional power flow proportion of the node pair exceeds the preset power flow threshold (e.g. 5%) and the bilinear measure exceeds the preset similarity threshold (e.g. 0.75), the coupled node pair is marked as a strong spatio-temporal correlation node pair. Strong spatio-temporal correlation indicates that these nodes have a close relationship in energy flow and time response and should be cooperatively scheduled.

[0074] The nodes with strong spatiotemporal correlation are analyzed by clustering to divide the scheduling subdomains. An affinity matrix is constructed based on the embedding representation of all nodes, and the matrix elements represent the similarity between node pairs. A diffusion kernel process is applied to the affinity matrix, and the diffusion matrix is obtained by iteratively propagating similarity information. The diffusion kernel parameter is usually set to 0.5, and the iteration number is 10. Based on the diffusion matrix, spectral decomposition is performed to extract the first few eigenvectors to form the spectral embedding. The dimension of the spectral embedding can be determined according to the cumulative contribution rate, and the number of eigenvectors that explain 85% of the variance is usually selected. The spectral clustering algorithm is applied to the node pairs with strong spatiotemporal correlation, and nodes with similar spectral embedding are divided into the same category. According to the clustering results, the energy nodes containing the node pairs with strong spatiotemporal correlation are divided into the same scheduling subdomain.

[0075] In this embodiment, by constructing a spatiotemporal correlation graph with energy nodes as vertices and coupled node pairs as edges, and combining power change gradient and phase feature for multi-dimensional representation, the amplitude change, phase synchronization and delay characteristics of energy transfer can be comprehensively reflected in the feature space, realizing deep modeling of complex energy interaction mode. The multi-dimensional feature fusion is introduced through Hadamard product weighting and nonlinear mapping mechanism, so that the node embedding representation can describe both local energy transfer strength and global time sequence dependence feature, significantly improving the completeness and discriminability of feature expression

[0076] In an optional implementation,

[0077] A boundary interaction protocol is configured for each scheduling subdomain, and the causal inference method is used to identify the causal influence strength of the interaction variables in the boundary interaction protocol on the scheduling subdomains, including:

[0078] For each scheduling subdomain, the coupling relationship parameters between the energy nodes in the current scheduling subdomain and the adjacent scheduling subdomain are taken as boundary interaction variables, and a boundary interaction protocol is configured. A multivariate time series of the boundary interaction variables and the state data in the scheduling subdomain is established, and a vector autoregression fitting is performed to extract a residual sequence. A time-varying precision matrix is calculated based on the residual sequence, an adaptive lasso regularization constraint is applied to the time-varying precision matrix, and a sparse precision matrix is obtained by coordinate descent iteration. The variable pairs corresponding to the non-zero elements in the sparse precision matrix are extracted as direct causal edges to construct a time-varying causal graph;

[0079] The historical observation data corresponding to the energy node is externally treated based on the direct causal edge in the time-varying causal graph to obtain a control sample set, a tendency score of the boundary interaction variable to the state data in the scheduling sub-domain is calculated, the observation data of the control sample set is weighted based on the inverse of the tendency score to obtain weighted observation data, counterfactual intervention is applied to the boundary interaction variable, and a difference value of a power change gradient corresponding to the state data in the scheduling sub-domain before and after intervention is calculated based on the weighted observation data to obtain a time difference component, and a group difference component is obtained by twice-differentiating the control sample set based on the time difference component, and the group difference component is output as the causal influence strength.

[0080] For each scheduling sub-domain, the energy node coupling relationship parameter between the current scheduling sub-domain and the adjacent scheduling sub-domain is taken as the boundary interaction variable, and the corresponding boundary interaction protocol is configured. For example, for the renewable energy-storage collaborative sub-domain divided, the boundary interaction variables of the heat and power cogeneration sub-domain include the power transmission amount, transmission delay and conversion efficiency boundary between the two domains. The boundary interaction protocol defines the exchange frequency, exchange format and data unit of the variable, etc. In a typical configuration, the exchange frequency is set to 5 minutes, the exchange format is a predefined data structure, and the data structure contains information such as time stamp, node identification and corresponding power value.

[0081] After the boundary interaction variable is determined, a multivariate time series of the boundary interaction variable and the state data in the scheduling sub-domain is established, and vector autoregression fitting is performed. Taking the renewable energy-storage sub-domain as an example, it contains photovoltaic nodes, energy storage nodes and part of flexible load nodes, and the internal state data includes the power output of each node, the energy storage state, etc. The state data and the boundary interaction variable form a multi-dimensional vector time series, and the sequence length is usually 24 hours, that is, it contains 288 sampling points. The vector autoregression model is applied to the multivariate time series for fitting, and the lag order is set to 12, which is equivalent to considering the historical influence within 1 hour. Through model fitting, the residual sequence is extracted, that is, the difference between the actual observation value and the predicted value. For example, the power prediction value of the energy storage node at a certain time is 200kW, and the actual observation value is 220kW, and the residual is 20kW.

[0082] A time-varying precision matrix is calculated based on the extracted residual sequence. The precision matrix is the inverse of the residual sequence covariance matrix, reflecting the conditional independence relationship between variables. In the calculation process, a sliding window method is used to process non-stationary time series, and the window size is set to 60 sampling points, equivalent to 5 hours of data. For the residual data in each sliding window, the covariance matrix is calculated and inverted to obtain the precision matrix. The obtained time-varying precision matrix is subjected to adaptive lasso regularization constraint, and the regularization parameter is determined through cross-validation, usually set to between 0.01 and 0.1. After regularization, the iterative solution is performed through the coordinate descent algorithm, and the maximum iteration number is set to 1000 and the convergence threshold is set to 0.001. In the coordinate descent process, only one variable is optimized each time, and the other variables are fixed, and the iteration is repeated until convergence, obtaining a sparse precision matrix. From the sparse precision matrix, the variable pairs corresponding to the non-zero elements are extracted, which represent direct causal relationships and constitute a time-varying causal graph. For example, if the element value corresponding to photovoltaic power and energy storage charging power in the sparse precision matrix is 0.72, while many other elements are close to zero, it is considered that there is a direct causal edge between photovoltaic power and energy storage charging power.

[0083] Based on the constructed time-varying causal graph, the historical observation data of the energy node is subjected to exogenous treatment to obtain a control sample set. Exogenous treatment is to regard the boundary interactive variable as an exogenous variable, and eliminate its endogenous association with other variables. Specifically, for each pair of variables connected by a direct causal edge, by controlling other variables, the part affected only by the boundary interactive variable is extracted as a control sample. For example, for the energy storage charging power, after controlling other variables, the part affected only by the combined heat and power output is extracted as a control sample.

[0084] The propensity score of the boundary interactive variable to the state data in the scheduling sub-domain is calculated, and the propensity score represents the probability of a certain node receiving a specific intervention under the condition of other variables. The propensity score is estimated by logistic regression, and the input is all state variables in the sub-domain, and the output is a binary probability value representing the possibility of receiving intervention. For example, when the boundary interactive variable is the power output of combined heat and power, the estimated propensity score of its influence on the state of the energy storage node may be 0.65.

[0085] The observation data of the control sample set is weighted based on the reciprocal of the calculated propensity score to obtain weighted observation data. The reciprocal of the propensity score is used as a weight to balance the selection bias of the sample. For example, if the propensity score of a certain observation sample is 0.2, the weight is 5, and the influence of the sample in subsequent calculations is adjusted accordingly. Counterfactual intervention is applied to the boundary interaction variable, that is, it is assumed that the value of the boundary variable changes in a certain way, such as an increase of 100 kW in the combined heat and power generation. Based on the weighted observation data, the difference between the power change gradient of the state data in the scheduling sub-domain before and after the intervention is calculated to obtain the time difference component. For example, the power gradient of the energy storage node before the intervention is 10 kW / min, and the power gradient after the intervention is 15 kW / min, then the time difference component is 5 kW / min.

[0086] Based on the time difference component, the control sample set is subjected to second difference to obtain an inter-group difference component. The second difference compares the differences between the intervention group and the control group, and eliminates the influence of common trends. For example, the average time difference component of the intervention group is 5 kW / min, and the average time difference component of the control group is 2 kW / min, then the inter-group difference component is 3 kW / min, indicating the pure causal effect of the boundary interaction variable. The inter-group difference component is output as the causal influence strength, which is used for subsequent cross-domain collaborative scheduling decisions.

[0087] In this embodiment, by setting the coupling relationship parameter between the scheduling sub-domains as the boundary interaction variable and combining the state data within the scheduling sub-domains to construct a multivariate time series, the time sequence dependence relationship between the cross-sub-domain energy transmission and the internal state evolution can be fully captured. By introducing vector autoregression and time-varying precision matrix modeling, fine-grained modeling and high-time-efficiency expression of the cross-domain dynamic coupling characteristics are realized. The time-varying precision matrix is constrained by adaptive lasso regularization, so that the causal graph structure remains sparse and interpretable in high-dimensional states, effectively suppressing noise interference and pseudo-correlation effects. Based on the time-varying causal graph, exogenous processing and propensity weighting calculation are performed, which can distinguish between endogenous variables and external intervention effects, and improve the scientificity and reliability of causal inference.

[0088] In an alternative embodiment,

[0089] The information transmission weight in the boundary interaction protocol is modified according to the causal influence strength to obtain a plurality of enhanced scheduling sub-domains, comprising:

[0090] For each scheduling sub-domain, the sensitivity coefficient of the boundary interaction variable to the state data within the scheduling sub-domain is calculated based on the causal influence strength. The information transmission weight of each boundary interaction variable in the boundary interaction protocol is nonlinearly mapped to an adjustment amplitude through the sensitivity coefficient. The modified information transmission weight is obtained by weighting and fusing the adjustment amplitude and the information transmission weight.

[0091] Based on the modified information transmission weight, the information transmission delay between adjacent scheduling sub-domains is calculated, and the boundary interaction variables in the boundary interaction protocol are weighted and aggregated according to the modified information transmission weight to obtain inter-sub-domain interaction characteristics. The inter-sub-domain interaction characteristics are time-aligned through the information transmission delay and spliced with the state data in the scheduling sub-domain to obtain a plurality of scheduling sub-domains.

[0092] For each divided scheduling sub-domain, the sensitivity coefficient of the boundary interaction variable to the state data in the scheduling sub-domain is calculated based on the previously obtained causal influence strength. The sensitivity coefficient represents the response degree of the state data in the scheduling sub-domain to the change of the boundary interaction variable, and the calculation method is the causal influence strength divided by the change amount of the boundary interaction variable. For example, in a renewable energy-storage collaborative sub-domain, if the combined heat and power increases by 100 kilowatts as a boundary interaction variable, resulting in a change of 3 kilowatts per minute in the charging power gradient of the energy storage, the corresponding sensitivity coefficient is 0.03. The sensitivity coefficient of each state data in the sub-domain to all boundary interaction variables is calculated to form a sensitivity coefficient matrix. In a typical scenario, the dimension of the sensitivity coefficient matrix may be 10x15, indicating the sensitivity of 10 boundary interaction variables to 15 internal state variables.

[0093] The information transmission weight of each boundary interaction variable in the boundary interaction protocol is nonlinearly mapped by the sensitivity coefficient to obtain an adjustment amplitude. The nonlinear mapping is realized by using a sigmoid function to convert the sensitivity coefficient into an adjustment amplitude between 0 and 1. The sigmoid function is applied to the sensitivity coefficient, and appropriate scaling parameters are set so that high sensitivity corresponds to a larger adjustment amplitude and low sensitivity corresponds to a smaller adjustment amplitude. For example, when the sensitivity coefficient is 0.03, the adjustment amplitude obtained by nonlinear mapping may be 0.55; when the sensitivity coefficient is 0.08, the adjustment amplitude may be 0.65. For very low sensitivity coefficients such as 0.001, the adjustment amplitude is close to 0.5, indicating that almost no adjustment is needed.

[0094] The adjustment amplitude and the original information transmission weight are weighted and fused to obtain a modified information transmission weight. The weighted fusion is realized by linear combination, and the modified information transmission weight is equal to the original weight multiplied by an adjustment factor, and the adjustment factor is a first function of the adjustment amplitude. For example, for a boundary interaction variable with an original information transmission weight of 0.7, if the adjustment amplitude is 0.55, the adjustment factor is 1.05, and the modified information transmission weight is 0.735. For all boundary interaction variables, such modification is performed to obtain a set of modified information transmission weights. The information transmission weight reflects the actual influence degree of the boundary interaction variable on the state in the scheduling sub-domain, and more accurately guides the information exchange process.

[0095] Based on the modified information transmission weight, the information transmission delay between adjacent scheduling subdomains is calculated. The information transmission delay refers to the time required for information to propagate from the source subdomain to the target subdomain, and is inversely proportional to the modified information transmission weight. In the specific calculation, the actual transmission delay is obtained by dividing the standard transmission time by the modified information transmission weight. For example, if the standard transmission time is set to 300 milliseconds and the modified information transmission weight is 0.735, the actual transmission delay is about 408 milliseconds. The transmission delay of all interaction channels between subdomains is calculated to form a delay matrix.

[0096] The boundary interaction variables in the boundary interaction protocol are weighted and aggregated according to the modified information transmission weight to obtain the interaction characteristics between subdomains. In the weighted aggregation process, each boundary interaction variable is multiplied by its modified information transmission weight and then summed to obtain the comprehensive interaction characteristics. For example, for the interaction between the combined heat and power subdomain and the renewable energy-storage subdomain, the boundary interaction variables include combined heat and power output, transmission loss rate, and heat-to-power ratio, and the modified information transmission weights are 0.735, 0.628, and 0.591, respectively. Therefore, the comprehensive interaction characteristics are the weighted sum of these three variables according to the corresponding weights.

[0097] The interaction characteristics between subdomains are time-aligned through the information transmission delay, and are spliced with the state data within the scheduling subdomain to obtain an enhanced scheduling subdomain. Time alignment compensates for the information lag caused by transmission delay, ensuring the consistency of interaction data and internal state in time. For each interaction characteristic, time translation is performed according to its transmission delay to align it with the internal state of the receiving subdomain in time. For example, if the transmission delay from the combined heat and power subdomain to the energy storage subdomain is 408 milliseconds, the interaction characteristics of the combined heat and power are time-shifted forward by 408 milliseconds to achieve time alignment. The time-aligned interaction characteristics are spliced with the internal state data of the subdomain to form an extended state vector, which constitutes the enhanced scheduling subdomain.

[0098] In this embodiment, by calculating the sensitivity coefficients of boundary interaction variables based on the causal influence strength, and performing nonlinear mapping and weighted correction on the information transmission weight, the key interaction variables most sensitive to the internal state changes of the scheduling subdomain can be automatically identified and strengthened, realizing fine-grained regulation of information flow. Based on the modified information transmission weight, the information transmission delay calculation and time alignment mechanism are introduced, so that the interaction characteristics between subdomains can be kept synchronized in the time dimension, avoiding the problems of scheduling mismatch and information feedback deviation caused by time lag accumulation. Through the splicing and fusion of the time-aligned interaction characteristics and the internal state data of the scheduling subdomain, the collaborative modeling of cross-domain energy and information in a unified spatiotemporal feature space can be realized, providing more accurate, continuous, and interpretable input basis for subsequent global scheduling optimization and dynamic prediction.

[0099] In an optional embodiment,

[0100] The actual response data of the energy nodes in each enhanced dispatch sub-domain is collected and compared with the expected response data to obtain a response deviation, the response deviation is sequentially matched with historical dispatch scenario data to obtain a time evolution pattern of the response deviation and a future impact trajectory corresponding to the response deviation is predicted, which comprises:

[0101] The actual power output of the energy nodes in each enhanced dispatch sub-domain is collected as actual response data, and a preset target power output is taken as expected response data. The response deviation of the actual response data and the expected response data is calculated, singular spectrum analysis is performed to obtain a reconstruction component and a residual component, the Lyapunov index of the reconstruction component is calculated, and the response deviation sequence is obtained by splicing the reconstruction component and the residual component.

[0102] A probability distribution is constructed by extracting boundary interaction variables from the pre-acquired historical dispatch scenario data, the maximum average difference between the boundary interaction variables of the current dispatch scenario and the probability distribution is calculated, a historical scenario subset is constructed by extracting historical dispatch scenarios with a maximum average difference less than a preset difference threshold, a common latent space representation of the response deviation sequence and historical response deviation sequences in the historical scenario subset is extracted by kernel canonical correlation analysis, the Mahalanobis distance is calculated based on the common latent space representation to determine the sequence similarity, and the historical response deviation sequence with a sequence similarity exceeding a preset similarity threshold is taken as a matching sequence. The matching sequence is subjected to Bayesian structural time series decomposition to extract a trend evolution path and a periodic disturbance component and determine a time evolution pattern, the time evolution pattern is applied to the response deviation sequence, and a multi-step forward deduction is performed by integrated Kalman filtering to obtain the future impact trajectory.

[0103] The actual power output of the energy nodes in each enhanced dispatch sub-domain is collected as actual response data, and a preset target power output is taken as expected response data. In actual application, taking the renewable energy-storage sub-domain of a certain virtual power plant as an example, the actual power output of photovoltaic arrays, energy storage systems and other nodes every 5 minutes within 24 hours is collected to form a time sequence of 288 sampling points, such as the actual output of photovoltaic [0, 0, 0,..., 230, 245, 260,..., 0, 0] kW; at the same time, the target power output sequence in the dispatch plan is obtained, such as [0, 0, 0,..., 240, 250, 260,..., 0, 0] kW. The response deviation of the actual response data and the expected response data is calculated, that is, the difference value of the corresponding time points of the two sequences, to form a response deviation sequence, such as [0, 0, 0,..., -10, -5, 0,..., 0, 0] kW.

[0104] Singular spectrum analysis is performed on the response deviation sequence to extract reconstructed components and residual components. Singular spectrum analysis is implemented through four steps: embedding, decomposition, grouping, and diagonal averaging. The embedding process converts the one-dimensional sequence into a trajectory matrix using a window length of 60 (i.e., 5 hours of data); the decomposition process performs singular value decomposition on the trajectory matrix to obtain singular values and corresponding left and right singular vectors; the grouping process groups the components corresponding to the first 10 singular values as reconstructed components and the remaining components as residual components; and the diagonal averaging process converts the grouped matrix back to a one-dimensional sequence. For example, for the aforementioned response deviation sequence, the reconstructed components may be [-2, -3, -5,..., -9, -6, -1,..., 0, 0] kW and the residual components may be [2, 3, 5,..., -1, 1, 1,..., 0, 0] kW.

[0105] The Lyapunov exponent of the reconstructed components is calculated to assess the degree of chaos and predictability of the sequence. The Lyapunov exponent calculation process involves phase space reconstruction and estimation of the divergence rate of neighboring trajectories. Phase space reconstruction is performed using an embedding dimension of 6 and a time delay of 5, tracking the evolution of neighboring trajectory points and calculating the exponential divergence rate. For stable systems, the exponent is typically negative; for chaotic systems, the exponent is positive. For example, the calculated Lyapunov exponent is 0.05, indicating that the system has a slight chaotic characteristic. The calculated Lyapunov exponent is concatenated with the reconstructed components and residual components to form an enhanced response deviation sequence for subsequent analysis.

[0106] The boundary interaction variables are extracted from the pre-acquired historical scheduling scenario data to construct a probability distribution. The historical data contains scheduling records from the past 6 months, with 288 sampling points per day, totaling approximately 50,000 historical scenarios. The boundary interaction variables of each scenario are extracted, such as combined heat and power output, regional electricity price, and grid frequency, to construct a multi-dimensional probability distribution. The probability distribution can be constructed using the kernel density estimation method, using a Gaussian kernel function with a bandwidth of 0.15. The maximum average discrepancy between the boundary interaction variables of the current scheduling scenario and the probability distribution is calculated. The maximum average discrepancy is calculated by calculating the difference between the boundary variables of the current scenario and the historical distribution in each dimension, and taking the maximum value. For example, if the current combined heat and power output is 200 kW and the historical distribution has a mean of 180 kW and a standard deviation of 30 kW, the difference in this dimension is 0.67 standard deviations; similarly, calculate the differences in other dimensions, and take the maximum value as the maximum average discrepancy.

[0107] The historical scheduling scenarios with maximum average difference less than a preset difference threshold are extracted to construct a historical scenario subset. The preset difference threshold is set to 1.0 standard deviation, and the historical scenarios with a difference less than the threshold are selected. In actual application, about 2000 scenarios meeting the conditions can be selected from 50,000 historical scenarios to form the historical scenario subset. The common latent space representation of the response deviation sequence and the historical response deviation sequences in the historical scenario subset is extracted by kernel canonical correlation analysis. The kernel canonical correlation analysis uses a Gaussian kernel function, and the kernel parameter is set to 5.0 to extract a 10-dimensional common latent representation. This method finds the maximum correlation direction of the two groups of data in the nonlinear mapping space, effectively capturing the complex nonlinear relationship.

[0108] The Mahalanobis distance based on the common latent space representation is calculated to determine the sequence similarity. The Mahalanobis distance considers the variance and covariance of each dimension, and can more accurately measure the similarity in multi-dimensional space. The Mahalanobis distance of the current response deviation sequence and each sequence in the historical scenario subset in the latent space is calculated, and the smaller the distance, the higher the similarity. The historical response deviation sequences with a sequence similarity exceeding a preset similarity threshold are selected as matching sequences. The preset similarity threshold is 0.75, and if the similarity corresponding to the Mahalanobis distance exceeds 0.75, the sequence is included in the matching sequence set. In actual application, about 200 high-similarity matching sequences can be further selected from the 2000 historical scenarios.

[0109] The matching sequences are decomposed by Bayesian structural time series to extract the trend evolution path and periodic disturbance component and determine the time evolution pattern. The Bayesian structural time series decomposition divides the sequence into three parts: trend item, seasonal item and residual item. The trend item is represented by a local linear model, the seasonal item is represented by a Fourier series, and the residual item is represented by an autoregressive model. The decomposition parameters are obtained by Bayesian inference, and Gibbs sampling is used for 2000 iterations, with the first 500 iterations discarded as preheating. The trend evolution path obtained by decomposition reflects the long-term change direction of the system, and the periodic disturbance component reflects the short-term fluctuation rule, which together constitute the time evolution pattern.

[0110] The time evolution pattern is applied to the current response deviation sequence, and a multi-step forward deduction is performed by integrated Kalman filtering to obtain the future impact trajectory. The integrated Kalman filter is composed of multiple Kalman filters with different parameters, and the state vector of each filter contains 6 components. The process noise covariance is set as a diagonal matrix, and the observation noise is set as a scalar. Each filter independently performs state prediction and update, and the final result is obtained by weighted averaging, with the weights being inversely proportional to the prediction error of each filter. When performing multi-step forward deduction, the prediction step is set to 72 time points (i.e. 6 hours), and each prediction result is used as the input of the next time step to gradually deduce the complete future impact trajectory.

[0111] In this embodiment, by comparing the actual power output with the target power output, combining singular spectrum analysis to extract the reconstruction component and residual component in response deviation, and calculating Lyapunov index to characterize the chaotic characteristics and stability level of system response, the potential unstable mode and nonlinear disturbance source in the system can be identified in the time series dimension. Through the scene screening mechanism based on maximum average difference, samples similar to the current running state distribution are extracted from the historical scheduling scenes, which can effectively eliminate the matching deviation caused by scene drift, ensure the accuracy and applicability of historical experience migration, use kernel canonical correlation analysis to extract the common latent space representation of response deviation sequence and historical sequence, and measure the similarity through Mahalanobis distance, which realizes high-precision matching in multi-dimensional nonlinear feature space. Compared with traditional correlation analysis methods, it has stronger discriminability and generalization ability. Combined with the multi-step forward deduction mechanism of Bayesian structure time series decomposition and integrated Kalman filter, the future power response deviation can be dynamically predicted based on the fusion of trend evolution and periodic disturbance components, which not only improves the timeliness and robustness of prediction, but also provides an interpretable future trend judgment under uncertainty.

[0112] In an alternative embodiment,

[0113] When the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, the corresponding affected scheduling sub-domain is determined and boundary coupling relationship adjustment is performed, combined with the forward iteration solution of the rolling time domain optimization method, to obtain the modified scheduling instruction and execute, including:

[0114] It is judged whether the deviation amplitude of each prediction time in the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set. If it exceeds, the corresponding exceeding time is extracted as the prediction time, the energy node identifier corresponding to the prediction time is obtained, and the affected scheduling sub-domain is determined based on the energy node identifier. The boundary coupling relationship between the affected scheduling sub-domain and the adjacent scheduling sub-domain is extracted, the Jacobian matrix corresponding to the power exchange constraint in the boundary coupling relationship is constructed, and the directional derivative is calculated to obtain the sensitivity coefficient. Based on the sensitivity coefficient, the power exchange constraint is modified to obtain the optimized boundary coupling relationship.

[0115] The optimized boundary coupling relationship is used as a consensus constraint, the prediction time domain and control time domain of the rolling time domain optimization are set, the objective function is constructed in the prediction time domain, and the Pareto optimal control sequence in the control time domain is obtained by solving through the alternating direction multiplier method. The control quantity of the first period in the Pareto optimal control sequence is extracted as the current execution action. The current execution action is converted into a scheduling instruction and issued to the energy nodes in the affected scheduling sub-domain. The actual execution data of the energy nodes is collected to update the actual response data and reconstruct the prediction time domain for forward iteration to obtain the modified scheduling instruction and execute.

[0116] It is judged whether the deviation amplitude of each prediction time in the future influence trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set. For each prediction time, the deviation amplitude is compared with the pre-set allowable fluctuation range. The allowable fluctuation range is set according to the energy node type and system stability requirements, such as the allowable fluctuation range of photovoltaic power generation is ±5% of the rated power, the allowable fluctuation range of wind power generation is ±8%, and the allowable fluctuation range of energy storage system is ±10%. For example, the rated power of a certain photovoltaic node is 1000 kW, and its allowable fluctuation range is ±50 kW; if the prediction deviation of the future influence trajectory at the 30th minute is 62 kW, which exceeds the allowable range, the 30th minute is extracted as the exceeding time. If the allowable range is exceeded, the corresponding exceeding time is extracted as the prediction time. In actual application, 10 exceeding times may be extracted from 72 prediction points, such as the 30th, 45th, 60th, 75th, 90th, 120th, 150th, 180th, 240th and 360th minutes.

[0117] The energy node identifier corresponding to the prediction time is obtained, and the affected scheduling sub-domain is determined based on the energy node identifier. For each prediction time, the affected scheduling sub-domain is determined according to the energy node identifier corresponding to the deviation. For example, the node identifier corresponding to the 30th minute deviation is "PV-001" (photovoltaic array), which belongs to the renewable energy-storage sub-domain, so the sub-domain is determined as the affected scheduling sub-domain. Similar processing is performed for all exceeding times, which may involve multiple scheduling sub-domains. The boundary coupling relationship between the affected scheduling sub-domain and the adjacent scheduling sub-domain is extracted. The boundary coupling relationship includes power exchange constraints, time delay constraints and energy balance constraints, etc. Taking the boundary between the renewable energy-storage sub-domain and the combined heat and power sub-domain as an example, the power exchange constraints include parameters such as the upper limit of transmission capacity of 500 kW, the minimum power exchange amount of 50 kW and the transmission loss rate of 3%.

[0118] The Jacobian matrix corresponding to the power exchange constraint in the boundary coupling relationship is constructed, and the sensitivity coefficient is calculated by calculating the directional derivative. The Jacobian matrix describes the partial derivative relationship of the power exchange constraint to each variable, which is constructed by numerical difference method. Taking the renewable energy-storage sub-domain as an example, its power exchange constraint involves variables such as its own power generation, energy storage state and exchange power with the combined heat and power sub-domain. For each variable, the rate of change of the constraint function is calculated by small perturbation to construct the Jacobian matrix. For example, if the energy storage power is increased by 10 kW and the power exchange constraint is observed to change by 5 kW, the corresponding sensitivity is 0.5. The directional derivative of the constraint function is calculated in the direction of the future influence trajectory to obtain the sensitivity coefficient. The sensitivity coefficient represents the sensitivity of the boundary constraint in the direction of the prediction deviation, and the larger the value is, the more sensitive the constraint is.

[0119] The power exchange constraints are modified based on the sensitivity coefficients to obtain the optimized boundary coupling relationship. The modification method is to adjust the constraint parameters according to the sensitivity coefficients. The constraint parameters with high sensitivity are adjusted by a larger amount, and the constraint parameters with low sensitivity are adjusted by a smaller amount. The specific adjustment amount is proportional to the sensitivity coefficient and inversely proportional to the importance of the constraint. For example, for the upper limit of the transmission capacity, if the sensitivity coefficient is 0.8 and the importance is 0.9, the adjustment amount is about 8.9%, and the modified upper limit of the transmission capacity is 545 kilowatts; for the minimum power exchange amount, if the sensitivity coefficient is 0.3 and the importance is 0.7, the adjustment amount is about 4.3%, and the modified minimum power exchange amount is 48 kilowatts. Similarly, all constraint parameters are modified to obtain the optimized boundary coupling relationship.

[0120] The optimized boundary coupling relationship is used as a consensus constraint, and the prediction time domain and control time domain of the rolling time domain optimization are set. The prediction time domain is set to 6 hours, containing 72 time points; the control time domain is set to 30 minutes, containing 6 time points. The objective function is constructed within the prediction time domain. The objective function considers three aspects of economy, stability and environmental protection: economy is embodied in the minimization of operation cost, including generation cost, energy storage loss and dispatch adjustment cost; stability is embodied in the minimization of power deviation, calculating the square difference of predicted power and target power; environmental protection is embodied in the minimization of carbon emissions, calculating the total emissions according to the carbon emission coefficients of different energy types. The three objectives are combined into a composite objective function by weighted summation, and the weights are dynamically adjusted according to the dispatch strategy, with typical values of economy 0.4, stability 0.4 and environmental protection 0.2.

[0121] The Pareto optimal control sequence within the control time domain is obtained by the alternating direction multiplier method. The alternating direction multiplier method decomposes the original problem into subproblems, and the overall optimal solution is obtained by iteration. The algorithm parameters are set as follows: the augmented Lagrange multiplier is 1.5, the step size is 0.01, the maximum iteration number is 300, and the convergence threshold is 0.001. For each subproblem, the quadratic programming method is used to solve, considering various constraint conditions, including physical limitations of energy nodes, optimized boundary coupling relationship and system balance constraints, etc. Through multiple iterations, the algorithm converges to the Pareto optimal solution, which is the best balance point among the three objectives. The Pareto optimal control sequence contains control variables at each time point within the control time domain, such as generation power, energy storage charging and discharging power, and load adjustment amount, etc.

[0122] The control quantity of the first time period in the Pareto optimal control sequence is extracted as the current execution action. For example, the control quantity of the first time period in the control sequence is the control quantity of the next 5 minutes, including the photovoltaic inverter power setting value 250 kW, the energy storage charging power 80 kW, and the adjustable load power 150 kW. The current execution action is converted into a scheduling instruction and issued to the energy nodes in the affected scheduling sub-domain. The conversion process involves instruction formatting, safety checks, and priority allocation. The formatted instruction contains information such as target node identification, execution time, control quantity value, and validity period. The instruction is issued to the target nodes, such as photovoltaic inverters, energy storage management systems, and load controllers, through a communication network.

[0123] The actual execution data of the energy nodes is collected to update the actual response data. The collection frequency is set to once per minute, including parameters such as power, voltage, and frequency. For example, 5 minutes after the instruction is issued, the actual output power of the photovoltaic is 248 kW, the energy storage charging power is 82 kW, and the adjustable load power is 148 kW. The actual execution data is compared with the expected data to calculate the new response deviation and update the actual response data. The prediction time domain is reconstructed for forward iteration to obtain the modified scheduling instruction and execute it. The prediction time domain is shifted forward by one control interval (5 minutes), the time window is discarded, and a new prediction period is added. Based on the updated response data and boundary conditions, the future impact trajectory prediction, boundary coupling relationship optimization, and rolling horizon optimization are performed again to obtain a new control sequence. The control quantity of the first time period in the new control sequence is extracted as the modified scheduling instruction and executed.

[0124] In this embodiment, by monitoring the deviation amplitude of each prediction time in the future impact trajectory and comparing it with the allowable fluctuation range in the coupling relationship set, the potential instability period and the affected sub-domain can be identified in advance, the prediction-driven active scheduling intervention is realized, the sensitivity coefficients calculated by the Jacobian matrix and directional derivative are used to accurately quantify the influence of the boundary energy interaction on the system stability, the adaptive optimization of the boundary coupling relationship is realized at the control level, and the modified boundary coupling relationship is used as a consensus constraint to solve the rolling horizon optimization framework and the alternating direction multiplier method, which can dynamically balance the system performance, energy utilization rate, and constraint stability in the prediction time domain, obtain the Pareto optimal control sequence, and ensure the global coordination and optimal allocation among different energy targets.

[0125] Figure 2 The sensitivity optimization multi-energy coordination flowchart for the virtual power plant multi-energy collaborative scheduling and real-time optimization control method of the embodiment of the application.

[0126] In a second aspect of the embodiment of the application, a virtual power plant multi-energy collaborative scheduling and real-time optimization control system is provided, which comprises:

[0127] The first unit is configured to acquire state data of energy resources and identify coupling relationships between energy nodes to obtain a coupling relationship set.

[0128] The second unit is configured to construct a space-time correlation graph according to the coupling relationship set and calculate an embedding representation of each energy node, identify a strong space-time correlation node pair based on the embedding representation, divide energy nodes containing the strong space-time correlation node pair to the same scheduling sub-domain through spectral clustering, configure a boundary interaction protocol for each scheduling sub-domain, identify the causal influence strength of an interaction variable on the scheduling sub-domain in the boundary interaction protocol through a causal inference method, and modify the information transmission weight in the boundary interaction protocol according to the causal influence strength to obtain a plurality of enhanced scheduling sub-domains.

[0129] The third unit is configured to acquire actual response data and expected response data of energy nodes in each enhanced scheduling sub-domain, compare the actual response data with the expected response data to obtain a response deviation, perform sequence matching on the response deviation and historical scheduling scene data to obtain a time sequence evolution mode of the response deviation and predict a future influence trajectory corresponding to the response deviation, determine an affected scheduling sub-domain when the future influence trajectory exceeds an allowable fluctuation range corresponding to the coupling relationship set, and perform boundary coupling relationship adjustment, and perform forward iteration solving through a rolling time domain optimization method to obtain a modified scheduling instruction and execute the same.

[0130] In a third aspect, an electronic device is provided, including:

[0131] A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0132] In a fourth aspect, a computer-readable storage medium is provided, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0133] The present application can be a method, device, system and / or computer program product. The computer program product can include a computer readable storage medium having computer readable program instructions stored therein, which are used to perform various aspects of the present application.

[0134] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement to part or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for multi-energy coordinated dispatch and real-time optimization control of virtual power plants, characterized in that, include: Acquire energy resource status data and identify coupling relationships between energy nodes to obtain a set of coupling relationships; A spatiotemporal correlation graph is constructed based on the set of coupling relationships, and the embedded representation of each energy node is calculated. Based on the embedded representation, strong spatiotemporal correlation node pairs are identified, and energy nodes containing strong spatiotemporal correlation node pairs are divided into the same scheduling subdomain through spectral clustering. A boundary interaction protocol is configured for each scheduling subdomain. The causal inference method is used to identify the causal influence strength of the interaction variables in the boundary interaction protocol on the scheduling subdomains. Based on the causal influence strength, the information transmission weight in the boundary interaction protocol is corrected to obtain multiple enhanced scheduling subdomains. The actual response data of energy nodes in each enhanced scheduling subdomain is collected and compared with the expected response data to obtain the response deviation. The response deviation is then matched with historical scheduling scenario data to obtain the temporal evolution pattern of the response deviation and predict the future impact trajectory corresponding to the response deviation. When the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, the corresponding affected scheduling subdomain is determined and the boundary coupling relationship is adjusted. The forward iterative solution is then performed using the rolling time domain optimization method to obtain the corrected scheduling instruction and execute it. By acquiring energy resource status data and identifying the coupling relationships between energy nodes, a set of coupling relationships is obtained, including: Collect power output data, energy conversion data, and power exchange data between resources for each energy resource within a continuous time window as status data; A power transfer matrix is ​​constructed based on the power exchange data between resources in the state data, and singular value decomposition is performed to obtain singular values ​​and corresponding singular vectors. The dominant path is identified from the singular vectors corresponding to singular values ​​that exceed the preset singular value threshold, and adjacent energy node pairs in the dominant path are selected as candidate node pairs. For each node pair in the candidate node pair, the power transfer between the source node and the target node is calculated as a percentage of the total power output of the source node based on the power output data. If the percentage exceeds a preset threshold, the current node pair is marked as a coupled node pair. Extract the power transfer amount and transfer delay corresponding to each coupled node pair from the power transfer matrix, extract the conversion efficiency boundary from the energy conversion data, and encapsulate the power transfer amount, conversion efficiency boundary and transfer delay as coupling relationship parameters and node pair identifiers into a set of coupling relationships. Based on the set of coupling relationships, a spatiotemporal correlation graph is constructed and the embedding representation of each energy node is calculated. Based on the embedding representation, strong spatiotemporal correlation node pairs are identified, and energy nodes containing strong spatiotemporal correlation node pairs are grouped into the same scheduling subdomain through spectral clustering, including: Extract the coupled node pairs from the set of coupling relationships as the edge topology, use the energy nodes as the vertex topology, construct the spatiotemporal correlation graph using the pre-acquired coupling relationship parameters as edge attributes, and calculate the power change gradient and phase characteristics corresponding to the state data. For each vertex in the spatiotemporal correlation graph, an initial representation is constructed based on the power change gradient, phase features, and edge attributes. The initial representations of adjacent energy nodes are multiplied by the power transfer amount using the Hadamard product to obtain the power representation. This power representation is then weighted element-wise by the inverse of the pre-acquired transfer delay. A nonlinear mapping is applied to obtain the projection representation. Based on the cosine distance between the projection representation and the initial representation, normalized aggregation and skip connections are performed to obtain the embedding representation. Calculate the bilinear metric between embedded representations, extract the bidirectional power flow ratio of the coupled node pair within a preset time window, and mark the coupled node pair as a strongly spatiotemporally related node pair when the bidirectional power flow ratio exceeds a preset power flow threshold and the bilinear metric exceeds a preset similarity threshold. Based on the embedding representation, an affinity matrix is ​​constructed and a diffusion kernel is applied to obtain a diffusion matrix. Based on the diffusion matrix, spectral decomposition is performed to obtain spectral embedding, and spectral clustering is performed on the strong spatiotemporal correlated node pairs. According to the clustering results, energy nodes containing strong spatiotemporal correlated node pairs are divided into the same scheduling subdomain.

2. The method according to claim 1, characterized in that, A boundary interaction protocol is configured for each scheduling subdomain, and the strength of the causal influence of the interaction variables in the boundary interaction protocol on the scheduling subdomains is identified using a causal inference method, including: For each scheduling subdomain, the energy node coupling relationship parameters between the current scheduling subdomain and adjacent scheduling subdomains are used as boundary interaction variables and a boundary interaction protocol is configured. A multivariate time series of the boundary interaction variables and the state data within the scheduling subdomain is established and vector autoregression is performed to fit and extract the residual sequence. A time-varying precision matrix is ​​calculated based on the residual sequence. An adaptive lasso regularization constraint is applied to the time-varying precision matrix and a sparse precision matrix is ​​obtained by iteratively solving the coordinate descent problem. Variable pairs corresponding to non-zero elements are extracted from the sparse precision matrix as direct causal edges to construct a time-varying causal graph. Based on the direct causal edges in the time-varying causal graph, the historical observation data corresponding to the energy node are exogenously processed to obtain a control sample set. The tendency score of the boundary interaction variable on the state data in the scheduling subdomain is calculated. The observation data of the control sample set is weighted based on the reciprocal of the tendency score to obtain weighted observation data. Counterfactual intervention is applied to the boundary interaction variable, and the difference between the power change gradients corresponding to the state data in the scheduling subdomain before and after the intervention is calculated based on the weighted observation data to obtain the time difference component. The control sample set is then subjected to a second difference based on the time difference component to obtain the inter-group difference component, and the inter-group difference component is output as the causal influence strength.

3. The method according to claim 1, characterized in that, Based on the correction of the information transmission weights in the boundary interaction protocol according to the causal influence strength, multiple enhanced scheduling subdomains are obtained, including: For each scheduling subdomain, the sensitivity coefficient of the boundary interaction variable to the state data within the scheduling subdomain is calculated based on the causal influence intensity. The information transmission weight of each boundary interaction variable in the boundary interaction protocol is nonlinearly mapped using the sensitivity coefficient to obtain the adjustment range. The adjustment range is then weighted and fused with the information transmission weight to obtain the corrected information transmission weight. Based on the corrected information transmission weight, the information transmission delay between adjacent scheduling subdomains is calculated. The boundary interaction variables in the boundary interaction protocol are weighted and aggregated according to the corrected information transmission weight to obtain the inter-subdomain interaction features. The inter-subdomain interaction features are time-aligned using the information transmission delay and concatenated with the state data within the scheduling subdomain to obtain multiple enhanced scheduling subdomains.

4. The method according to claim 1, characterized in that, The actual response data of energy nodes within each enhanced scheduling subdomain is collected and compared with the expected response data to obtain the response deviation. The response deviation is then matched with historical scheduling scenario data to obtain the temporal evolution pattern of the response deviation and predict the future impact trajectory corresponding to the response deviation, including: The actual power output of energy nodes in each enhanced scheduling subdomain is collected as actual response data, and the preset target power output is used as expected response data. The response deviation between the actual response data and the expected response data is calculated and singular spectrum analysis is performed to obtain the reconstruction component and residual component. The Lyapunov exponent of the reconstruction component is calculated and concatenated with the reconstruction component and residual component to obtain the response deviation sequence. Boundary interaction variables are extracted from pre-acquired historical scheduling scenario data to construct a probability distribution. The maximum average difference between the boundary interaction variables of the current scheduling scenario and the probability distribution is calculated. Historical scheduling scenarios with a maximum average difference less than a preset difference threshold are extracted to construct a subset of historical scenarios. Kernel canonical correlation analysis is used to extract the common latent space representation of the response deviation sequence and the historical response deviation sequences in the subset of historical scenarios. Mahalanobis distance is calculated based on the common latent space representation to determine the sequence similarity. Historical response deviation sequences with a sequence similarity exceeding a preset similarity threshold are used as matching sequences. Bayesian structure time series decomposition is performed on the matching sequences to extract the trend evolution path and periodic disturbance components and determine the temporal evolution mode. The temporal evolution mode is applied to the response deviation sequence and multi-step forward extrapolation is performed through integrated Kalman filtering to obtain the future influence trajectory.

5. The method according to claim 1, characterized in that, When the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, the corresponding affected scheduling subdomain is determined and boundary coupling relationship adjustment is performed. A forward iterative solution is then obtained using a rolling time-domain optimization method to obtain and execute the corrected scheduling instruction, which includes: Determine whether the deviation amplitude of each predicted moment in the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set. If it does, extract the corresponding exceeding moment as the prediction moment, obtain the energy node identifier corresponding to the prediction moment, and determine the affected scheduling subdomain based on the energy node identifier. Extract the boundary coupling relationship between the affected scheduling subdomain and the adjacent scheduling subdomain, construct the Jacobian matrix corresponding to the power exchange constraint in the boundary coupling relationship, and calculate the directional derivative to obtain the sensitivity coefficient. Based on the sensitivity coefficient, correct the power exchange constraint to obtain the optimized boundary coupling relationship. Using the optimized boundary coupling relationship as a consensus constraint, a prediction time domain and a control time domain are set for rolling time domain optimization. An objective function is constructed in the prediction time domain and solved using the alternating direction multiplier method to obtain the Pareto optimal control sequence in the control time domain. The control quantity of the first time period in the Pareto optimal control sequence is extracted as the current execution action. The current execution action is converted into a scheduling instruction and issued to the energy nodes in the affected scheduling subdomain. The actual execution data of the energy nodes is collected to update the actual response data, and the prediction time domain is reconstructed for forward iteration to obtain the corrected scheduling instruction and execute it.

6. A virtual power plant multi-energy collaborative dispatch and real-time optimization control system, used to implement the method of any one of claims 1-5, characterized in that, include: The first unit is used to acquire energy resource status data and identify the coupling relationships between energy nodes to obtain a set of coupling relationships. The second unit is used to construct a spatiotemporal correlation graph based on the set of coupling relationships and calculate the embedded representation of each energy node. Based on the embedded representation, it identifies strong spatiotemporal correlation node pairs and divides the energy nodes containing strong spatiotemporal correlation node pairs into the same scheduling subdomain through spectral clustering. It configures a boundary interaction protocol for each scheduling subdomain, identifies the causal influence strength of the interaction variables in the boundary interaction protocol on the scheduling subdomains through a causal inference method, and corrects the information transmission weight in the boundary interaction protocol according to the causal influence strength to obtain multiple enhanced scheduling subdomains. The third unit is used to collect the actual response data of energy nodes in each enhanced scheduling subdomain and compare it with the expected response data to obtain the response deviation. The response deviation is then matched with historical scheduling scenario data to obtain the temporal evolution pattern of the response deviation and predict the future impact trajectory corresponding to the response deviation. When the future impact trajectory exceeds the allowable fluctuation range corresponding to the coupling relationship set, the corresponding affected scheduling subdomain is determined and the boundary coupling relationship is adjusted. The forward iterative solution is performed using the rolling time domain optimization method to obtain the corrected scheduling instruction and execute it.

7. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 5.

8. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 5.

Images