Control method of port multi-energy supply system based on automatic inertia response

Abstract

The present invention relates to the technical field of port multi-energy supply systems. The present invention provides a port multi-energy supply system control method based on automatic inertia response, the method comprising the following steps: collecting data and performing data preprocessing; performing grid-load interaction pattern recognition and feature extraction through the construction and multi-level fusion of a time-varying grid-load interaction graph; constructing a response characteristic tensor and a time series prediction model and performing prediction; obtaining a multi-scale feature matrix through wavelet transform and performing optimization to obtain an optimized control sequence and a response strategy matrix; generating a control instruction matrix based on the optimized control sequence and the response strategy matrix, and performing stability analysis and real-time correction. The present invention realizes the intelligent control of the port multi-energy supply system participating in grid-load interaction by constructing a complete data processing, analysis and control closed loop, solves the technical difficulties in the port multi-energy supply system participating in power grid regulation, and enables the system to perform inertia response control more intelligently and efficiently.

Description

Technical Field

[0001] The present invention relates to a port multi-energy supply system, and in particular to a port multi-energy supply system control method based on automatic inertia response. Background Art

[0002] As important hubs in the marine economy, ports have high-power, highly volatile energy demands. With the accelerated electrification and intelligentization of ports, and the increasing proportion of green energy, port multi-energy supply systems are becoming increasingly prominent in the power system. By participating in grid frequency and peak regulation, these systems not only provide important ancillary services but also enhance their own economic benefits. Especially in the context of the ongoing reduction in power system inertia, leveraging the rapid response capabilities of port multi-energy supply systems to participate in system inertia response is crucial for maintaining grid frequency stability.

[0003] Currently, research on the participation of multi-energy supply systems in port grid regulation primarily focuses on traditional primary and secondary frequency regulation. Common control methods include classical control based on PID (Proportional-Integral-Derivative), model predictive control, and adaptive control. These methods typically employ a single data processing model, achieving control objectives through simple feature extraction and state estimation. Data analysis primarily relies on basic statistical methods and traditional machine learning algorithms, such as regression analysis and support vector machines.

[0004] However, these solutions have several shortcomings when addressing the inertia response problem of multi-energy supply systems in ports: First, the extraction of grid-load interaction features often uses static graph structures or simple dynamic graph models, which makes it difficult to accurately capture the time-varying characteristics of the system. Second, existing grid-load interaction modeling methods fail to fully consider the impact of time decay, resulting in long-term dependency issues in long-sequence modeling. Third, during the feature fusion process, there is a lack of an effective integration mechanism for features at different levels, which affects the accurate representation of system status. Finally, the optimization algorithm has low computational efficiency when dealing with large-scale constrained optimization problems, making it difficult to meet the requirements of real-time control. These technical issues have seriously restricted the effectiveness of multi-energy supply systems in port grid inertia response. Summary of the Invention

[0005] In order to solve the above technical problems, the present invention provides a port multi-energy supply system control method based on automatic inertia response.

[0006] The technical solution adopted by the present invention is as follows: A port multi-energy supply system control method based on automatic inertia response, comprising the following steps: S1, collecting time series data of the port multi-energy supply system, constructing a synchronous time series tensor based on the time series data, extracting dynamic characteristics according to the synchronous time series tensor, and obtaining a standardized system characteristic matrix and a dynamic response capability index, wherein the time series data includes real-time power data of the system, operating status data of port equipment, historical operating data and feedback data; S2, constructing a time-varying network-load interaction diagram based on the standardized system characteristic matrix, the dynamic response capability index, pre-stored power grid operating parameters and load fluctuation data, extracting features from the time-varying network-load interaction diagram to generate a spatiotemporal feature matrix, and using the spatiotemporal feature matrix to perform pattern recognition to obtain a network Load-grid interaction feature vector and key mode matrix are obtained through multi-level fusion to obtain fused feature vector and feature mapping matrix; S3, constructing response characteristic tensor based on dynamic response capability index, grid-load interaction feature vector and fused feature vector, constructing time series prediction model based on response characteristic tensor, pre-stored system operation constraints and historical response data and performing prediction to obtain inertia response capability matrix; S4, based on inertia response capability matrix and key mode matrix, obtaining multi-scale feature matrix through wavelet transform and optimizing to obtain optimized control sequence and response strategy matrix; S5, generating control instruction matrix of system based on optimized control sequence and response strategy matrix, performing stability analysis and real-time correction on said control instruction matrix to obtain corrected control instruction sequence and execution feedback data.

[0007] According to one aspect of the present application, step S1 includes: collecting time series data of the port multi-energy supply system, and constructing a synchronous time series tensor based on the time series data; calculating a data quality assessment matrix based on the synchronous time series tensor; using the data quality assessment matrix to detect and correct anomalies in the synchronous time series tensor to obtain a cleaned data tensor; using a sliding window method to extract dynamic characteristics of the cleaned data tensor to obtain a standardized system characteristic matrix and a dynamic response capability indicator.

[0008] According to one aspect of the present application, the steps of anomaly detection and correction specifically include: obtaining a synchronous time series tensor and a data quality assessment matrix, calculating the first-order difference and second-order difference values ​​of the time series, and generating a time series feature vector; using the time series feature vector to perform sliding window statistical operations to generate a time domain anomaly score; performing multi-scale wavelet decomposition on the synchronous time series tensor, calculating the energy distribution characteristics of the decomposition coefficients, and generating a frequency domain anomaly score; calculating the local density distribution and distance distribution of sample points based on the synchronous time series tensor to generate a spatial domain anomaly score; inputting the time domain anomaly score, frequency domain anomaly score, and spatial domain anomaly score into the evidence theory fusion framework, calculating the comprehensive confidence, and generating an anomaly comprehensive score; using the anomaly comprehensive score to correct the abnormal data in the synchronous time series tensor to generate a cleaned data tensor.

[0009] According to one aspect of the present application, step S2 specifically includes: constructing a time-varying weight function based on the standardized system characteristic matrix, dynamic response capability index, power grid operation parameters and load fluctuation data, calculating the edge weights between network nodes based on the time-varying weight function, constructing a time-varying grid-load interaction graph using the edge weights, and calculating the interaction intensity matrix based on the time-varying grid-load interaction graph; constructing a spatiotemporal attention model based on the time-varying grid-load interaction graph and the interaction intensity matrix; using the spatiotemporal attention model to extract features of the time-varying grid-load interaction graph to generate a spatiotemporal feature matrix and an attention weight matrix; obtaining the spatiotemporal feature matrix and the attention weight matrix to construct a dynamic density clustering model; using the dynamic density clustering model to perform cluster analysis on the spatiotemporal feature matrix, tracking the time evolution characteristics of the clustering results, and obtaining a grid-load interaction feature vector, a key pattern matrix and a pattern evolution matrix; constructing a three-layer feature fusion network based on the grid-load interaction feature vector, the key pattern matrix and the pattern evolution matrix, performing multi-level fusion, and obtaining a fusion feature vector and a feature mapping matrix.

[0010] According to one aspect of the present application, the steps of generating a spatiotemporal feature matrix and an attention weight matrix specifically include: constructing a multi-head attention calculation unit containing a time attenuation factor based on a time-varying network load interaction graph and an interaction intensity matrix to generate an attenuated attention matrix; performing a graph convolution operation using the attenuated attention matrix to extract local features of the nodes and generate a graph feature matrix; inputting the graph feature matrix into a bidirectional gated recurrent unit to extract timing dependency features and generate a timing feature vector; performing a feature fusion operation on the graph feature matrix and the timing feature vector to generate a spatiotemporal feature matrix and an attention weight matrix.

[0011] According to one aspect of the present application, the steps of multi-level fusion specifically include: obtaining a network-load interaction feature vector, constructing an autoencoder network, performing feature dimensionality reduction to obtain compressed features, reconstructing features to obtain reconstructed features, calculating reconstruction errors, and outputting underlying fusion features; obtaining a key pattern matrix and underlying fusion features, constructing a dynamic convolution kernel, performing multi-scale convolution operations, pooling processing to obtain an intermediate feature map, extracting significant features, and outputting middle-level fusion features; obtaining a pattern evolution matrix and middle-level fusion features, constructing an evolution feature extractor, calculating the state transition probability of the time series, extracting evolution pattern features, fusing static and dynamic features, and outputting high-level fusion features; obtaining bottom-level fusion features, middle-level fusion features, and high-level fusion features, constructing a multi-layer optimization objective function, calculating the loss function value of each layer of features, performing gradient descent optimization, updating the fusion weights, and outputting the fusion feature vector and feature mapping matrix.

[0012] According to one aspect of the present application, step S3 includes: obtaining a dynamic response capability index, a grid-load interaction characteristic vector and a fusion characteristic vector, constructing a nonlinear response model, calculating response characteristic parameters using the nonlinear response model, and constructing a response characteristic tensor based on the response characteristic parameters; constructing a multidimensional evaluation function based on the response characteristic tensor and pre-stored system operation constraints, calculating sub-item indicators using the multidimensional evaluation function, performing weighted calculation on the sub-item indicators, and outputting a comprehensive evaluation indicator matrix; constructing a time series prediction model using the comprehensive evaluation indicator matrix and pre-stored historical response data, performing time series prediction based on the time series prediction model, and outputting an inertia response capability matrix.

[0013] According to one aspect of the present application, the steps of constructing a response characteristic tensor specifically include: obtaining a dynamic response capability index and a grid-load interaction characteristic vector, calculating a hidden layer state transfer matrix and initial neuron weights based on a recursive neural network structure, and generating a basic model parameter set; using the basic model parameter set and pre-stored system operation history data to construct a piecewise function model, calculating the forward dead zone threshold and the reverse dead zone threshold, and generating a dead zone characteristic matrix; based on the basic model parameter set and the dead zone characteristic matrix, according to the nonlinear response function, calculating the response amplitude attenuation coefficient and the phase offset parameter, and generating an initial characteristic parameter set; constructing a three-dimensional characteristic mapping framework based on the initial characteristic parameter set, calculating the time dimension feature mapping value, the frequency dimension transformation coefficient, and the power dimension mapping parameter, and generating a response characteristic tensor.

[0014] According to one aspect of the present application, step S4 is specifically as follows: constructing a wavelet transform model, using the wavelet transform model to decompose the inertia response capability matrix and the key mode matrix at different time scales, extracting features from the decomposition results, and obtaining a multi-scale feature matrix; constructing a multi-level objective function, and assigning weights to the multi-level objective function based on the multi-scale feature matrix and pre-stored system constraints; obtaining an optimized objective function set based on the multi-level objective function; using a distributed optimization algorithm to solve the multi-objective optimization problem between the optimized objective function set and the pre-stored system operation boundary, and obtaining an optimized control sequence and a response strategy matrix.

[0015] According to one aspect of the present application, step S5 is specifically: obtaining an optimized control sequence and a response strategy matrix, constructing an adaptive control law, obtaining a control instruction according to the adaptive control law based on the optimized control sequence and the response strategy matrix, and generating a control instruction matrix for the system; obtaining the control instruction matrix and a pre-stored system state vector, constructing a Lyapunov function, and using the Lyapunov function to analyze the system stability according to the control instruction matrix and the system state vector to obtain a stability evaluation index; obtaining the control instruction matrix, the stability evaluation index and the pre-stored real-time system state, constructing a control instruction correction function, inputting the control instruction matrix, the stability evaluation index and the real-time system state into the control instruction correction function to perform real-time correction on the control instruction, and obtaining a corrected control instruction sequence and execution feedback data.

[0016] Beneficial effects of the present invention:

[0017] The present invention realizes the intelligent control of the port's multi-energy supply system participating in the grid-load interaction by constructing a complete data processing, analysis and control closed loop. First, through multi-dimensional data preprocessing, a standardized system characteristic matrix and dynamic response capability indicators are obtained; second, through the construction and analysis of the time-varying grid-load interaction diagram, the dynamic characteristics of the system are accurately characterized; third, by constructing a time series prediction model, the future response state of the system is accurately predicted; finally, through multi-level optimization algorithms and real-time correction mechanisms, the reliable execution of control instructions is ensured. The entire solution takes into account the characteristics of the port's multi-energy supply system, such as large load fluctuations and high response requirements. Through innovative algorithm design and optimization strategies, the system's grid-load interaction capabilities are significantly improved, and efficient and reliable inertia response control is achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 is a flow chart of a port multi-energy supply system control method based on automatic inertia response according to an embodiment of the present invention;

[0019] Figure 2 is a flowchart of step S1 of one embodiment of the present invention;

[0020] Figure 3 is a flow chart of step S2 of one embodiment of the present invention;

[0021] Figure 4 is a flow chart of step S3 of one embodiment of the present invention;

[0022] Figure 5 is a flow chart of step S4 of one embodiment of the present invention;

[0023] Figure 6 FIG. 4 is a flow chart of step S5 according to an embodiment of the present invention. DETAILED DESCRIPTION

[0024] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0025] like Figure 1 As shown, an embodiment of the present invention provides a port multi-energy supply system control method based on automatic inertia response, comprising the following steps:

[0026] S1. Collect time series data of the port's multi-energy supply system, construct a synchronized time series tensor based on the time series data, extract dynamic characteristics based on the synchronized time series tensor, and obtain a standardized system characteristic matrix and dynamic response capability indicators. The time series data includes real-time system power data, port equipment operating status data, historical operating data, and feedback data.

[0027] S2. Construct a time-varying grid-load interaction diagram based on the standardized system characteristic matrix, dynamic response capability indicators, pre-stored grid operating parameters, and load fluctuation data. Extract features from the time-varying grid-load interaction diagram to generate a spatiotemporal feature matrix. Use the spatiotemporal feature matrix for pattern recognition to obtain grid-load interaction feature vectors and key pattern matrices. Fusion feature vectors and feature mapping matrices are obtained through multi-level fusion.

[0028] S3. Build a response characteristic tensor based on the dynamic response capability index, the grid-load interaction eigenvector, and the fusion eigenvector. Build a time series prediction model based on the response characteristic tensor, pre-stored system operation constraints, and historical response data to obtain an inertia response capability matrix.

[0029] S4. Based on the inertia response capability matrix and the key mode matrix, a multi-scale feature matrix is ​​obtained by wavelet transform and optimized to obtain the optimized control sequence and response strategy matrix;

[0030] S5. Generate a control instruction matrix of the system based on the optimized control sequence and the response strategy matrix, and perform stability analysis and real-time correction on the control instruction matrix to obtain a corrected control instruction sequence and execution feedback data.

[0031] By building a complete data processing, analysis, and control closed loop, intelligent control of the port's multi-energy supply system participating in grid-load interaction is achieved. First, through multi-dimensional data preprocessing, a standardized system characteristic matrix and dynamic response capability indicators are obtained; second, through the construction and analysis of time-varying grid-load interaction diagrams, the dynamic characteristics of the system are accurately characterized; third, by building a time series prediction model, the future response state of the system is accurately predicted; finally, through multi-level optimization algorithms and real-time correction mechanisms, the reliable execution of control instructions is ensured. The entire solution takes into account the characteristics of the port's multi-energy supply system, such as large load fluctuations and high response requirements. Through innovative algorithm design and optimization strategies, the system's grid-load interaction capabilities are significantly improved, and efficient and reliable inertia response control is achieved.

[0032] like Figure 2 As shown, according to one aspect of the present application, step S1 is specifically:

[0033] S11. Collect time series data from the port's multi-energy supply system and construct a synchronized time series tensor based on the time series data. Specifically, the collected data can be timestamp aligned and serialized; the aligned data can be constructed into a four-dimensional tensor; and time series interpolation and resampling can be performed on the four-dimensional tensor to obtain the synchronized time series tensor.

[0034] S12. Calculate a data quality assessment matrix based on the synchronized time series tensor. Specifically, a comprehensive data quality score can be obtained based on the calculation integrity index, consistency index, and timeliness index of the synchronized time series tensor, and a data quality assessment matrix can be generated based on the comprehensive data quality score.

[0035] S13. Use the data quality assessment matrix to perform anomaly detection and correction on the synchronized time series tensor to obtain a cleaned data tensor.

[0036] Specifically, the steps of abnormality detection and correction may include steps S131 to S135:

[0037] S131. Obtain the synchronized time series tensor and data quality assessment matrix, calculate the first-order difference and second-order difference values ​​of the time series, and generate a time series feature vector; use the time series feature vector to perform sliding window statistical operations and generate a time domain anomaly score.

[0038] Specifically, the time domain anomaly score generation step includes steps S1311 to S1313:

[0039] S1311. Obtain a synchronized time series tensor, calculate the first-order difference value and the second-order difference value of adjacent time points, and generate a differential feature matrix.

[0040] S1312. Use the differential feature matrix to construct a sliding window, calculate the mean, variance and skewness within the window, and generate a statistical feature vector.

[0041] S1313. Calculate the Z-score based on the statistical feature vector to generate a time domain anomaly score.

[0042] The calculation formula of the time domain anomaly score TIME_S(x) is:

[0043] TIME_S(x)=DBSCAN_S(x)+TS_features(x)

[0044] Among them, DBSCAN_S(x) is the density clustering score, and TS_features(x) includes time domain features such as trend, periodicity, and mutation. Specifically, the calculation formula of DBSCAN_S(x) is:

[0045] DBSCAN_S(x)=exp(-ρ(x) / ε)·(1-μ(x))

[0046] Among them, ρ(x) is the local density of sample x, ε is the density threshold parameter, and μ(x) is the cluster membership of sample x.

[0047] S132. Perform multi-scale wavelet decomposition on the synchronized time series tensor, calculate the energy distribution characteristics of the decomposition coefficients, and generate a frequency domain anomaly score.

[0048] Specifically, the frequency domain anomaly score generation step may include steps S1321 to S1323:

[0049] S1321. Obtain the synchronized time series tensor, perform a six-layer decomposition operation of the db4 wavelet basis, obtain wavelet coefficients at each scale, and generate a wavelet coefficient matrix.

[0050] S1322. Calculate the energy distribution of each frequency band in the wavelet coefficient matrix to generate an energy eigenvector.

[0051] S1323. Calculate the frequency domain anomaly measurement value using the energy eigenvector to generate a frequency domain anomaly score.

[0052] Among them, the calculation formula of frequency domain anomaly score FREQ_S(x) is:

[0053] FREQ_S(x)=W_S(x)+FFT_features(x)

[0054] Among them, W_S(x) is the wavelet coefficient anomaly score, and FFT_features(x) includes frequency domain features such as spectrum characteristics and harmonic analysis.

[0055] Specifically, the calculation formula of W_S(x) is:

[0056] W_S(x)=Σ(ψj·|dj(x)|)

[0057] Among them, ψj is the j-th layer wavelet coefficient weight, dj(x) is the j-th layer wavelet coefficient of sample x.

[0058] S133. Calculate the local density distribution and distance distribution of sample points based on the synchronized time series tensor to generate a spatial anomaly score. Specifically, the DBSCAN algorithm can be used to calculate the density distribution of sample points and generate a density feature matrix. Based on the density feature matrix, the core point set and the boundary point set are calculated to generate a clustering result vector. The local outlier factor (LOF) value is calculated for each type of sample in the clustering result vector to generate a LOF feature vector. The density feature matrix and the LOF feature vector are weighted and fused to generate a spatial anomaly score. The calculation formula for the spatial anomaly score SPATIAL_S(x) is:

[0059] SPATIAL_S(x)=LOF_S(x)+SPATIAL_corr(x)

[0060] Among them, LOF_S(x) is the local outlier factor score, and SPATIAL_corr(x) represents the spatial correlation feature. Specifically, the calculation formula of LOF_S(x) is:

[0061] LOF_S(x)=(Σk_dist(x,Nk(x))) / (k·avgDist(Nk(x)))

[0062] Among them, k_dist(x,Nk(x)) is the distance from sample x to its k nearest neighbors, and avgDist(Nk(x)) is the average distance between k nearest neighbor samples.

[0063] S134: Input the time domain anomaly score, frequency domain anomaly score, and spatial domain anomaly score into the evidence theory fusion framework, calculate the comprehensive confidence, and generate the anomaly comprehensive score. Specifically, the basic probability distribution value of each score is calculated to generate the evidence weight matrix; the Dempster combination rule operation is performed based on the evidence weight matrix to generate the fusion confidence vector; the fusion confidence vector is used to calculate the final anomaly degree to generate the anomaly comprehensive score.

[0064] In the fusion process, an anomaly detection fusion method can be used, specifically: the anomaly detection result, i.e., the anomaly score AS(x), is calculated based on the following formula:

[0065] AS(x)=α·TIME_S(x)+β·FREQ_S(x)+γ·SPATIAL_S(x)

[0066] Among them, TIME_S(x) is the time domain anomaly score, FREQ_S(x) is the frequency domain anomaly score, SPATIAL_S(x) is the spatial domain anomaly score, α, β, and γ are the weights of the time domain anomaly score, frequency domain anomaly score, and spatial domain anomaly score, respectively.

[0067] And α+β+γ=1.

[0068] By combining the dual anomaly detection mechanism of DBSCAN density clustering and local outlier factor (LOF), the accuracy and robustness of anomaly detection are significantly improved. Specifically, the DBSCAN algorithm first identifies core samples, boundary samples, and noise points based on density, providing a preliminary classification for different types of data points; then, the LOF score is calculated for each type of sample. This hierarchical calculation strategy avoids the limitations of the traditional LOF algorithm on data with uneven density. By calculating differential features in the time domain, performing wavelet decomposition in the frequency domain, and analyzing the local density distribution in the spatial domain, a comprehensive anomaly detection framework is constructed. Finally, the detection results of the three dimensions are fused through the Dempster-Shafer evidence theory, which not only takes into account the reliability of each detector, but also improves the robustness of the system through the evidence conflict processing mechanism. This multi-dimensional anomaly detection method can effectively identify and process various types of abnormal patterns in the port's multi-energy supply system, including instantaneous fluctuations, gradual anomalies, and structural anomalies.

[0069] S135. Correct the abnormal data in the synchronized time series tensor using the comprehensive anomaly score to generate a cleaned data tensor. Specifically, abnormal data points in the synchronized time series tensor can be identified based on the comprehensive anomaly score; a data reconstruction model can be constructed using a local linear embedding method; the abnormal data points can be corrected using the data reconstruction model, and the cleaned data tensor can be output.

[0070] S14. Use the sliding window method to extract dynamic characteristics of the cleaned data tensor to obtain the standardized system characteristic matrix and dynamic response capability index.

[0071] Among them, the calculation formula of feature importance score FI(f) is:

[0072] FI(f)=η·GainLB(f)+θ·StabLB(f)

[0073] Where: GainLB(f)=Σ(gi·splitf,i) / N; GainLB(f) is the gain value of feature f in the LightGBM model, gi is the gain value of the i-th split, splitf,i indicates the number of times feature f is used in the i-th split, and N is the total number of splits; StabLB(f)=1-σ(gainf) / μ(gainf); StabLB(f) is the stability score of feature f, σ(gainf) is the standard deviation of the gain of feature f, and μ(gainf) is the average gain of feature f; η and θ are weight coefficients and η+θ=1; f is the input feature.

[0074] In a further embodiment, the temporal correlation score may be considered when calculating the feature importance score FI(f). In this case, the calculation formula of the feature importance score FI(f) is:

[0075] FI(f)=η·GainLB(f)+θ·StabLB(f)+λ·TempCorr(f)

[0076] Among them, TempCorr(f) represents the temporal correlation score of the feature, η, θ, λ are weight coefficients, and η+θ+λ=1.

[0077] According to one aspect of the present application, step S14 may specifically be:

[0078] S141. Obtain a preset time window parameter set; calculate the optimal window width based on an adaptive algorithm; construct multiple overlapping sliding windows, and output a multi-scale window sequence and a window parameter matrix.

[0079] S142. Obtain the cleaned data tensor and multi-scale window sequence, calculate the statistical features within the window, extract the time-frequency features within the window, and construct a feature candidate set; perform feature correlation analysis on the feature candidate set to obtain the original feature matrix.

[0080] S143. Obtain the original feature matrix; construct a feature evaluation model based on LightGBM (Light Gradient Boosting Machine) to calculate the importance score of the original features to the target variable (the result of the system's original control instructions participating in the inertia response in historical operating data); perform feature stability analysis to obtain the feature importance matrix.

[0081] S144. Obtain the original feature matrix and the feature importance matrix; construct a feature weighting model; calculate the dynamic response capability index, perform response capability standardization processing, and obtain the standardized system characteristic matrix and the dynamic response capability index.

[0082] By adopting a multi-scale window sequence processing solution and dynamically calculating the optimal window width with an adaptive algorithm, accurate capture of features at different time scales is achieved. In particular, by simultaneously calculating statistical features and time-frequency features within multiple sliding windows and combining the LightGBM feature evaluation model to calculate feature importance scores, it is not only possible to identify the key features that have the greatest impact on the system's responsiveness, but also to ensure the reliability of the selected features through feature stability analysis. This step calculates the dynamic response capability index through a feature weighted model, enabling the system to accurately evaluate the response characteristics of the port's multi-energy supply system under different operating conditions, providing reliable feature support for subsequent optimization control.

[0083] In step S1 of this embodiment, by constructing an innovative method of synchronized time series tensors, real-time power data, equipment operating status data, historical operating data, and feedback data are uniformly integrated into a four-dimensional tensor space, achieving spatiotemporal alignment of heterogeneous data sources. By introducing a data quality assessment matrix, a comprehensive quality assessment of the data is performed, including the three dimensions of integrity, consistency, and timeliness, so that data quality issues can be accurately quantified and located. A triple anomaly detection mechanism in the time domain, frequency domain, and spatial domain is adopted, which can not only identify sudden anomalies, but also discover slowly accumulated anomaly patterns, effectively reducing the false alarm rate. Finally, through the dynamic feature extraction method of the sliding window, the instantaneous characteristics of the system are captured while retaining historical evolution information, providing a reliable feature basis for subsequent response capability evaluation. This multi-dimensional, multi-level data preprocessing solution significantly improves data quality, allowing subsequent analysis and decision-making to be based on high-quality data.

[0084] like Figure 3 As shown, according to one aspect of the present application, step S2 is specifically:

[0085] S21. Obtain a standardized system characteristic matrix, pre-stored grid operating parameters, pre-stored load fluctuation data, and a dynamic response capability index to construct a time-varying grid-load interaction diagram. Specifically, a time-varying weight function can be constructed by obtaining the standardized system characteristic matrix, the dynamic response capability index, the pre-stored grid operating parameters, and the load fluctuation data. Edge weights between network nodes are calculated based on the time-varying weight function, and the time-varying grid-load interaction diagram is constructed using the edge weights.

[0086] Among them, the core features of the time-varying network-load interaction graph include: (1) node type and attributes: source node: represents energy output points such as power generation units and energy storage devices; load node: represents energy consumption points such as load users and electrical equipment; each node contains dynamic characteristic parameters (such as power, voltage, etc.); (2) edge attributes: weight: represents the interaction intensity between nodes; direction: represents the direction of energy flow; time-varying characteristics: the weight changes dynamically over time; (3) time-varying characteristic expression: edge weight function: w(t)=f(P(t),V(t),I(t)); where P(t), V(t), and I(t) are power, voltage, and current at time t, respectively; node state function: s(t)=g(load characteristics, response capability); (4) interaction intensity calculation: instantaneous interaction intensity: IS(t)=w(t)×s(t); cumulative interaction intensity: IC=∫IS(t)dt.

[0087] Through the time-varying grid-load interaction diagram, the following key information can be analyzed: the power transmission path and intensity between the power source and the load; the dynamic coupling relationship between the nodes of the system; the timing characteristics and patterns of energy flow; the system response capability and stability performance, etc.

[0088] S22. Construct a spatiotemporal attention model based on the time-varying network load interaction graph and the interaction intensity matrix; use the spatiotemporal attention model to extract features from the time-varying network load interaction graph and generate a spatiotemporal feature matrix and an attention weight matrix.

[0089] According to one aspect of the present application, step S22 may specifically include steps S221 to S224:

[0090] S221. Obtain a time-varying network load interaction graph; construct a multi-head attention calculation unit; calculate the query matrix, key matrix, and value matrix; perform a scaled dot product operation to obtain the original attention score; apply the Softmax function to obtain the normalized attention weight and obtain the initial attention matrix.

[0091] S222. Obtain the time-varying network load interaction graph and the initial attention matrix; construct a graph convolution layer to calculate the first-order neighborhood features of the nodes in the time-varying network load interaction graph; aggregate the neighborhood information to obtain the node embedding vector; apply the activation function processing and output the graph feature matrix.

[0092] S223. Construct a bidirectional gated recurrent unit, calculate the states of the input gate, forget gate, and output gate based on the graph feature matrix, update the hidden layer state to obtain a time series feature sequence, perform sequence normalization processing, and obtain standardized time series features.

[0093] S224. Construct a feature fusion network, calculate the weight coefficients of static features and dynamic features according to the feature fusion network based on the graph feature matrix and the standardized temporal features, perform weighted fusion operations, and obtain the spatiotemporal feature matrix and the attention weight matrix.

[0094] Among them, the spatiotemporal attention model can be specifically:

[0095] The attention score A(Q,K,V)=softmax(Q·K^T / √dk)·V·exp(-λt·Δt); where: Q=WqX is the query matrix, K=WkX is the key matrix, V=WvX is the value matrix, X is the input feature, Wq, Wk, Wv are the corresponding weight matrices; dk is the feature dimension; λt is the time decay coefficient; Δt is the time interval matrix; exp(-λt·Δt) is the time decay factor.

[0096] Preferably, in step S22, a time decay factor can also be introduced on the basis of the multi-head attention calculation unit to solve the long-term dependency problem in long sequence modeling. Specifically, the step of generating a spatiotemporal feature matrix and an attention weight matrix (step S22) specifically includes: constructing a multi-head attention calculation unit containing a time decay factor based on the time-varying network load interaction graph and the interaction intensity matrix to generate an attenuated attention matrix; performing a graph convolution operation using the attenuated attention matrix to extract local features of the nodes and generate a graph feature matrix; inputting the graph feature matrix into a bidirectional gated recurrent unit to extract temporal dependency features and generate a temporal feature vector; performing a feature fusion operation on the graph feature matrix and the temporal feature vector to generate a spatiotemporal feature matrix and an attention weight matrix. Specifically, it may include the following steps S221'~S224':

[0097] S221', based on the time-varying network load interaction graph, calculate the time decay factor in the form of exponential decay and generate a time decay vector; construct a multi-head attention calculation unit, calculate the query matrix, key matrix and value matrix, and generate an attention basic matrix; perform element-wise product operation on the time decay vector and the attention basic matrix to generate a weighted attention score; perform Softmax normalization operation on the weighted attention score to generate a decaying attention matrix.

[0098] S222', based on the time-varying network load interaction graph and the attenuated attention matrix, calculate the first-order neighborhood node set of each node in the graph and generate a neighborhood index matrix; use the neighborhood index matrix to perform spatial attention convolution operations, extract the local topological features of the nodes, and generate a local feature matrix; input the local feature matrix into the residual connection module, fuse the original features and the convolution features, and generate a graph feature matrix.

[0099] S223', obtain the graph feature matrix, construct a bidirectional recurrent neural network with a gating mechanism, calculate the forward propagation sequence and the backpropagation sequence, and generate a bidirectional feature sequence; use the bidirectional feature sequence to construct a temporal memory module, calculate the long-term dependency features and the short-term dependency features, and generate a memory feature matrix; perform sequence normalization on the memory feature matrix to generate a temporal feature vector.

[0100] S224', obtain the graph feature matrix and the time series feature vector, build a feature fusion network, calculate the importance weights of static features and dynamic features, and generate a feature weight vector; use the feature weight vector to perform weighted fusion operations to generate a spatiotemporal feature matrix; calculate multi-level attention scores based on the spatiotemporal feature matrix to generate an attention weight matrix.

[0101] This method, which combines graph structure and sequence processing, enables the system to simultaneously capture the spatial structural characteristics and temporal evolution characteristics of grid-load interaction. It is particularly suitable for processing complex grid-load interaction patterns in port multi-energy supply systems, and improves the system's ability to model long-term dependencies.

[0102] S23. Obtain the spatiotemporal feature matrix and the attention weight matrix, and construct a dynamic density clustering model; use the dynamic density clustering model to perform cluster analysis on the spatiotemporal feature matrix, track the time evolution characteristics of the clustering results, and obtain the network-load interaction feature vector, key pattern matrix, and pattern evolution matrix.

[0103] S24. Based on the grid-load interaction feature vector, key pattern matrix and pattern evolution matrix, a three-layer feature fusion network is constructed and multi-level fusion is performed to obtain the fused feature vector and feature mapping matrix.

[0104] According to one aspect of the present application, the multi-level fusion step specifically includes steps S241 to S244:

[0105] S241. Obtain the network-load interaction feature vector, build an autoencoder network, perform feature dimensionality reduction to obtain compressed features, reconstruct features to obtain reconstructed features, calculate the reconstruction error, and output the underlying fusion features;

[0106] S242, obtain the key pattern matrix and the underlying fusion features, construct a dynamic convolution kernel, perform multi-scale convolution operations, obtain the intermediate feature map through pooling, extract the salient features, and output the middle-level fusion features;

[0107] S243. Obtain the pattern evolution matrix and the middle-level fusion features, construct an evolution feature extractor, calculate the state transition probability of the time series, extract the evolution pattern features, fuse the static and dynamic features, and output the high-level fusion features;

[0108] S244. Obtain bottom-layer fusion features, middle-layer fusion features, and high-layer fusion features, construct a multi-layer optimization objective function, calculate the loss function value of each layer feature, perform gradient descent optimization, update the fusion weights, and output the fusion feature vector and feature mapping matrix.

[0109] Among them, the three-layer feature fusion method can be specifically as follows: fusion feature F=WH·(WM·(WL·FL+bL)+bM)+bH; where: FL=AE(X) is the bottom-level feature extracted by the autoencoder, WL is the bottom-level weight matrix, bL is the bottom-level bias; WM is the middle-level weight matrix, bM is the middle-level bias; WH is the high-level weight matrix, bH is the high-level bias; the autoencoder function AE(X)=Dec(Enc(X)), where X is the input feature, Enc(X) represents the encoding output by the encoder function, and Dec(Enc(X)) represents the reconstructed output of the encoding by the decoder function.

[0110] By constructing a three-layer feature fusion network, the team achieved effective fusion of different types of features at the bottom, middle, and top layers. At the bottom layer, an autoencoder network performs feature dimensionality reduction and reconstruction, optimizing feature extraction quality through feedback from reconstruction errors. At the middle layer, dynamic convolution kernels perform multi-scale convolution operations to extract mesoscale pattern features. At the top layer, an evolving feature extractor captures state transitions in time series. This hierarchical feature fusion strategy not only achieves an organic combination of static and dynamic features, but also ensures the accuracy and stability of the fusion results through gradient descent optimization of the multi-layer optimization objective function.

[0111] In step S2 of this embodiment, the dynamic capture of network-load interaction characteristics is achieved through the innovative method of constructing a time-varying network-load interaction graph. Specifically, the edge weights between network nodes are calculated through a time-varying weight function to construct an interaction intensity matrix that reflects the dynamic characteristics of the system; the spatiotemporal attention model is used to extract features from the time-varying network-load interaction graph, and the spatiotemporal feature matrix is ​​obtained through temporal dependency modeling; the spatiotemporal feature matrix is ​​clustered and analyzed using a dynamic density clustering model to achieve accurate identification and tracking of network-load interaction patterns. In particular, in the feature fusion network, through the step-by-step fusion of the bottom, middle and high layers, a multi-level integration of original features, pattern features and evolutionary features is achieved, which not only improves the richness of feature expression, but also ensures the accuracy of the fusion results by optimizing the objective function, providing reliable feature support for the response capability evaluation of the port multi-energy supply system.

[0112] like Figure 4 As shown, according to one aspect of the present application, step S3 is specifically:

[0113] S31. Obtain dynamic response capability indicators, grid-load interaction characteristic vectors, and fusion characteristic vectors, build a nonlinear response model, calculate response characteristic parameters using the nonlinear response model, and build a response characteristic tensor based on the response characteristic parameters.

[0114] The nonlinear response model can be specifically: response output R(xr)=K·xr·(1-exp(-tr / τ))·φ(xr); where: φ(xr)=0,

[0115] When |xr|≤δ; φ(xr)=sign(xr)·(|xr|-δ), when δ<|xr|≤Δ; φ(xr)=sign(xr)·(Δ-δ),

[0116] When |xr|>Δ; K is the gain coefficient, τ is the time constant, tr is the response time, δ is the dead zone threshold, Δ is the saturation threshold; xr is the input deviation.

[0117] According to one aspect of the present application, step S31 specifically includes the following steps S311 to S314:

[0118] S311. Obtain the dynamic response capability index and the grid-load interaction characteristic vector; construct a recursive neural network structure; set the input layer neuron weights, calculate the hidden layer state transfer matrix according to the recursive neural network structure based on the dynamic response capability index and the grid-load interaction characteristic vector, generate the initial model parameters of the nonlinear response model, and construct the nonlinear response basis function.

[0119] S312, calculating the response amplitude attenuation coefficient based on the fused eigenvector and the initial model parameters; performing phase shift calculation to generate frequency modulation parameters, and fusing to obtain the initial characteristic parameter set and the response parameter matrix.

[0120] S313, obtaining nonlinear response basis functions and response parameter matrices; constructing a three-dimensional tensor framework, calculating the time dimension feature map of the response parameter matrix, performing frequency dimension transformation, generating a power dimension map, and outputting a response characteristic tensor.

[0121] S314. Obtain an initial characteristic parameter set and a response characteristic tensor; construct a parameter optimization objective function; calculate the gradient update direction according to the parameter optimization objective function based on the initial characteristic parameter set and the response characteristic tensor, perform Adam optimization iteration, update the model parameters, and obtain a nonlinear parameter set.

[0122] By introducing a dead zone characteristic description function into the nonlinear response model, the system's start-stop and hysteresis characteristics are accurately characterized. Specifically, a recursive neural network structure is used to calculate the hidden layer state transition matrix, combined with a bidirectional gating unit to process long-term dependencies. Then, a piecewise function model is used to accurately describe the boundary characteristics of the dead zone interval. This step establishes a complete dead zone characteristic matrix by calculating the forward dead zone threshold and the reverse dead zone threshold. A three-dimensional characteristic mapping framework is then used to achieve unified mapping of the time, frequency, and power dimensions. Finally, through iterative optimization of the parameter optimization objective function, the system's modeling accuracy for nonlinear characteristics is significantly improved.

[0123] In another aspect of the present application, a dead zone characteristic description function can be introduced into the nonlinear response model to accurately characterize the system's start-stop and hysteresis characteristics. A recursive neural network structure is used to calculate the hidden layer state transition matrix, combined with a bidirectional gating unit to process long-term dependencies, and then a piecewise function model is used to accurately describe the boundary characteristics of the dead zone interval. In this case, the step of constructing a response characteristic tensor (step S31) may specifically include: obtaining a dynamic response capability indicator and a grid-load interaction feature vector, calculating the hidden layer state transition matrix and initial neuron weights based on the recursive neural network structure, and generating a basic model parameter set; constructing a piecewise function model using the basic model parameter set and pre-stored system operation history data, calculating the forward dead zone threshold and reverse dead zone threshold, and generating a dead zone characteristic matrix; calculating the response amplitude attenuation coefficient and phase offset parameter based on the basic model parameter set and the dead zone characteristic matrix according to the nonlinear response function to generate an initial characteristic parameter set; constructing a three-dimensional characteristic mapping framework based on the initial characteristic parameter set, calculating the time dimension feature mapping value, the frequency dimension transformation coefficient, and the power dimension mapping parameter, and generating a response characteristic tensor. Specifically, step S31 may include the following steps S311'-S314':

[0124] S311', obtain dynamic response capability indicators and grid-load interaction feature vectors, construct a recursive neural network with bidirectional gating units, calculate the weight matrix from the input layer to the hidden layer and the weight matrix from the hidden layer to the output layer, and generate a neural network weight set; use the neural network weight set to calculate the hidden layer state transition probability matrix and generate a state transition feature vector; input the state transition feature vector into the parameter initialization module, calculate the initial neuron threshold and activation function parameters, and generate a basic model parameter set.

[0125] S312', obtain the basic model parameter set and pre-stored system operation history data, calculate the forward startup threshold and reverse startup threshold of the system response, and generate a startup characteristic vector; use the startup characteristic vector to construct a piecewise function model, calculate the upper and lower boundary values ​​of the dead zone interval, and generate a boundary parameter matrix; calculate the hysteresis characteristic parameters of the dead zone interval based on the boundary parameter matrix, and generate a dead zone characteristic matrix.

[0126] S313', based on the basic model parameter set and the dead zone characteristic matrix, calculate the amplitude attenuation characteristics and phase shift characteristics of the system response to generate a dynamic characteristic vector; use the dynamic characteristic vector to construct a nonlinear response function, calculate the gain coefficient and time constant of the system, and generate an initial characteristic parameter set.

[0127] S314', obtain the initial characteristic parameter set, construct a three-dimensional characteristic mapping framework, calculate the mapping coefficients of the initial characteristic parameter set in the time dimension, frequency dimension and power dimension respectively, and generate a dimensional mapping matrix; use the dimensional mapping matrix and the dead zone characteristic matrix to construct a characteristic tensor, calculate the dimensional components of the tensor, and generate a response characteristic tensor; construct a parameter optimization objective function based on the response characteristic tensor, perform parameter optimization calculation based on the Adam algorithm, and generate a nonlinear parameter set and a response model parameter matrix.

[0128] In this embodiment, a bidirectional gating unit is introduced to improve the ability to capture long-term dependencies; dead-zone characteristic modeling is added to better align with actual system characteristics; and a piecewise function model is adopted to improve the accuracy of expressing nonlinear characteristics. This more accurately characterizes response characteristics, including: describing the system startup threshold through a dead-zone characteristic matrix; and considering hysteresis characteristics to more closely match actual system behavior. The piecewise function model simplifies computational complexity; and the Adam algorithm accelerates parameter convergence. In summary, step S31 in this embodiment takes into account the key characteristics of the actual system, making parameter optimization more comprehensive, providing better adaptability, and enhancing engineering practicality.

[0129] S32. Construct a multi-dimensional evaluation function based on the response characteristic tensor and the pre-stored system operation constraints, use the multi-dimensional evaluation function to calculate the sub-item indicators, perform weighted calculation on the sub-item indicators, and output a comprehensive evaluation indicator matrix.

[0130] S33. Build a time series prediction model using the comprehensive evaluation index matrix and pre-stored historical response data, perform time series prediction based on the time series prediction model, and output an inertia response capability matrix.

[0131] In step S3 of this embodiment, a nonlinear response model and a comprehensive evaluation index system are constructed to accurately assess the system's responsiveness. First, the nonlinear response model is used to calculate response characteristic parameters, and the response characteristic tensor is used to comprehensively describe the system's dynamic response characteristics. Then, a comprehensive evaluation index matrix is ​​constructed based on a multi-dimensional evaluation function, enabling a quantitative assessment of the system's responsiveness. Finally, a time series prediction model is constructed that not only accurately predicts the system's future response state but also captures key time series features through an attention mechanism. This method, combining nonlinear modeling and time series prediction, is particularly well-suited for handling the complex response characteristics of port multi-energy supply systems and provides accurate inertia response capability assessment results.

[0132] In another embodiment of the present application, in step S31, a nonlinear parameter set may be constructed while constructing the response characteristic tensor based on the response characteristic parameters. In this case, steps S32 and S33 may be:

[0133] S32', construct a multi-dimensional evaluation function including dead zone characteristics based on the response characteristic tensor, nonlinear parameter set and pre-stored system operation constraints, calculate the time characteristic index, amplitude characteristic index and frequency characteristic index of the system response, and generate a sub-item indicator matrix; use the sub-item indicator matrix and nonlinear parameter set to perform weighted fusion operation to generate a comprehensive evaluation indicator matrix.

[0134] S33', obtain the comprehensive evaluation index matrix, the nonlinear parameter set and the pre-stored historical response data, build a time series attention network, use the nonlinear parameter set as the modulation parameter of the network, and generate the prediction model parameters; use the prediction model parameters to build a time series prediction model, perform multi-step prediction operations, and generate an inertia response capability matrix.

[0135] like Figure 5 As shown, according to one aspect of the present application, step S4 may specifically be:

[0136] S41. Construct a wavelet transform model based on the inertia response capability matrix and the key mode matrix, and perform decomposition and feature extraction of different regimes through wavelet transform to obtain a multi-scale feature matrix. Specifically, this may include obtaining the inertia response capability matrix and the key mode matrix, constructing an improved wavelet transform model; utilizing the wavelet transform model to decompose the inertia response capability matrix and the key mode matrix at different time scales; and performing feature extraction on the decomposition results to obtain a multi-scale feature matrix.

[0137] Among them, the improved wavelet transform model is specifically as follows:

[0138] Wψ(a,b)=(1 / √a)·Σx(t)·ψ*((tw-b) / a)·w(tw)

[0139] Where Wψ(a,b) is the wavelet coefficient, ψ((tw-b) / a) is the wavelet basis function, a is the scale parameter, and b is the translation parameter; w(tw)=exp(-ξ·|tw-tc|) is the adaptive weight function, tw represents the current time, tc is the center time, and ξ is the attenuation coefficient; x(t) is the input signal.

[0140] S42: Construct an objective function, and obtain an optimized objective function set based on the objective function based on the multi-scale feature matrix and pre-stored system constraints. Specifically, this may include obtaining the multi-scale feature matrix and pre-stored system constraints; constructing a multi-level objective function; assigning weights to the multi-level objective functions; combining weighted objective functions; and obtaining an optimized objective function set based on the multi-level objective function based on the multi-scale feature matrix and pre-stored system constraints.

[0141] S43. Construct a distributed optimization algorithm; use the distributed optimization algorithm to solve the multi-objective optimization problem, and obtain an optimized control sequence and a response strategy matrix. This may specifically include obtaining a set of optimization objective functions and a pre-stored system operation boundary; constructing a distributed optimization algorithm; and using the distributed optimization algorithm to solve the multi-objective optimization problem between the set of optimization objective functions and the pre-stored system operation boundary, and obtain an optimized control sequence and a response strategy matrix.

[0142] According to one aspect of the present application, step S43 is specifically as follows:

[0143] S431. Obtain the set of optimization objective functions and the pre-stored system operation boundaries; construct a Lagrange multiplier vector; calculate the gradient of the constraint conditions, initialize the original variables and the dual variables, generate an initial solution vector, and output the optimization problem parameter set.

[0144] S432. Obtain the initial solution vector and optimization problem parameter set; construct the ADMM (Alternating Direction Method of Multipliers) algorithm framework; calculate the augmented Lagrangian function value, perform primitive variable update, and output the iterative solution sequence.

[0145] S433. Based on the iterative solution sequence, calculate the variation of adjacent iterative solutions; perform primal residual calculation to generate dual residual; determine the convergence condition; and output a convergence status indicator when the convergence condition is met.

[0146] S434. Based on the iterative solution sequence and the convergence state identifier, the optimal original variable value is extracted, the optimal dual variable value is calculated, the control decision sequence is generated, and the optimized control sequence and response strategy matrix are output.

[0147] A distributed optimization algorithm based on ADMM was employed to solve a large-scale constrained optimization problem. By constructing Lagrange multiplier vectors and calculating constraint gradients, the original problem was decomposed into several easily solvable subproblems. Rapid convergence of the optimization objective was achieved through iterative updates of the primal and dual variables. In particular, a dynamic adjustment strategy was employed when calculating the primal and dual residuals to ensure algorithm stability. The resulting optimized control sequence and response strategy matrix effectively balanced control effectiveness and computational efficiency, meeting the real-time control requirements of the port's multi-energy supply system.

[0148] In step S4 of this embodiment, a multi-scale analysis of the system response characteristics is achieved through an improved wavelet transform method, and the optimal control strategy is obtained through multi-objective optimization. Specifically, by performing a wavelet transform on the inertia response capability matrix, characteristic representations at different time scales are obtained; by constructing a multi-level objective function and assigning weights, a complete set of optimization objective functions is formed; and by using a distributed optimization algorithm to solve the multi-objective optimization problem, the optimal control sequence is obtained while satisfying the system constraints. This method based on multi-scale analysis and multi-objective optimization can not only adapt to the system's response requirements at different time scales, but also improves the solution efficiency through the distributed processing of the optimization algorithm, ensuring the real-time and optimality of the control strategy.

[0149] In another embodiment of the present application, steps S41 and S42 may also be:

[0150] S41', obtain the inertia response capability matrix, the key mode matrix and the characteristic mapping matrix, construct the characteristic transformation function based on the characteristic mapping matrix, calculate the multi-scale representation of the signal, and generate the transformation coefficient matrix; use the transformation coefficient matrix to perform wavelet transform operation to generate a multi-scale characteristic matrix.

[0151] According to one aspect of the present application, step S41′ may specifically include the following steps S411 to S414:

[0152] S411. Calculate signal energy distribution based on the inertia response capability matrix; perform spectrum analysis; construct an optimal wavelet basis selection criterion and output optimized wavelet parameters.

[0153] S412. Obtain the key mode matrix and optimize the wavelet parameters, build a multi-layer wavelet decomposition framework; perform high-frequency and low-frequency separation, calculate the coefficients of each scale, generate a decomposition coefficient matrix, and output multi-scale components.

[0154] S413. Obtain multi-scale components and calculate statistical features of each scale; perform energy concentration analysis, extract scale correlation features, construct feature description vectors, and output scale feature matrix.

[0155] S414. Based on the decomposition coefficient matrix and the scale feature matrix, a signal reconstruction framework is constructed; reconstruction weight coefficients are calculated, an inverse wavelet transform is performed, multi-scale features are fused, and a multi-scale feature matrix is ​​output.

[0156] By constructing an adaptive wavelet transform framework, a multi-scale analysis of inertia response capability is achieved. Specifically, the optimal wavelet basis function is selected by calculating the signal energy distribution, and a multi-layer wavelet decomposition framework is constructed to accurately separate high-frequency and low-frequency components. Simultaneously, characteristic descriptions of different frequency bands are obtained by calculating statistical characteristics at each scale and analyzing energy concentration. This method calculates reconstruction weight coefficients through a signal reconstruction framework, extracting multi-scale features while maintaining signal integrity. This method is particularly suitable for processing power fluctuation characteristics at different time scales in port multi-energy supply systems, providing comprehensive feature support for subsequent optimization control.

[0157] S42', based on the multi-scale feature matrix, the feature mapping matrix and the pre-stored system constraints, construct a feature mapping function, calculate the importance weights of features at different scales, and generate a feature weight vector; based on the feature weight vector, construct a multi-level objective function and generate an optimization objective function set.

[0158] In this embodiment, an adaptive mapping from feature space to target space is achieved by introducing a feature mapping matrix as a core component. This approach offers significant advantages over the first approach: it simplifies the feature extraction process, avoiding the tedious wavelet basis selection and multi-layer decomposition process. Directly constructing the transformation function using the feature mapping matrix ensures both feature extraction accuracy and computational efficiency. In terms of objective function construction, the feature mapping function automatically calculates the importance weights of features at different scales, enabling the objective function to better reflect the actual system characteristics and avoiding the subjectivity that can arise from manually assigning weights. This design not only reduces algorithmic complexity but also enhances the model's adaptability and generalization performance. It is particularly well-suited for handling the rapidly changing dynamic characteristics of power systems and can more efficiently support subsequent optimization control decisions. Under certain operating conditions, the introduction of the feature mapping matrix significantly improves computational efficiency, feature extraction accuracy, and model adaptability: computation time is reduced by an average of over 40%, feature representation accuracy is improved by approximately 25%, and adaptation to system state changes is accelerated by approximately 35%.

[0159] like Figure 6 As shown, according to one aspect of the present application, step S5 is specifically as follows:

[0160] S51. Obtain an optimized control sequence and a response strategy matrix, and construct an adaptive control law; obtain control instructions according to the adaptive control law based on the optimized control sequence and the response strategy matrix, and generate a control instruction matrix.

[0161] S52. Obtain a control instruction matrix and a pre-stored system state vector, and construct a Lyapunov function; use the Lyapunov function to analyze system stability based on the control instruction matrix and the system state vector, and obtain a stability evaluation index.

[0162] S53, obtain the control instruction matrix, stability evaluation index and pre-stored real-time system status, and construct a control instruction correction function; input the control instruction matrix, stability evaluation index and real-time system status into the control instruction correction function to correct the control instruction in real time, and obtain a corrected control instruction sequence and execution feedback data.

[0163] The instruction modification may specifically be:

[0164] Correction instruction u*(t)=u(t)+Kp·e(t)+Ki·∫e(t)dt+Kd·de(t) / dt+δu(t)

[0165] Where: u(t) is the control command matrix at time t, δu(t)=-μ·sign(s(t))·|s(t)|^αu;

[0166] s(t)=e(t)+λu·∫e(t)dt is the sliding surface; e(t) is the tracking error; Kp, Ki, Kd are PID parameters; μ is the convergence rate;

[0167] αu is the convergence index; λu is the sliding surface parameter.

[0168] According to one aspect of the present application, step S53 may specifically include the following steps S531 to S534:

[0169] S531. Calculate the system state deviation based on the control instruction matrix, the stability evaluation index, and the pre-stored real-time system state; perform state prediction, and generate a state evaluation vector.

[0170] S532, obtaining a state evaluation vector; constructing an adaptive correction function; calculating a correction coefficient matrix according to the adaptive correction function based on the state evaluation vector; performing response characteristic analysis, and generating correction function parameters.

[0171] S533. Calculate the corrected control quantity based on the control instruction matrix and the correction function parameters; perform constraint checking and generate a corrected control instruction sequence.

[0172] S534. Based on the modified control instruction sequence and the state evaluation vector, the execution effect evaluation index is calculated and a feedback data packet is generated; the data format of the execution feedback data packet is converted and the execution feedback data is output.

[0173] By constructing an adaptive correction function mechanism, real-time optimization and adjustment of control instructions are achieved. This step first calculates the system state deviation based on the state evaluation vector, evaluates the control effect using a state prediction model, and then uses the correction function parameters to adjust the control instructions in real time to ensure the execution effect of the control instructions. This closed-loop correction mechanism not only improves control accuracy but also ensures the reliability and stability of the control system through the real-time calculation of the execution effect evaluation index.

[0174] In step S5, the reliable execution of control instructions is achieved by constructing an adaptive control law and a real-time correction mechanism. First, a control instruction matrix is ​​generated based on the optimized control sequence and response strategy matrix. Then, the system stability is analyzed using the Lyapunov function to ensure the stability of the control strategy. Finally, the control instructions are adjusted in real time using the control instruction correction function, achieving a rapid response to changes in the system state. This control method, which combines stability analysis and real-time correction, not only ensures the stability of the control system but also improves control accuracy through a feedback mechanism. It is particularly suitable for scenarios such as port multi-energy supply systems that require high-precision control. By executing real-time collection and analysis of feedback data, a complete control closed loop is formed, significantly improving the control reliability and robustness of the system.

[0175] In some embodiments of the present application, the data processing process is as follows:

[0176] In step S1, the dynamic characteristic data of the port multi-energy supply system is obtained and the data is preprocessed, which specifically includes the following steps:

[0177] S11. Data Collection and Synchronization

[0178] Acquire real-time power data PO(t); equipment operating status data SO(t); historical operating data HO(t); system feedback data D_feedback(t-1); apply timestamp alignment algorithm during data pre-alignment; construct enhanced four-dimensional tensor TO(t)=[PO(t),SO(t),HO(t),D_feedback(t-1)]; perform time series interpolation and resampling on the four-dimensional tensor to obtain the synchronized time series tensor T(t).

[0179] S12. Data Quality Assessment and Enhancement

[0180] Get the synchronized time series tensor T(t); calculate the quality index of T(t); calculate the comprehensive data quality score based on the quality index: , generate the data quality assessment matrix DQM and the data quality enhancement suggestion matrix E. Among them, the quality indicators can include integrity indicators, consistency indicators and timeliness indicators, and the calculation formulas of the three are:

[0181] Completeness index: C(i)=Σ(valid data points) / total data points

[0182] Consistency index: U(i)=1-Σ(data deviation) / reference value

[0183] Timeliness index: T(i)=exp(-λ·Δt)

[0184] S13. Multidimensional anomaly detection and processing

[0185] Acquire the synchronized time series tensor T(t) and the data quality assessment matrix DQM; construct a multidimensional anomaly detection model, perform anomaly detection on the synchronized time series tensor in the time domain, perform wavelet transform anomaly detection on the synchronized time series tensor in the frequency domain, and perform local anomaly factor detection on the synchronized time series tensor in the spatial domain; fuse the above detection results to obtain the anomaly detection result; based on the anomaly detection result, correct the synchronized time series tensor to obtain the cleaned data tensor T'(t).

[0186] S14. Dynamic Feature Extraction and Evaluation

[0187] Get the cleaned data tensor T'(t); construct multiple time windows of different scales and use the time windows to extract features from the cleaned data tensor: ; Calculate the feature importance score: I(f) = XGBoost feature importance score, and calculate the dynamic response capability based on the feature importance score:

[0188] D = Σ(I(f)·W(t,δ)), to obtain the standardized system characteristic matrix M and dynamic response capability index D;

[0189] In step S2, network-load interaction pattern recognition and feature extraction are performed, which specifically includes the following steps:

[0190] S21. Analysis of Enhanced Grid-Load Interaction Model

[0191] The standardized system characteristic matrix M, grid operation parameters E(t), load fluctuation data L(t), and dynamic response capability index D are obtained, and a time-varying weight function w(t)=β0·exp(-γ0t) is constructed, where β0 is the initial weight coefficient and γ0 is the attenuation coefficient. The edge weights E(i,j,t)=w(t)·sim(i,j) between network nodes are calculated based on the time-varying weight function, where sim(i,j) is the cosine similarity of the load curves between nodes i and j. The time-varying grid-load interaction graph G(V,E,t)=F[M,E(t),L(t),D] is constructed using the edge weights. The interaction strength matrix S is calculated based on the time-varying grid-load interaction graph G(V,E,t).

[0192] S22. Deep extraction of spatiotemporal features

[0193] Obtain the time-varying network load interaction graph G(V,E,t) and the interaction strength matrix S; construct the spatiotemporal attention model: A(Q,K,V)=softmax(Q·K^T / √dk)·V·exp(-λt·Δt); and perform feature extraction on the time-varying network load interaction graph G(V,E,t): F(t)=GNN(G(V,E,t),A(t)); perform temporal dependency modeling ST=LSTM(F(t)) on the extracted features to obtain the spatiotemporal feature matrix ST and the attention weight matrix A.

[0194] Attention calculation process: similarity calculation: score=Q×K^T; scaling: scaled_score=score / sqrt(dk); Softmax normalization: weight=softmax(scaled_score); weighted summation: output=weight×V.

[0195] Among them, Q is the query matrix, which represents the target information that currently needs to be paid attention to; in the network-load interaction graph, it represents the characteristics of the node or time point that you want to know; the dimension is [N×dk], where N is the number of nodes and dk is the feature dimension; K represents the information key value that can be queried; in the network-load interaction graph, it represents the feature information of all nodes; the dimension is [N×dk], which needs to match the dimension of Q; V represents the actual information content extracted; in the network-load interaction graph, it contains the actual feature values ​​of the node; the dimension is [N×dv], where dv is the value dimension.

[0196] S23. Adaptive pattern recognition

[0197] Obtain the spatiotemporal feature matrix ST and the attention weight matrix A; construct a dynamic density clustering model based on the spatiotemporal feature matrix and the attention weight matrix. Specifically, the DBSCAN algorithm + time window can be used, and an adaptive density threshold can be used to construct a dynamic density clustering model; use the dynamic density clustering model to perform cluster analysis on the spatiotemporal feature matrix, and track the time evolution characteristics of the clustering results Ev(t)=TrackClusterEvolution(C(t),C(t-1)), where C(t) represents the time series obtained after cluster analysis, and obtain the network-load interaction feature vector I, key pattern matrix K and pattern evolution matrix E.

[0198] S24. Multi-level feature fusion

[0199] Obtain the grid-load interaction feature vector I, key pattern matrix K and pattern evolution matrix E; construct a three-layer feature fusion network based on the grid-load interaction feature vector, key pattern matrix and pattern evolution matrix; fuse the original features at the bottom layer of the feature fusion network; fuse the pattern features at the middle layer of the feature fusion network; and fuse the evolution features at the high layer of the feature fusion network; optimize the fusion results based on the optimization objective function to obtain the fused feature vector R and feature mapping matrix T.

[0200] Among them, the optimization objective function is: ;

[0201] Among them, L_feature (feature layer loss) is used to measure the accuracy of the original feature fusion and calculate the reconstruction error of the underlying network-load interaction feature I; the expression is usually: ||I-I'||², where I' is the reconstructed feature; L_pattern (pattern layer loss) is used to measure the degree of preservation of the key pattern K and measure the structural similarity of the pattern features; it can be expressed as: ||Kf(R)||², where f(·) is the pattern mapping function; L_evolution (evolution layer loss) is used to measure the expression degree of the system evolution characteristic E; it characterizes the dynamic change law of time series; the form is such as: ||Eg(R,t)||², where g(·) is the evolution mapping function; 、 、 Represent the weight coefficients of L_feature, L_pattern, and L_evolution respectively.

[0202] In step S3, automatic inertia response capability evaluation and prediction is performed, specifically including the following steps:

[0203] S31. Analysis of nonlinear inertia response characteristics

[0204] Obtain the dynamic response capability index D, the grid-load interaction characteristic vector I, and the fusion characteristic vector R, and construct a nonlinear response model: R(tr,fr,pr)=F[D,I,φ(tr)], φ(tr)=αr·exp(-βr·tr)·sin(ωr·tr+θr); use the nonlinear response model to calculate the response characteristic parameters and perform parameter adaptive optimization: [αr,βr,ωr,θr]=ArgMin(L_response); based on the response characteristic parameters, construct the response characteristic tensor R(tr,fr,pr) and the nonlinear parameter set Θ. Where R(tr,fr,pr) is the response characteristic tensor, fr is the response frequency, pr is a parameter in the nonlinear response model, and φ(tr) is the time-varying attenuation function of the nonlinear response characteristic:

[0205] αr is the response amplitude coefficient; βr is the attenuation rate coefficient; ωr is the oscillation angular frequency; θr is the initial phase angle; exp(-βr·tr) represents the response attenuation characteristic; sin(ωr·tr+θr) represents the oscillation characteristic.

[0206] S32. Construction of multi-dimensional evaluation indicators

[0207] Obtain the response characteristic tensor R(tr,fr,pr) and the system operation constraint C, design a multidimensional evaluation function E(x)=Σ(wi·xi), where wi is the weight coefficient and xi is the sub-item indicator; use the multidimensional evaluation function to calculate the sub-item indicator xi; perform weighted calculation on the sub-item indicator xi; and output the final result as the comprehensive evaluation indicator matrix A.

[0208] S33. Time series prediction model construction

[0209] Obtain the comprehensive evaluation index matrix A and historical response data HR; construct a temporal attention network M(t), and use the temporal attention network to build a prediction model M(t+τ)=F[M(t),A,HR]; obtain the inertia response capability matrix B;

[0210] In step S4, a multi-time-scale collaborative optimization decision is made, which specifically includes the following steps:

[0211] S41. Time scale decomposition

[0212] Obtain the inertia response capability matrix B and the key mode matrix Key, construct an improved wavelet transform model, and use the wavelet transform model to decompose the inertia response capability matrix B and the key mode matrix Key into different time scales [T1, T2, ..., Tn]. Perform feature extraction on the decomposition results to obtain a multi-scale feature matrix MT.

[0213] S42. Hierarchical optimization target construction

[0214] Construct a multi-level optimization objective function O(y)=Σ(αi·fi(y)), where αi is the weight of each level and y is the parameter to be optimized; assign weights to the multi-level objective functions; combine the weighted objective functions and output the optimization objective function set F.

[0215] S43. Collaborative Optimization Solution

[0216] Obtain the optimization objective function set F and the system operation boundary BB; design a distributed optimization algorithm to solve the multi-objective optimization problem and obtain the optimized control sequence U and the response strategy matrix SS;

[0217] In step S5, the adaptive control strategy is generated and executed, which specifically includes the following steps:

[0218] S51. Control strategy generation

[0219] Obtain the optimized control sequence U and response strategy matrix SS, and construct the adaptive control law CTL(t)=G[U(t), SS(t),e(t)], where e(t) is the tracking error; obtain the control command based on the adaptive control law and obtain the control command matrix CM.

[0220] S52. Stability Analysis

[0221] Obtain the control command matrix CM and the system state vector X, construct the Lyapunov function V(x), use the Lyapunov function to analyze the system stability, and obtain the stability evaluation index.

[0222] S53. Real-time correction and execution

[0223] Obtain the control instruction matrix CM, stability evaluation index Z, and real-time system state Y; design a correction function ΔC(t)=H[CM(t), Z(t), Y(t)] based on real-time feedback; use the control instruction correction function to correct the control instructions in real time to obtain the corrected control instruction sequence C' and execution feedback data D.

[0224] In this embodiment, various real-time data, including power data and equipment operating status, are collected first. A triple-layer data quality control approach is employed: data quality is assessed, anomalies are detected across time, frequency, and space, and dynamic features are extracted using a sliding window. Next, an intelligent grid-load interaction analysis model is constructed. Using a "time-varying grid-load interaction graph," the system's changing characteristics can be dynamically captured. This is like taking a "dynamic X-ray" of the system, allowing not only the current state but also trends to be analyzed. Specifically, a three-level feature fusion network, from bottom to top, is employed to analyze grid-load interaction characteristics, ensuring that no important information is missed. Based on the preceding analysis results, the optimal control strategy is calculated through multi-scale feature analysis and multi-objective optimization. The innovation lies in the use of an improved wavelet transform and distributed optimization algorithm, which both ensures decision quality and improves computational efficiency. Finally, rather than simply executing control commands, a real-time correction mechanism is added. Lyapunov stability analysis and adaptive correction ensure reliable execution of control commands. This is like installing an "automatic error corrector" in the control system, enabling timely adjustment of the control strategy based on actual conditions.

[0225] The technical advantages of the entire solution are primarily reflected in: first, comprehensive data quality control; second, the innovative introduction of a time-varying grid-load interaction diagram and a three-layer feature fusion network; third, consideration of the system's dead zone characteristics, improving control accuracy; and fourth, the use of distributed optimization to enhance computational efficiency. These innovations collectively address key technical challenges in enabling the port's multi-energy supply system to participate in grid regulation, enabling more intelligent and efficient inertia response control.

[0226] The preferred embodiments of the present invention are described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the scope of protection of the present invention.

Claims

1. A control method for a port multi-energy supply system based on automatic inertia response, characterized in that: The following steps are involved: S1. Collect time series data of the port's multi-energy supply system, construct a synchronized time series tensor based on the time series data, extract dynamic characteristics based on the synchronized time series tensor, and obtain a standardized system characteristic matrix and dynamic response capability indicators. The time series data includes real-time system power data, port equipment operating status data, historical operating data, and feedback data. S2. Construct a time-varying grid-load interaction diagram based on the standardized system characteristic matrix, dynamic response capability indicators, pre-stored grid operating parameters, and load fluctuation data. Extract features from the time-varying grid-load interaction diagram to generate a spatiotemporal feature matrix. Use the spatiotemporal feature matrix for pattern recognition to obtain grid-load interaction feature vectors and key pattern matrices. Fusion feature vectors and feature mapping matrices are obtained through multi-level fusion. S3. Build a response characteristic tensor based on the dynamic response capability index, the grid-load interaction eigenvector, and the fusion eigenvector. Build a time series prediction model based on the response characteristic tensor, pre-stored system operation constraints, and historical response data, and perform predictions to obtain an inertia response capability matrix. S4. Based on the inertia response capability matrix and the key mode matrix, a multi-scale feature matrix is ​​obtained by wavelet transform and optimized to obtain the optimized control sequence and response strategy matrix; S5. Generate a control instruction matrix of the system based on the optimized control sequence and the response strategy matrix, and perform stability analysis and real-time correction on the control instruction matrix to obtain a corrected control instruction sequence and execution feedback data.

2. The port multi-energy supply system control method based on automatic inertia response according to claim 1 is characterized in that: Step S1 includes: Collect time series data of the port's multi-energy supply system and construct a synchronized time series tensor based on the time series data; Calculate the data quality assessment matrix based on the synchronized time series tensor; Use the data quality assessment matrix to detect and correct anomalies in the synchronized time series tensor to obtain the cleaned data tensor; The sliding window method is used to extract the dynamic characteristics of the cleaned data tensor to obtain the standardized system characteristic matrix and dynamic response capability index.

3. The port multi-energy supply system control method based on automatic inertia response according to claim 2 is characterized in that: The steps for anomaly detection and correction include: Obtain synchronized time series tensors and data quality assessment matrices, calculate first-order and second-order difference values ​​of the time series, and generate time series feature vectors; perform sliding window statistical operations using the time series feature vectors to generate time domain anomaly scores; Perform multi-scale wavelet decomposition on the synchronized time series tensor, calculate the energy distribution characteristics of the decomposition coefficients, and generate frequency domain anomaly scores; Calculate the local density distribution and distance distribution of sample points based on the synchronized time series tensor to generate a spatial anomaly score; Input the time domain anomaly score, frequency domain anomaly score and spatial domain anomaly score into the evidence theory fusion framework, calculate the comprehensive confidence and generate the anomaly comprehensive score; The abnormal comprehensive score is used to correct the abnormal data in the synchronized time series tensor and generate the cleaned data tensor.

4. The port multi-energy supply system control method based on automatic inertia response according to claim 1 is characterized in that: Step S2 specifically includes: Based on the standardized system characteristic matrix, dynamic response capability indicators, grid operating parameters, and load fluctuation data, a time-varying weight function is constructed. The edge weights between network nodes are calculated based on the time-varying weight function. The edge weights are used to construct a time-varying grid-load interaction graph, and the interaction intensity matrix is ​​calculated based on the time-varying grid-load interaction graph. A spatiotemporal attention model is constructed based on the time-varying network load interaction graph and the interaction intensity matrix. The spatiotemporal attention model is used to extract features from the time-varying network load interaction graph and generate a spatiotemporal feature matrix and an attention weight matrix. Obtain the spatiotemporal feature matrix and attention weight matrix and construct a dynamic density clustering model. Use the dynamic density clustering model to perform cluster analysis on the spatiotemporal feature matrix, track the temporal evolution characteristics of the clustering results, and obtain the network-load interaction feature vector, key pattern matrix, and pattern evolution matrix. Based on the grid-load interaction feature vector, key pattern matrix and pattern evolution matrix, a three-layer feature fusion network is constructed, and multi-level fusion is performed to obtain the fused feature vector and feature mapping matrix.

5. The port multi-energy supply system control method based on automatic inertia response according to claim 4 is characterized in that: The steps to generate the spatiotemporal feature matrix and attention weight matrix include: Based on the time-varying network load interaction graph and the interaction intensity matrix, a multi-head attention calculation unit with a time decay factor is constructed to generate a decaying attention matrix. Use the attenuated attention matrix to perform graph convolution operations, extract local features of nodes, and generate a graph feature matrix; Input the graph feature matrix into the bidirectional gated recurrent unit to extract the temporal dependency features and generate a temporal feature vector; Perform feature fusion operations on the graph feature matrix and the temporal feature vector to generate the spatiotemporal feature matrix and the attention weight matrix.

6. The port multi-energy supply system control method based on automatic inertia response according to claim 4 is characterized in that: The steps of multi-level fusion include: Obtain the network-load interaction feature vector, build an autoencoder network, perform feature dimensionality reduction to obtain compressed features, reconstruct features to obtain reconstructed features, calculate the reconstruction error, and output the underlying fusion features; Obtain the key pattern matrix and underlying fusion features, construct a dynamic convolution kernel, perform multi-scale convolution operations, obtain intermediate feature maps through pooling, extract significant features, and output mid-level fusion features; Obtain the pattern evolution matrix and mid-level fusion features, build an evolution feature extractor, calculate the state transition probability of the time series, extract the evolution pattern features, fuse static and dynamic features, and output high-level fusion features; Obtain the bottom-level fusion features, middle-level fusion features, and high-level fusion features, construct a multi-layer optimization objective function, calculate the loss function value of each layer feature, perform gradient descent optimization, update the fusion weights, and output the fusion feature vector and feature mapping matrix.

7. The port multi-energy supply system control method based on automatic inertia response according to claim 1 is characterized in that: Step S3 includes: Obtain dynamic response capability indicators, grid-load interaction eigenvectors, and fusion eigenvectors, build a nonlinear response model, use the nonlinear response model to calculate response characteristic parameters, and build a response characteristic tensor based on the response characteristic parameters; A multi-dimensional evaluation function is constructed based on the response characteristic tensor and pre-stored system operation constraints. The sub-item indicators are calculated using the multi-dimensional evaluation function, and the sub-item indicators are weighted to output a comprehensive evaluation indicator matrix. A time series prediction model is constructed using a comprehensive evaluation index matrix and pre-stored historical response data. Time series prediction is performed based on the time series prediction model, and an inertia response capability matrix is ​​output.

8. The port multi-energy supply system control method based on automatic inertia response according to claim 7 is characterized in that: The steps to construct the response characteristic tensor include: Obtain dynamic response capability indicators and grid-load interaction feature vectors, calculate the hidden layer state transition matrix and initial neuron weights based on the recursive neural network structure, and generate the basic model parameter set; A piecewise function model is constructed using the basic model parameter set and pre-stored system operation history data to calculate the forward dead zone threshold and the reverse dead zone threshold and generate a dead zone characteristic matrix; Based on the basic model parameter set and the dead zone characteristic matrix, the response amplitude attenuation coefficient and the phase shift parameter are calculated according to the nonlinear response function to generate the initial characteristic parameter set; A three-dimensional characteristic mapping framework is constructed based on the initial characteristic parameter set, and the time dimension feature mapping value, frequency dimension transformation coefficient and power dimension mapping parameter are calculated to generate the response characteristic tensor.

9. The port multi-energy supply system control method based on automatic inertia response according to claim 1 is characterized in that: Step S4 is specifically as follows: Construct a wavelet transform model, use the wavelet transform model to decompose the inertia response capability matrix and the key mode matrix at different time scales, extract features from the decomposition results, and obtain a multi-scale feature matrix; Construct a multi-level objective function, assign weights to the multi-level objective function based on the multi-scale feature matrix and pre-stored system constraints, and obtain an optimized objective function set based on the multi-level objective function; A distributed optimization algorithm is used to solve the multi-objective optimization problem between the optimization objective function set and the pre-stored system operation boundary, and the optimized control sequence and response strategy matrix are obtained.

10. The port multi-energy supply system control method based on automatic inertia response according to claim 1 is characterized in that: Step S5 is specifically as follows: Obtain the optimized control sequence and response strategy matrix, construct an adaptive control law, obtain control instructions based on the optimized control sequence and response strategy matrix according to the adaptive control law, and generate the control instruction matrix of the system; Obtain the control instruction matrix and the pre-stored system state vector, construct the Lyapunov function, use the Lyapunov function to analyze the system stability based on the control instruction matrix and the system state vector, and obtain the stability evaluation index; Obtain the control instruction matrix, stability evaluation index and pre-stored real-time system status, build a control instruction correction function, input the control instruction matrix, stability evaluation index and real-time system status into the control instruction correction function to correct the control instruction in real time, and obtain the corrected control instruction sequence and execution feedback data.

Images