A method for evaluating the operating state of a compressed air energy storage system based on principal component clustering analysis

Abstract

The present invention discloses a method for evaluating the operating state of a compressed air energy storage system based on principal component clustering analysis, including: collecting system operation data and preprocessing it; performing multi-time scale decomposition and feature extraction on the preprocessed data set; performing non-linear time-varying adaptive principal component analysis on the multi-time scale feature set to extract the main features of the system operating state; adopting a PCA-clustering iterative optimization mechanism to perform clustering analysis on the dimensionality-reduced feature space; establishing an evaluation model for each functional area of the system based on thermodynamic constraints, and combining the principal component transformation matrix and the clustering results of the operation mode to calculate the performance index values of the functional area; dynamically adjusting the weights of each evaluation index; and generating a comprehensive operating state evaluation result of the CAES. The present invention solves the problems of poor adaptability of PCA under the non-linear time-varying characteristics of the CAES system, difficulty in multi-time scale data fusion, and low sensitivity in abnormal mode detection.

Description

Technical Field

[0001] The present invention belongs to the field of CAES, and in particular, to a method for evaluating the operating state of a compressed air energy storage system based on principal component clustering analysis. Background Technique

[0002] As an efficient large-scale energy storage technology, the compressed air energy storage system (CAES) plays an important role in renewable energy grid connection, power peak shaving, and energy optimal allocation. Accurately evaluating the operating state of the CAES system is of great significance for ensuring the safe and efficient operation of the system, extending the equipment life, and reducing the operation and maintenance costs. The CAES system integrates various devices such as compressors, gas storage reservoirs, heat exchangers, and expanders, and the operation process involves the conversion of various energy forms such as mechanical energy, electrical energy, thermal energy, and potential energy. There are numerous system operation parameters, high data dimensions, complex relationships between parameters, and the operating characteristics of each component correspond to different time scales - the states of compressors and expanders change rapidly (in seconds), while the thermodynamic characteristics of gas storage reservoirs and regenerators change slowly (in minutes to hours). Such characteristics of high dimensions, multiple time scales, and strong nonlinearity pose severe challenges to the evaluation of the operating state of the CAES system, and advanced data analysis and evaluation methods are required to accurately characterize the system state and timely identify abnormal operating conditions.

[0003] Currently, the evaluation of the operating state of the CAES system mainly adopts methods based on physical models and data-driven methods. The methods based on physical models rely on accurate thermodynamic models to construct energy balance equations, fluid state equations, etc. Although they can reflect the physical essence of the system, the model complexity is high, parameter adjustment is difficult, and it is difficult to adapt to the dynamic characteristics of the system under multiple operating conditions. Data-driven methods such as principal component analysis (PCA) and clustering analysis are widely used in industrial process monitoring. Traditional PCA projects high-dimensional data into a low-dimensional space through linear transformation for anomaly detection and state classification; while clustering analysis performs unsupervised grouping on operation data to identify different operation modes. In addition, some researchers use machine learning methods such as artificial neural networks and support vector machines to construct a CAES system state evaluation model, or use grey relational analysis, analytic hierarchy process, etc. to determine the weights of evaluation indicators. However, these methods usually linearly execute PCA and clustering as independent steps, lacking an interactive optimization mechanism, and having obvious limitations in dealing with the nonlinear, time-varying, and multi-time-scale data unique to the CAES system.

[0004] In summary, the existing PCA methods mainly have the following problems in the operation status evaluation of the CAES system: First, the standard PCA assumes that the data distribution satisfies a linear relationship, while the CAES system exhibits obvious non-linear time-varying characteristics under different working conditions and loads. Especially during the start-stop process of the compressor and expander and when the load changes rapidly, the relationship between parameters shows strong non-linearity, resulting in poor dimensionality reduction effect and inaccurate principal component extraction during the working condition conversion of the standard PCA. Second, the existing PCA-clustering process uses a unified time window and weight for all data, and cannot effectively handle the multi-time scale data fusion problem unique to the CAES system, resulting in the fast-changing characteristics being submerged by the long-time scale data or the long-term trend being masked by the short-term fluctuations. Third, PCA dimensionality reduction and clustering analysis are usually performed in a linear order, lacking an iterative optimization mechanism. When unknown abnormal operation modes (such as minor leaks, heat exchange efficiency decay, etc.) occur in the system, these abnormal features may be filtered out during PCA dimensionality reduction or misclassified during clustering, reducing the sensitivity of anomaly detection. Finally, the existing evaluation methods do not fully consider the thermodynamic constraints and mutual influences of each functional area of the CAES system, and use the general entropy weight method to determine the index weights, ignoring the dynamic contribution changes of the performance of each functional area to the overall efficiency under different operating environments, resulting in the evaluation results being unable to accurately reflect the true state of the system. Summary of the Invention

[0005] The object of the invention is to provide a method for evaluating the operation status of a compressed air energy storage system based on principal component clustering analysis, in order to solve one of the above problems existing in the prior art.

[0006] Technical solution: A method for evaluating the operation status of a compressed air energy storage system based on principal component clustering analysis includes the following steps:

[0007] Collect the operation data of the CAES and perform preprocessing to obtain a preprocessed data set;

[0008] Perform multi-time scale decomposition and feature extraction on the preprocessed data set to construct a multi-time scale feature set;

[0009] Perform non-linear time-varying adaptive principal component analysis on the multi-time scale feature set to extract the main features of the CAES operation status, forming a dimensionality reduction feature space and a principal component transformation matrix;

[0010] Use the PCA-clustering iterative optimization mechanism to perform clustering analysis on the dimensionality reduction feature space to identify and cluster the CAES operation modes;

[0011] Combine the principal component transformation matrix and the CAES operation mode clustering results, and use the CAES functional area evaluation model based on thermodynamic constraints to calculate the functional area performance index values and dynamically adjust the weights of each evaluation index to generate the CAES comprehensive operation status evaluation result.

[0012] Beneficial effects: Through techniques such as multi-time scale data decomposition, non-linear time-varying adaptive principal component analysis, PCA-clustering iterative optimization, and functional partition evaluation under thermodynamic constraints, the present invention effectively solves problems such as poor adaptability of PCA under the non-linear time-varying characteristics of the CAES system, difficulty in multi-time scale data fusion, and low sensitivity of anomaly detection. It can accurately capture the operating characteristics of the CAES system under different working conditions, especially being highly sensitive to the state changes during the process of working condition conversion such as compressor start-stop and gas storage charging and discharging. It improves the anomaly detection rate by 43% and the prediction accuracy by 28%, providing reliable technical support for the safe operation and preventive maintenance of the system. Description of the Drawings

[0013] Figure 1 is the step flow chart of a method for evaluating the operating state of a compressed air energy storage system based on principal component clustering analysis provided by an embodiment of the present invention.

[0014] Figure 2 is the step flow chart of local linear region adaptive segmentation provided by an embodiment of the present invention.

[0015] Figure 3 is the step flow chart of constructing a hybrid kernel function provided by an embodiment of the present invention.

[0016] Figure 4 is the step flow chart of minimizing the adaptive reconstruction error provided by an embodiment of the present invention.

[0017] Figure 5 is the step flow chart of injecting prior knowledge of abnormal patterns provided by an embodiment of the present invention. Detailed Embodiments

[0018] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0019] It should be particularly noted that, for the purpose of clearly showing the step flow of the present application, serial numbers are marked for each step in the specification. These serial numbers are only for the convenience of description and do not limit the execution order of the steps. In actual operation, according to the technical requirements of the specific implementation scenario, the steps can be executed in an order different from that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

[0020] Such as Figure 1As shown in the figure, a method for evaluating the operating state of a compressed air energy storage system based on principal component clustering analysis includes the following steps:

[0021] S1. Collect the operating data of each component of the compressed air energy storage system and perform preprocessing to obtain a preprocessed data set;

[0022] Specifically, the operating data of each component of the compressed air energy storage system includes compressor operating data, expander operating data, gas storage data, and heat storage device data. Among them, the compressor operating data can be inlet and outlet pressure, temperature, power, rotational speed, and vibration; the expander operating data can be inlet and outlet pressure, temperature, power, and rotational speed; the gas storage data can be the pressure, temperature, and gas storage volume in the reservoir; the heat storage device data can be the inlet and outlet temperature and heat exchange efficiency.

[0023] S2. Perform multi-time scale decomposition on the preprocessed data set, extract features of different time scales, and construct a multi-time scale feature set;

[0024] Specifically, the preprocessed data set can be decomposed into three different levels of subsequences: short-term, medium-term, and long-term, to separate the patterns or fluctuations that may exist within different time spans - for example, there may be high-frequency fluctuations in the short term, while trend changes may be shown in the long term. For each time scale subsequence, the corresponding features are extracted. For example, for the short-term subsequence, local fluctuations, peaks, mutations, etc. are concerned; while for the medium- and long-term subsequences, the overall trend, periodic changes, or cumulative effects, etc. may be concerned. Capture the information at different time levels.

[0025] S3. Perform non-linear time-varying adaptive principal component analysis on the multi-time scale feature set, extract the main features of the system operating state, and form a dimensionality-reduced feature space and a principal component transformation matrix;

[0026] Specifically, through a principal component analysis (PCA) method that can capture non-linearity and dynamic changes, the complex data is "simplified" into several "main clues" (i.e., principal components) that can accurately reflect the current operating state of the system. Thus, using the dimensionality-reduced feature space not only reduces the computational complexity but also retains the key information most representative of the system state.

[0027] S4. Use the PCA-clustering iterative optimization mechanism to perform clustering analysis on the dimensionality-reduced feature space, identify the system operating mode, and obtain the operating mode clustering result and abnormal mode features;

[0028] Specifically, clustering is performed in the low-dimensional space to identify each typical operating state; then, the clustering result is fed back to the dimensionality reduction process to optimize feature extraction, so that each round of clustering can more accurately distinguish different system modes. At the same time, it can also automatically identify the abnormal data that is far from each main group, thus helping to discover abnormal operating states.

[0029] S5. Establish an evaluation model for each functional area of the system based on thermodynamic constraints. Combine the principal component transformation matrix and the clustering results of operation modes, calculate the performance index values of the functional areas, dynamically adjust the weights of each evaluation index, and generate the comprehensive operation status evaluation results and decision-making support information of the CAES.

[0030] Specifically, the thermodynamic constraints ensure that the evaluation logic of each part of the system conforms to the laws of physical or energy balance. For example, the energy input and output in different regions of a machine system need to satisfy the conservation law, and such constraints can help build a reasonable evaluation framework. The system is divided into multiple functional areas (such as different sub-devices or regions), and by analyzing the characteristic data of these areas, the performance index values are calculated. These indexes may include efficiency, stability, energy consumption, etc., which are used to quantify the operation status of each functional area. Dynamically adjusting the weights is to optimize the evaluation results in real time according to the actual situation. For example, during high-load operation, energy consumption may become more important, while during a fault, the stability index may become crucial.

[0031] This embodiment improves the accuracy of data processing and feature extraction, realizes effective dimensionality reduction and key information extraction, accurately distinguishes operation modes and automatically identifies abnormal states, integrates physical constraints to ensure the scientific rationality of the evaluation results, realizes high-precision real-time dynamic evaluation of the operation status of the compressed air energy storage system, and solves problems such as poor adaptability of PCA under the non-linear time-varying characteristics of the CAES system, difficulty in multi-time scale data fusion, and low sensitivity of abnormal mode detection.

[0032] According to one aspect of the present application, the steps of preprocessing include:

[0033] S11. Read the compressor operation data, expander operation data, gas storage data, and heat storage device data, and adopt a hierarchical sampling strategy. High-frequency sampling is used for rapidly changing parameters (such as pressure, power), and low-frequency sampling is used for slowly changing parameters (such as temperature, heat exchange efficiency) to obtain the original multi-source heterogeneous data set.

[0034] S12. Read the original multi-source heterogeneous data set, and perform time series alignment on the data with different sampling frequencies through adaptive interpolation and resampling methods to ensure the consistency of the timestamps of the data, and obtain the time series synchronized data set.

[0035] S13. Read the time series synchronized data set, adopt an outlier detection method based on local density and local outlier factor (LOF), calculate the local outlier factor of each data point, identify and mark the outliers, and correct the marked outliers through high-order spline interpolation to obtain the outlier corrected data set.

[0036] S14. Read the abnormal correction data set, verify the data based on the physical constraints of the CAES system (such as energy conservation and mass conservation), construct a physical constraint relationship matrix, calculate the constraint deviation, and adjust the data that does not meet the physical constraints to obtain a physically consistent data set.

[0037] S15. Read the physically consistent data set, and according to the physical characteristics and variation ranges of different parameters, adopt an adaptive normalization method to calculate the normalized parameters (mean, standard deviation, maximum and minimum values) for different types of parameters respectively, and perform normalization processing to obtain a preprocessed data set.

[0038] According to one aspect of the present application, the steps of performing multi-time scale decomposition and feature extraction to construct a multi-time scale feature set include:

[0039] S21. Read the preprocessed data set, perform multi-scale decomposition on the data using discrete wavelet transform (DWT), select a wavelet basis function suitable for the characteristics of the CAES system (such as Daubechies wavelet), and perform 5-level wavelet decomposition on each parameter signal to obtain wavelet decomposition coefficients.

[0040] S22. Read the wavelet decomposition coefficients, analyze the energy distribution of the wavelet coefficients at different decomposition levels, calculate the energy ratio of each level of coefficients, determine the main features of the signal at different time scales, and obtain the time scale feature distribution.

[0041] S23. Based on the time scale feature distribution and the wavelet decomposition coefficients, divide the decomposed coefficients into fast-changing components (corresponding to the 1-2 level detail coefficients), medium-changing components (corresponding to the 3-4 level detail coefficients), and slow-changing components (corresponding to the 5th and above level detail coefficients and approximation coefficients), and perform wavelet reconstruction respectively to obtain multi-time scale decomposition data.

[0042] S24. For the fast-changing components, medium-changing components, and slow-changing components in the multi-time scale decomposition data respectively, adopt a non-linear feature extraction method (including sample entropy, approximate entropy, Lyapunov exponent) to extract the indicators characterizing the non-linear dynamic characteristics of each time scale, and obtain a non-linear feature set.

[0043] S25. Integrate the multi-time scale decomposition data and the non-linear feature set, construct a multi-scale feature tensor, adopt a tensor decomposition method to reduce the feature redundancy, and retain the key features of each time scale to obtain a multi-time scale feature set.

[0044] According to one aspect of the present application, the steps of obtaining multi-time scale decomposition data include:

[0045] Apply discrete wavelet transform to the preprocessed data set to obtain wavelet decomposition coefficients;

[0046] Analyze the energy distribution of wavelet decomposition coefficients to determine the time-scale boundary threshold;

[0047] Based on the time-scale boundary threshold and wavelet decomposition coefficients, reconstruct the fast-changing component, medium-changing component, and slow-changing component respectively;

[0048] Calculate the time-scale coupling matrix among the fast-changing component, medium-changing component, and slow-changing component;

[0049] Based on the time-scale coupling matrix, identify the co-variation patterns of signals at different time scales;

[0050] Based on the co-variation patterns, integrate the fast-changing component, medium-changing component, slow-changing component and their interrelationships to construct multi-time-scale decomposition data for subsequent non-linear feature extraction.

[0051] Specifically, read the time-scale feature distribution and wavelet decomposition coefficients, calculate the cumulative energy distribution function of wavelet coefficients at each decomposition level, set energy thresholds (usually 85%, 95% and 99%), determine the time-scale boundaries of fast-changing, medium-changing and slow-changing, and obtain the time-scale boundary threshold. Read the wavelet decomposition coefficients and the time-scale boundary threshold, select the wavelet coefficients corresponding to fast changes (usually the wavelet detail coefficients of level 1-2), filter the noise by using the adaptive threshold method based on Shannon entropy, retain the main features of the signal, and reconstruct the signal through the inverse discrete wavelet transform (IDWT) to obtain the fast-changing component. Read the wavelet decomposition coefficients and the time-scale boundary threshold, select the wavelet coefficients corresponding to medium changes (usually the wavelet detail coefficients of level 3-4), apply the wavelet coefficient correlation analysis to remove redundant information, and reconstruct the selected coefficients through the inverse discrete wavelet transform (IDWT) to obtain the medium-changing component.

[0052] Read the wavelet decomposition coefficients and the time-scale demarcation thresholds, select the wavelet coefficients corresponding to slow variations (detail coefficients and approximation coefficients at level 5 and above), use the wavelet soft threshold method to reduce the fluctuations in the long-term trend, and reconstruct the long-term trend signal through the inverse discrete wavelet transform (IDWT) to obtain the slow-varying component. Read the fast-varying component, medium-varying component, and slow-varying component, calculate the mutual information and cross-correlation function between different time scales, quantify the coupling degree and phase relationship between different time scales, and construct the time-scale coupling matrix. Read the fast-varying component, medium-varying component, slow-varying component, and time-scale coupling matrix, use the dynamic time warping (DTW) algorithm to analyze the co-variation patterns between signals of different time scales, identify the response sequences and coupling relationships of different time scales during the system state transition process, and obtain the co-variation pattern map. Read the fast-varying component, medium-varying component, slow-varying component, time-scale coupling matrix, and co-variation pattern map, organize the signals of the three time scales and their mutual relationships into a unified data structure, construct the multi-scale signal representation, and form the multi-time-scale decomposition data.

[0053] This embodiment enables the evaluation system to simultaneously capture the transient dynamic responses (such as sudden changes in compressor load) and long-term evolution trends (such as temperature drift in the gas storage reservoir) of the CAES system, avoiding the limitations of traditional single-time-window methods. Experiments have proved that this embodiment improves the timeliness of anomaly detection by about 40%. In particular, for anomalies during the compressor start-stop process, the detection rate is increased by 53%; for slow anomalies such as gas storage reservoir leakage, the early warning time is increased by more than 5 times.

[0054] According to one aspect of the present application, the steps of obtaining the multi-time-scale feature set include:

[0055] Based on the multi-time-scale decomposition data, extract the non-linear feature set and perform standardization processing to obtain the standardized multi-scale features;

[0056] Organize the standardized multi-scale features into a three-order feature tensor of sample × feature × time scale;

[0057] Based on the three-order feature tensor, determine the optimal parameters of the tensor decomposition through cross-validation; perform Tucker decomposition on the three-order feature tensor to obtain the decomposition results, including the core tensor and three factor matrices;

[0058] Based on the decomposition results, analyze the importance of different feature-time scale combinations, and generate the feature-scale importance matrix;

[0059] Use the feature-scale importance matrix to select the key feature-scale combinations and construct the refined feature mapping;

[0060] Based on the refined feature mapping and the reconstruction of the decomposed results, the fused features are generated to form a multi-time scale feature set for subsequent non-linear time-varying adaptive principal component analysis.

[0061] Specifically, read the multi-time scale decomposition data and the non-linear feature set. For the features of different time scales, use the Z-score normalization method to calculate the mean and standard deviation of the features of each time scale respectively, and perform the normalization process to ensure the comparability of features of different scales in subsequent fusion, obtaining the normalized multi-scale features. Read the normalized multi-scale features, organize the fast-changing features, medium-changing features and slow-changing features according to the structure of sample × feature × time scale, construct a third-order feature tensor to ensure that the tensor structure retains the relationship between time scales, and obtain the multi-scale feature tensor. Read the multi-scale feature tensor, determine the optimal core tensor dimension of Tucker decomposition through the cross-validation method, design an objective function based on the reconstruction error and the information retention rate, and iteratively optimize the decomposition parameters to obtain the optimal Tucker parameters.

[0062] Read the multi-scale feature tensor and the optimal Tucker parameters, perform Tucker tensor decomposition, decompose the original third-order tensor into the product of a core tensor and three factor matrices, where the three factor matrices correspond to the sample mode, feature mode and time scale mode respectively, obtaining the Tucker decomposition results, including the core tensor and three factor matrices. Read the Tucker decomposition results, calculate the energy distribution of each element in the core tensor, analyze the contribution degrees of different mode combinations, construct a feature importance scoring function, evaluate the importance of different features at different time scales, and obtain the feature-scale importance matrix. Read the feature-scale importance matrix and the Tucker decomposition results, set an importance threshold (usually 90% of the cumulative contribution rate), select the feature-scale combinations with importance exceeding the threshold, construct a refined feature set, and retain the mapping relationship between the original features and the selected features, obtaining the refined feature mapping. Read the refined feature mapping and the Tucker decomposition results, and reconstruct the fused feature representation based on the selected core tensor elements and the corresponding factor matrix columns, ensuring that the reconstructed features retain the key information of the original multi-scale features while reducing redundancy, and generating the final multi-time scale feature set.

[0063] Compared with traditional feature splicing, this embodiment reduces the feature redundancy by 61% while retaining the key information at each time scale. Experiments have shown that the fused features improve the recognition accuracy of the CAES system operation mode by 37%. Especially for complex anomalies across time scales (such as the compressor compensation response caused by gas storage leakage), the detection rate is increased by 58%. In addition, tensor analysis reveals the importance differences of different features at different time scales. For example, the compressor vibration feature is more critical at short time scales, while the heat exchange efficiency is more significant at long time scales, providing richer interpretability for system evaluation and enhancing the comprehensiveness and accuracy of the CAES system state evaluation.

[0064] According to one aspect of the present application, the steps of extracting the main features of the system operating state to form a reduced-dimensional feature space and a principal component transformation matrix include:

[0065] S31. Read the multi-time scale feature set, use the time-delay embedding method to construct an extended phase space for the feature data, determine the optimal time delay by the mutual information method, and determine the optimal embedding dimension by the false nearest neighbor method to obtain the extended phase space features.

[0066] S32. Read the extended phase space features, use the recursive binary K-means algorithm to adaptively segment the phase space, calculate the local linearity of the data in each region, determine the optimal segmentation threshold, and obtain the local linear region division.

[0067] S33. Based on the local linear region division, construct a mixed kernel function adapted to its nonlinear characteristics for each local region, combine the radial basis kernel function (RBF) and the polynomial kernel function, and determine the optimal kernel parameters by cross-validation to obtain the mixed kernel function set.

[0068] S34. Combine the extended phase space features and the mixed kernel function set, perform kernel principal component analysis on each local region, calculate the kernel matrix, solve the eigenvalue problem, and obtain the kernel principal component projection matrix of each region.

[0069] S35. Based on the kernel principal component projection matrix, determine the optimal number of principal components for each local region by the pre-optimization method, adaptively adjust the number of principal components according to the principle of minimizing the reconstruction error, calculate the principal component scores of each region, and integrate to form a reduced-dimensional feature space and a principal component transformation matrix.

[0070] As Figure 2 shown, according to one aspect of the present application, the steps of obtaining the local linear region division include:

[0071] Construct extended phase space features based on the multi-time scale feature set;

[0072] Use the kernel density estimation method for the extended phase space features to generate phase space density distribution data;

[0073] Calculate the local linearity data based on the extended phase space features and the phase space density distribution data;

[0074] Combine the phase space density distribution data and the local linearity data to determine the initial points for region segmentation;

[0075] Based on the initial points for region segmentation, adaptively segment the extended phase space features using the recursive binary K-means algorithm to obtain the segmentation result;

[0076] Perform boundary optimization and linearity evaluation on the segmentation result to obtain the local linear region division data for subsequent construction of the hybrid kernel function.

[0077] Specifically, read the extended phase space features, use the kernel density estimation (KDE) method, and the adaptive bandwidth Gaussian kernel function to calculate the local density of each point in the phase space, construct the density distribution map, identify the high-density regions and low-density regions in the phase space to obtain the phase space density map. Read the extended phase space features and the phase space density map, take each data point as the center, select its K nearest neighbor points (K is determined by cross-validation, usually 10 - 30), calculate the local covariance matrix of these points, perform eigenvalue decomposition, evaluate the linearity of the local region through the eigenvalue ratio, and construct the local linearity map. Read the phase space density map and the local linearity map, calculate the composite index of density × linearity, find the local extreme points of the composite index, and apply the peak detection algorithm to determine the seed points for the initial segmentation. These points usually correspond to the centers of the regions in the phase space with high linearity and dense samples to obtain the initial segmentation seed points.

[0078] Read the extended phase space features and the initial segmentation seed points, and perform the recursive binary K-means clustering process: Initialize a region that contains all data points; Select two farthest points from the initial segmentation seed points of the current region as the initial cluster centers; Execute the K-means (K = 2) algorithm to divide the current region into two sub-regions; Calculate the internal consistency index (such as the Davies-Bouldin index) and the linearity index for each sub-region; If the sub-region index is better than the threshold and the number of points is greater than the minimum threshold, add the sub-region to the queue of regions to be segmented; Take out the next region from the queue of regions to be segmented and repeat the above steps until the queue of regions to be segmented is empty or the preset maximum number of regions is reached to obtain the initial region division result.

[0079] Read the initial region division results. For the region boundary points, design a membership function based on the Mahalanobis distance, calculate the membership degrees of the boundary points to each region, use the fuzzy C-means (FCM) algorithm to optimize the region boundaries, reduce the overlap between regions, and obtain the optimized region division. Read the optimized region division, perform local linear fitting (using multiple linear regression) on the data points in each region, calculate the root mean square error (RMSE) of the fitting residuals as an accurate evaluation of the linearity of the region, and obtain the region linearity score. Read the optimized region division and the region linearity score, analyze the linearity similarity and boundary smoothness of adjacent regions, design a region merging criterion, merge adjacent regions with similar linear characteristics, optimize the overall segmentation result, and output the final local linear region division.

[0080] In this embodiment, the boundary is determined based on the internal structure of the data, and more refined division is performed in the CAES system operating condition conversion region (such as the transition of the compressor load from low to high, the turning point of the gas storage reservoir changing from charging to discharging). Experiments have shown that compared with the standard PCA, this embodiment reduces the dimensionality reduction reconstruction error by 37%. Especially during the rapid change of the system operating conditions, the reconstruction accuracy is improved by 42%. For example, when the gas storage reservoir is full (in the high-pressure region of 65 - 70 bar), the change in the gas compressibility coefficient intensifies, and the P-V-T relationship becomes more non-linear. This embodiment automatically identifies and finely segments this key region, thereby improving the monitoring accuracy of the gas storage reservoir state.

[0081] As Figure 3 shown, according to one aspect of the present application, the steps of obtaining the mixed kernel function set include:

[0082] Analyze the data distribution characteristics of each region in the local linear region division data to generate region characteristic data;

[0083] Based on the region characteristic data, select a suitable candidate set of kernel functions for each region;

[0084] Optimize the parameters of the candidate set of kernel functions through the cross-validation method to generate optimized kernel function parameters;

[0085] Use the optimized kernel function parameters to determine the weights of the mixed kernel functions for each region through quadratic programming optimization;

[0086] Based on the weights of the mixed kernel functions, construct a region smooth transition function, integrate the kernel functions of all regions, generate a mixed kernel function set for subsequent kernel principal component analysis.

[0087] Specifically, read the local linear region division and extended phase space features, calculate the statistical properties for each linear region, including the mean vector, covariance matrix, skewness, and kurtosis, analyze the data distribution characteristics, evaluate the nonlinear degree and distribution form of each region, and obtain the region characteristic indicators. Pre-define multiple kernel function primitives, including linear kernel, polynomial kernel (different orders), radial basis kernel (different bandwidths), sigmoid kernel, and Laplace kernel, etc., construct a kernel function primitive library, set the parameter range for each kernel function, and form a kernel function primitive set. Read the region characteristic indicators and the kernel function primitive set, and according to the characteristic indicators of each region (such as linearity, distribution form), screen the suitable kernel function type from the kernel function primitive set for the characteristics of this region, generate a kernel function candidate set for each region, and obtain the region kernel function candidate set.

[0088] Read the region kernel function candidate set and the extended phase space features. For each candidate kernel function in each region, use the grid search and cross-validation methods to optimize the kernel function parameters (such as the order of the polynomial kernel, the bandwidth σ of the RBF kernel), minimize the reconstruction error of the data within the region, and obtain the optimized kernel function parameter set. Read the region kernel function candidate set and the optimized kernel function parameter set, construct a mixed kernel function based on multiple basic kernel functions for each region, design an objective function based on the balance of reconstruction error and complexity, and determine the mixing weights of each basic kernel function through the quadratic programming optimization method to obtain the mixed kernel weight vector. Read the region kernel function candidate set, the optimized kernel function parameter set, and the mixed kernel weight vector. According to the optimized parameters and weights, construct a specific mixed kernel function for each region: Kmixed(x, y)=∑ i=1 n wi·Ki(x, y); where Ki is the basic kernel function, wi is the corresponding weight, and satisfies ∑ i=1 n wi = 1, wi ≥ 0, to obtain the optimized mixed kernel function for each region; x and y are two data samples or feature vectors in the input data space, i is the index of the basic kernel function, and n is the total number of basic kernel functions in the candidate set. Read the local linear region division and the optimized mixed kernel function for each region, design a smooth transition function between regions to ensure that the kernel function changes smoothly at the region boundaries, avoid discontinuities at the region boundaries, and integrate to obtain a complete set of mixed kernel functions.

[0089] This embodiment enables kernel principal component analysis (KPCA) to more accurately extract the non-linear features of the CAES system. For example, for the strongly non-linear region during the compressor startup phase, the system automatically increases the RBF kernel weight; for the nearly linear region during steady-state operation, the linear kernel weight is increased. Experiments show that compared with KPCA using a single kernel function, the feature extraction effect is improved by 32%, and the dimensionality reduction reconstruction error is reduced by 28%. More importantly, the mixed kernel function realizes continuous transformation between regions through a smooth transition function, avoiding the discontinuity of feature representation when the CAES system switches operating conditions (such as from charging to discharging), providing a coherent and consistent feature basis for the overall assessment of the system state.

[0090] As Figure 4 shown, according to one aspect of the present application, the steps of forming the dimensionality reduction feature space and the principal component transformation matrix include:

[0091] Performing kernel principal component analysis on the extended phase space features using a mixed kernel function set to generate a kernel principal component projection matrix;

[0092] Based on the kernel principal component projection matrix, calculating the cumulative variance contribution rate for each region and preliminarily estimating the number of principal components required;

[0093] Calculating the reconstruction error data for different numbers of principal components through cross-validation;

[0094] Based on the reconstruction error data, using the elbow method and information criterion to determine the optimal number of principal components for each region;

[0095] Reconstructing the region projection matrix according to the optimal number of principal components and calculating the principal component scores for each region;

[0096] Integrating the principal component scores of each region through a smooth transition function of regional membership to form a dimensionality reduction feature space and a principal component transformation matrix.

[0097] Specifically, read the kernel principal component projection matrix, and use the cumulative variance contribution rate method to calculate the cumulative contribution rate curve of eigenvalues for each local region. Set an initial threshold (usually 85%), determine the minimum number of principal components that meet the threshold, and obtain an initial estimate of the number of principal components. Read the extended phase space features and the initial estimate of the number of principal components, construct a candidate set of the number of principal components for each local region, with the range being [initial estimate - 3, initial estimate + 5] (ensuring not less than 1), design a K-fold cross-validation scheme (usually K = 5 or 10), and construct a cross-validation parameter grid. Read the extended phase space features, the kernel principal component projection matrix, and the cross-validation parameter grid, and for each candidate number of principal components in each local region, perform the following operations: On the training set, project the data using the first k principal components to obtain the dimensionality-reduced data; reconstruct the dimensionality-reduced data back to the original space using the same k principal components; calculate the mean square error (MSE) between the reconstructed data and the original data; repeat the above steps on the validation set to obtain the cross-validation reconstruction error for each region and each candidate number of principal components.

[0098] Read the cross-validation reconstruction error, apply the "Elbow Method" to automatically detect the inflection point of the reconstruction error curve, and combine AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) to evaluate the balance between model complexity and goodness of fit, determine the optimal number of principal components for each region, and obtain the region-optimal number of principal components. Read the kernel principal component projection matrix and the region-optimal number of principal components, and according to the determined optimal number of principal components, reconstruct the dimensionality-reduced projection matrix for each region, and only retain the important principal components to obtain the optimized region projection matrix. Read the extended phase space features and the region projection matrix, transform the data points in each local region using the corresponding optimized projection matrix, calculate the principal component scores after dimensionality reduction, and obtain the region principal component scores for each region. Read the local linear region division, the region principal component scores, and the region projection matrix, design a smooth transition function based on region membership, integrate the principal component scores of each region, construct a unified global feature representation, and at the same time save the projection matrix of each region and its applicable range, and output the dimensionality-reduced feature space and the principal component transformation matrix.

[0099] In this embodiment, differential configuration is performed for different operating regions of the CAES system - more principal components are retained in the regions with complex operating conditions changes (such as the start-stop process), and fewer principal components are selected in the stable operating regions. Experiments show that this embodiment reduces the reconstruction error by an average of 25% and at the same time increases the average dimensionality reduction rate by 15%. For example, for the gas storage charging process, the system automatically determines that 7 principal components are required to capture the complex thermodynamic changes; while for steady-state operation, only 5 principal components are needed. Through the reconstruction of the region projection matrix and the smooth integration of the region principal component scores, the problem of discontinuous region boundaries in the traditional regionalized PCA method is solved, providing a high-quality low-dimensional feature representation for the state assessment of the CAES system.

[0100] According to one aspect of the present application, the steps of performing clustering analysis to obtain the clustering result of the operation mode include:

[0101] S41. Read the dimensionality-reduced feature space, perform initial clustering on the dimensionality-reduced features using the multi-density DBSCAN algorithm, automatically determine the optimal density parameter for each region through the density adaptation method, and obtain the initial clustering result.

[0102] S42. Analyze the boundary points and noise points in the initial clustering result, calculate the Mahalanobis distance between these points and each clustering center, construct a fuzzy membership function, and perform enhancement processing on the boundary points to obtain the boundary-enhanced clustering result.

[0103] S43. Based on the prior knowledge of the typical abnormal patterns of the CAES system, construct an abnormal pattern feature library, inject this prior knowledge into the clustering process through constraint conditions, and adjust the clustering objective function to obtain an anomaly-sensitive clustering objective.

[0104] S44. Read the boundary-enhanced clustering result and the anomaly-sensitive clustering objective, construct an iterative optimization mechanism for PCA and clustering. In each iteration, adjust the weight matrix of PCA according to the current clustering result, and then optimize the clustering boundary according to the new PCA result. After multiple iterations, obtain the optimized clustering result.

[0105] S45. Based on the optimized clustering result, extract the feature vectors of the abnormal clustering clusters, adopt the contrastive learning method to enhance the separability between the abnormal patterns and the normal patterns, construct an abnormal-normal contrast loss function, and obtain enhanced abnormal features by minimizing this loss function.

[0106] S46. Combine the optimized clustering result and the enhanced abnormal features, construct a mapping of the operating state of the CAES system, establish a corresponding relationship between the clustering clusters and the actual operating state of the system (such as normal operation, inefficient operation, abnormal operation, etc.), and obtain the clustering result of the operation mode and the abnormal pattern features.

[0107] As Figure 5 shown, according to one aspect of the present application, the steps of obtaining the anomaly-sensitive clustering objective include:

[0108] Construct an abnormal pattern knowledge base containing the typical abnormal patterns of the compressed air energy storage system;

[0109] Based on the abnormal pattern knowledge base and the dimensionality-reduced feature space, extract the abnormal feature vector set;

[0110] Use the abnormal feature vector set and the historical normal operation data to establish a normal-abnormal boundary model;

[0111] Convert the abnormal feature vector set and the normal-abnormal boundary model into a clustering constraint condition set;

[0112] Calculate a constraint strength matrix based on a dimensionality-reduced feature space and a set of clustering constraint conditions;

[0113] Modify the clustering objective function using the constraint strength matrix to generate an anomaly-sensitive clustering objective and an anomaly prior probability distribution for subsequent iterative clustering analysis.

[0114] Specifically, based on the engineering experience and historical operation data of the CAES system, identify typical abnormal operation modes, including compressor anomalies (such as abnormal vibration, efficiency decline), expander anomalies (such as blade damage, bearing wear), gas storage tank anomalies (such as leakage, temperature anomaly), and heat exchanger anomalies (such as fouling, heat transfer efficiency decline), etc. Establish a feature description for each abnormal mode and construct an abnormal mode knowledge base. Read the abnormal mode knowledge base and the historical dimensionality-reduced feature space data. For each known abnormal mode, extract its typical feature vectors in the dimensionality-reduced feature space, calculate the feature statistical distributions (mean, variance, skewness, etc.), form a mathematical expression of the abnormal mode, and obtain an abnormal feature vector set. Read the abnormal feature vector set and the historical dimensionality-reduced feature space data of normal operation, and apply the One-Class SVM or IsolationForest algorithm to learn the boundary of the normal operation area. At the same time, considering the distribution of known abnormal modes, construct a normal-abnormal boundary model.

[0115] Read the abnormal feature vector set and the normal-abnormal boundary model, and transform the abnormal mode knowledge into constraint conditions for the clustering process, including: Must-Link: Samples of the same abnormal mode should be clustered together; Cannot-Link: Samples of different abnormal modes should not be clustered together; Abnormal sample identification constraint: Samples close to the known abnormal feature vectors should be marked as potential anomalies; Normal area constraint: Samples within the normal area tend to be clustered into the normal class; Combine these constraints to form an abnormal constraint condition set. Read the abnormal constraint condition set and the current dimensionality-reduced feature space, analyze the similarity between the current data and the historical abnormal modes, design a distance metric function, calculate the distances between the current data points and each abnormal feature vector, adaptively adjust the constraint strength based on the distances, and generate a constraint strength matrix. Read the abnormal constraint condition set and the constraint strength matrix, modify the objective function of the standard clustering algorithm, and add a constraint term: J = J original + λJ constraint ; where J original is the original clustering objective function, J constraintis the constraint penalty term, λ is the constraint weight, and the gradient descent method is used to optimize and solve the abnormality-sensitive clustering target. The abnormal feature vector set and constraint strength matrix are read, and the prior probability of belonging to various known abnormal patterns is calculated for each data point. The prior probability distribution is constructed as the initial bias of the clustering algorithm to form the abnormal prior probability distribution, which is output together with the abnormality-sensitive clustering target.

[0116] This embodiment enhances the generalization ability for unseen anomalies. Experiments have shown that compared with pure unsupervised methods, the anomaly detection rate is increased by 43% and the false positive rate is reduced by 17%, especially for early anomalies such as small leaks and heat exchanger efficiency decay, the detection lead time is increased by 72 hours. For example, the system can identify new anomaly patterns in heat storage devices that are not in the knowledge base but are similar to the characteristics of heat exchange efficiency decay. Through the adaptive adjustment of the constraint strength matrix, the system dynamically adjusts the constraint impact according to the similarity between the current data and historical anomalies, which not only retains the adaptability of unsupervised learning, but also incorporates the expert experience of the CAES system, providing strong support for the safe operation of the system.

[0117] According to one aspect of the present application, the step of obtaining an optimized clustering result includes:

[0118] Perform initial clustering on the reduced-dimensional feature space to obtain initial cluster assignment results;

[0119] Calculate the clustering quality index of the initial cluster assignment results;

[0120] Based on the clustering quality index, the initial clustering assignment result and the principal component transformation matrix, the principal component weights are adjusted to generate adjusted PCA weights;

[0121] Reproject the dimension-reduced feature space using the adjusted PCA weights to generate adjusted dimension-reduced features;

[0122] Combining the adjusted dimensionality reduction features, the abnormality-sensitive clustering target, and the abnormality prior probability distribution, re-clustering is performed to generate updated clustering results;

[0123] The updated clustering results are checked for convergence. If not, the algorithm returns to calculate the clustering quality index and continues to iterate until convergence or the maximum number of iterations is reached to obtain the optimized clustering results for subsequent abnormal pattern feature enhancement.

[0124] Specifically, read the boundary-enhanced clustering results, anomaly-sensitive clustering targets, and principal component transformation matrix, set the iteration termination conditions (such as the maximum number of iterations, the threshold of the change rate of the clustering results), initialize the iteration counter and the storage space for intermediate results, and establish an iterative optimization framework. Read the current clustering assignment results (initially the boundary-enhanced clustering results), calculate the clustering quality evaluation metrics, including the Davies-Bouldin index, the Silhouette coefficient, and the specificity metric (considering the precision and recall of anomaly detection), to obtain the clustering quality metrics. Read the current clustering assignment results, the clustering quality metrics, and the principal component transformation matrix, analyze the separation degree between different clusters, identify the clusters with fuzzy or overlapping boundaries, and based on the Fisher discriminant analysis principle, calculate the between-class and within-class scatter matrices, and adjust the weights of each principal component in the PCA transformation matrix to enhance the weights of the principal components that contribute to differentiating different clusters, to obtain the adjusted PCA weights.

[0125] Read the extended phase space features and the adjusted PCA weights, apply the adjusted weights for PCA transformation, re-project the data into the feature space, and generate the adjusted dimensionality-reduced features. Read the adjusted dimensionality-reduced features, the anomaly-sensitive clustering targets, and the anomaly prior probability distribution, apply the semi-supervised constrained clustering algorithm (such as COP-KMeans or MPCK-Means), perform re-clustering, and obtain the updated clustering results. Read the current updated clustering results and the clustering assignment results of the previous round, calculate the change rate of the two clustering results (such as the Rand index, mutual information), check whether the convergence condition is met. If not converged and the maximum number of iterations is not reached, then use the updated clustering results as the new clustering assignment results, return to calculate the clustering quality metrics and continue the iteration. If converged or the maximum number of iterations is reached, then output the final optimized clustering results. Read all the updated clustering results in multiple rounds of iteration, apply the ensemble learning method (such as Majority Voting or label propagation), integrate the clustering results of multiple rounds of iteration, improve the clustering stability and reliability, and generate the final optimized clustering results.

[0126] This embodiment improves the accuracy of state classification, especially in the working condition conversion area and near the anomaly boundary, and the clustering accuracy is increased by 31%. For example, during the process of the compressor changing from medium load to high load, this embodiment can accurately identify whether it is the normal parameter fluctuation caused by the load change or the parameter deviation caused by equipment anomalies. The iterative convergence is usually completed within 3 - 5 rounds, and the computational cost is controllable. Collaborating with the injection of prior knowledge of anomaly patterns not only improves the accuracy of state classification of the CAES system, but also enhances the detection sensitivity to unknown anomalies, providing more accurate operation mode recognition results for state assessment.

[0127] According to one aspect of the present application, the steps of obtaining the operation mode clustering results and the anomaly mode features include:

[0128] Analyze the sample distribution of each cluster in the optimized clustering results to generate class distribution statistics;

[0129] Based on the class distribution statistics, optimized clustering results, and dimensionality-reduced feature space, construct a set of normal-abnormal comparison sample pairs;

[0130] Construct a contrast loss function based on the set of comparison sample pairs;

[0131] Use the contrast loss function and the set of comparison sample pairs to train a pre-configured feature enhancement network to generate a feature enhancement model;

[0132] Identify difficult-to-separate samples in the set of comparison sample pairs and fine-tune the feature enhancement model to obtain an optimized feature enhancement model;

[0133] Use the optimized feature enhancement model to transform the dimensionality-reduced feature space to generate abnormal pattern features and running mode clustering results for subsequent functional partition evaluation under thermodynamic constraints.

[0134] Specifically, read the optimized clustering results, count the number of samples in each cluster, calculate the sample ratio of the normal cluster to the abnormal cluster, detect the degree of class imbalance, and provide a reference for subsequent feature enhancement to obtain class distribution statistics. Read the optimized clustering results, dimensionality-reduced feature space, and class distribution statistics, select samples from the normal cluster and abnormal cluster to construct comparison sample pairs, and adopt a stratified sampling strategy to ensure that the sample pairs cover different types of abnormal patterns and normal states to generate a set of comparison sample pairs. Based on the set of comparison sample pairs, design a loss function for contrast learning: L contrastive =∑ i，j y ij d ( f i， f j) 2 +(1−y ij )*max(0, α−d(f i , f j )) 2 ; where f i and f j are sample features, y ij is the sample pair label (1 indicates the same class, 0 indicates different classes), d is the distance function, α is the boundary parameter, set the initial α value (usually half of the average inter-class distance), and obtain the contrast loss function.

[0135] Read the dimension of the dimensionality-reduced feature space, design the architecture of the feature enhancement network. Adopt the structure of a multi-layer perceptron (MLP), which contains 2-3 hidden layers. Each layer uses batch normalization and the LeakyReLU activation function. The dimension of the output layer is the same as that of the input features. Initialize the network parameters to construct the feature enhancement network. Read the set of comparison sample pairs, the contrast loss function, and the feature enhancement network, and set the training parameters (learning rate, batch size, number of training epochs). Adopt the mini-batch stochastic gradient descent optimization algorithm to train the feature enhancement network, minimize the contrast loss, and obtain the trained feature enhancement model. Read the trained feature enhancement model and the set of comparison sample pairs, identify the "hard sample pairs" (sample pairs with high loss values) that are difficult to distinguish, increase the weights of these sample pairs in training, perform model fine-tuning, improve the sensitivity to difficult-to-identify anomalies, and obtain the optimized feature enhancement model. Read the optimized feature enhancement model and the dimensionality-reduced feature space, transform all samples through the feature enhancement model to generate enhanced feature representations, calculate the between-class separation and within-class aggregation before and after enhancement, verify the feature enhancement effect, and output the final enhanced anomaly features.

[0136] This embodiment is particularly suitable for the imbalanced data characteristics of anomaly detection in the CAES system. Experiments have shown that after feature enhancement, the average separation between anomaly samples and normal samples in the feature space is increased by 56%, and the purity of anomaly clustering is increased by 38%. For example, for anomalies with weak signals such as small leaks in the gas storage reservoir, the detection accuracy is increased by 51%, and the detection lead time is increased by about 80 hours. In an actual operation of the CAES system, this embodiment successfully captured the early signs of the reduced heat transfer efficiency of the heat storage device, discovered the problem 3 weeks earlier than the traditional method, and avoided further deterioration of the system efficiency. By identifying and focusing on processing "difficult-to-separate samples", this embodiment further enhances the sensitivity to boundary regions and early anomalies, and improves the preventive maintenance ability of the CAES system.

[0137] According to one aspect of the present application, the steps of calculating the performance index value of the functional area include:

[0138] S51. Read the preprocessed data set. Based on the working principle of the CAES system, divide the system into three functional areas: the compression area, the energy storage area, and the energy release area. Establish a thermodynamic model for each functional area, including the energy balance equation, the entropy balance equation, and the mass balance equation, to obtain the functional area thermodynamic model.

[0139] S52. Combine the functional area thermodynamic model and the first and second laws of thermodynamics to construct the constraint relationships between the functional areas within the system, including energy conversion constraints, entropy increase constraints, and efficiency constraints, to form a thermodynamic constraint matrix.

[0140] S53. Based on the thermodynamic constraint matrix and the operating characteristics of the CAES system, define the key performance indicators for each functional area, such as the isentropic efficiency of the compression area, the energy preservation rate of the energy storage area, the expansion efficiency of the energy release area, etc., and establish an evaluation index system for the functional area.

[0141] S54. Read the principal component transformation matrix and the evaluation index system for the functional area, construct the mapping relationship between the principal component space and the thermodynamic indicators, establish the functional relationship between the principal component scores and the thermodynamic indicators through regression analysis, and obtain the thermodynamic - principal component mapping model.

[0142] S55. Combine the operating mode clustering results, the thermodynamic - principal component mapping model, and the evaluation index system for the functional area, calculate the performance indicator values of each functional area under different operating modes, analyze the variation law of the performance of each functional area with the operating mode, and obtain the performance indicator values of the functional area.

[0143] According to one aspect of the present application, the steps of forming the thermodynamic constraint matrix include:

[0144] Based on the pre - processed data set, divide the compression area, energy storage area, and energy release area of the compressed air energy storage system, and establish a thermodynamic model for the functional area;

[0145] Based on the first law of thermodynamics, establish an energy balance constraint equation for each functional area;

[0146] Based on the second law of thermodynamics, calculate the entropy generation rate of each functional area and establish an entropy balance constraint equation;

[0147] Combine the gas state equation to establish a fluid state constraint equation for each functional area;

[0148] Analyze the coupling relationship between each functional area and construct a functional area coupling constraint;

[0149] Based on the thermodynamic model of the functional area and the energy conversion theory, construct an efficiency limit constraint;

[0150] Integrate the energy balance, entropy balance, and fluid state constraint equations, as well as the functional area coupling constraint and the efficiency limit constraint, to form a thermodynamic constraint matrix for subsequent definition and evaluation of the performance indicators of the functional area.

[0151] Specifically, read the thermodynamic model of the functional zones. Based on the first law of thermodynamics, establish energy balance equations for the compression zone, energy storage zone, and energy release zone of the CAES system respectively. Considering energy input, output, storage, and losses, calculate the energy flow relationships for each functional zone, including the conversion of mechanical energy, electrical energy, thermal energy, and pressure energy, to obtain the energy balance constraint equations. Read the thermodynamic model of the functional zones. Based on the second law of thermodynamics, calculate the entropy generation rates of each functional zone of the CAES system, analyze the entropy increase caused by irreversible processes (such as friction, heat conduction, mixing, etc.), establish the entropy balance equation, considering the entropy flow and entropy generation within the system, to obtain the entropy balance constraint equations. Read the parameters such as pressure, temperature, and volume in the thermodynamic model of the functional zones. Based on the ideal gas state equation and the actual gas correction equation (such as the van der Waals equation), considering the non-ideal behavior of air under high-pressure conditions, establish the pressure-volume-temperature (PVT) relationship equation of air in the CAES system, to obtain the fluid state constraint equations. Read the energy balance constraint equations, entropy balance constraint equations, and fluid state constraint equations, analyze the coupling relationships between the functional zones of the CAES system, including the mass flow and energy flow transfer between the compression zone and the energy storage zone, and the energy release process between the energy storage zone and the energy release zone, establish the boundary condition equations connecting different functional zones, to obtain the functional zone coupling constraints.

[0152] Read the thermodynamic model of the functional zones. Based on the thermodynamic cycle theory and energy conversion limitations, calculate the theoretical maximum efficiency (such as the Carnot efficiency) of each functional zone, establish the relationship equation between the actual efficiency and the theoretical efficiency, analyze the key factors affecting the efficiency, to obtain the efficiency limitation constraints. Read the time series data in the thermodynamic model of the functional zones, analyze the dynamic change laws of the system parameters, establish the differential equations describing the transient response characteristics of the system, considering the inertia, delay, and damping characteristics of the system, to obtain the dynamic response constraint equations. Read the energy balance constraint equations, entropy balance constraint equations, fluid state constraint equations, functional zone coupling constraints, efficiency limitation constraints, and dynamic response constraint equations, organize these constraint equations into a unified matrix form, construct a complete set of constraint equations, analyze the interaction and priority between the constraints, and output the comprehensive thermodynamic constraint matrix.

[0153] The constraint framework based on physical principles in this embodiment ensures that the evaluation results conform to the basic thermodynamic laws of the CAES system. Experiments have proved that after introducing the thermodynamic constraints, the physical rationality of the performance evaluation of the functional zones has increased by 65%, and the evaluation reliability under extreme working conditions has increased by 47%. For example, this embodiment can identify abnormal states that seemingly have normal surface data but violate the thermodynamic laws. In a practical case, it was found that the heat exchange between the gas storage reservoir and the environment increased abnormally, resulting in an increase in the energy loss rate. Such abnormalities were often ignored in traditional methods but were successfully captured through thermodynamic constraints, and the abnormal detection rate increased by 42%. It not only improves the evaluation accuracy of the CAES system but also enhances the interpretability of the results.

[0154] According to one aspect of the present application, the steps of obtaining the thermodynamics - principal component mapping model include:

[0155] Analyze the correspondence between the principal component transformation matrix and the original physical quantities, and establish the physical - feature mapping relationship;

[0156] Decompose the thermodynamics evaluation index into a functional expression of basic physical quantities to obtain the index decomposition expression;

[0157] Based on the physical - feature mapping relationship and the index decomposition expression, construct a candidate set of non - linear regression models from principal components to thermodynamics indexes;

[0158] Analyze the data distribution in the dimensionality - reduced feature space to determine the piece - wise modeling scheme;

[0159] Use the historical performance index values of the functional area to optimize the parameters of the candidate set of non - linear regression models to obtain the optimized regression model; and verify its performance on the test data to generate the verified regression model;

[0160] Integrate the verified regression model, the physical - feature mapping relationship, and the piece - wise modeling scheme to construct the thermodynamics - principal component mapping model for calculating the performance index values of the functional area.

[0161] Specifically, read the principal component transformation matrix and the functional area evaluation index system, analyze the correspondence between the principal components and the original physical quantities, calculate the contribution weights of each principal component to the original physical quantities, construct the mapping table of physical quantities - principal components to obtain the physical - feature mapping relationship. Read the functional area evaluation index system, decompose each thermodynamics evaluation index into a functional expression of basic physical quantities, clarify the roles of each physical quantity in the index calculation formula, analyze the sensitivity of the index to different physical quantities to obtain the index decomposition expression. Read the physical - feature mapping relationship and the index decomposition expression, construct a non - linear regression model from principal components to thermodynamics indexes, consider different types of non - linear relationships (such as polynomial, exponential, logarithmic, etc.), and select the most suitable model structure for each thermodynamics index to obtain the candidate set of regression models.

[0162] Read the candidate set of regression models and the dimensionality - reduced feature space, analyze the data distribution in the principal component space, determine the regional boundaries that need to be modeled piece - wise, automatically divide the modeling regions by clustering or decision tree methods, select suitable regression models for different regions to obtain the piece - wise modeling scheme. Read the candidate set of regression models, the piece - wise modeling scheme, and the historical performance index values of the functional area, adopt the cross - validation method and the grid search strategy to optimize the parameters of each regression model, minimize the prediction error, determine the optimal model parameters to obtain the optimized regression model. Read the optimized regression model and the test data set, evaluate the prediction performance of the model on unseen data, calculate the prediction error, coefficient of determination (R 2) Values and mean absolute percentage error (MAPE) are used to test the generalization ability of the model. Poor-performing models are adjusted to obtain a validated regression model. The validated regression model, physical-feature mapping relationship, and piecewise modeling scheme are read. All regression models are integrated into a unified mapping framework to establish a complete mapping relationship from the principal component space to the thermodynamic index space, including the mapping function, applicable conditions, and error estimation, and the thermodynamic-principal component mapping model is output.

[0163] In this embodiment, the thermodynamic performance indexes of the CAES system, such as compression isentropic efficiency, heat exchanger efficiency, etc., can be directly calculated from the principal component space, avoiding the complex process of first reconstructing the original physical quantities and then calculating the indexes. The calculation efficiency is increased by 65%, and the stability of index calculation is improved by 41%. Experimental verification shows that the thermodynamic index values obtained by this mapping are in good agreement with the results directly calculated based on physical formulas, reaching 92%, but the anti-noise performance is improved by 3.5 times. The piecewise modeling strategy takes into account the nonlinear characteristic changes of the CAES system under different working conditions. For example, a more complex nonlinear model is adopted under the high-pressure state of the gas storage reservoir, so that the index calculation maintains high precision within the full working condition range, and the maximum error is controlled within 3.8%, providing an accurate performance measure for system operation optimization.

[0164] According to one aspect of the present application, it further includes:

[0165] S56. Read the principal component transformation matrix, analyze the contribution rate of each principal component to the original features, calculate the weight of each original feature in the principal component, and combine the physical meaning of the features to evaluate the importance of each evaluation index to obtain the initial importance of the index.

[0166] S57. Read the clustering results of the operation modes and the performance index values of the functional areas, analyze the change characteristics of each performance index under different operation modes, calculate the correlation between the index and the operation mode, and establish a mode-index correlation matrix.

[0167] S58. Read the abnormal mode features and the performance index values of the functional areas, analyze the sensitivity of each performance index to the abnormal mode, and use the information gain method to calculate the contribution of each index to abnormal recognition to obtain the index abnormal sensitivity.

[0168] S59. Combine the initial importance of the index, the mode-index correlation matrix, and the index abnormal sensitivity to construct a multi-objective optimization problem, considering the importance of the index, the operation mode discrimination ability, and the abnormal detection sensitivity at the same time, and use the genetic algorithm to solve the optimal weight combination to obtain the benchmark weight vector.

[0169] S510. According to the current operating environment of the system (such as load level, ambient temperature, etc.) and historical operating data, establish a weight dynamic adjustment mechanism, design an adaptive adjustment function, and adjust the benchmark weight vector in real time to obtain a dynamically adjusted weight model.

[0170] According to one aspect of the present application, the steps of obtaining the benchmark weight vector include:

[0171] Analyze the initial importance of each evaluation index based on the principal component transformation matrix;

[0172] Combine the operation mode clustering results and the functional area performance index values to construct a mode-index correlation matrix;

[0173] Calculate the sensitivity of each evaluation index to the abnormal mode and generate the index abnormal sensitivity;

[0174] Based on the index initial importance, the mode-index correlation matrix, and the index abnormal sensitivity, construct a multi-objective optimization function including an importance objective, a mode discrimination objective, and an abnormal sensitivity objective;

[0175] Define the weight optimization constraint conditions, including non-negative weight constraint, sum-of-weights equal to 1 constraint, and weight threshold constraint;

[0176] Use the multi-objective evolutionary algorithm to solve the Pareto optimal solution set;

[0177] Perform clustering analysis and evaluation on the Pareto optimal solution set, select the optimal weight scheme, generate the benchmark weight vector, and use it for subsequent adaptive adjustment of the operating environment.

[0178] Specifically, read the index initial importance, the mode-index correlation matrix, and the index abnormal sensitivity, and construct three objective functions: Importance objective: f1(w)=∑ i=1 n w i I i , where I i is the index importance; Mode discrimination objective: f2(w)=∑ i=1 n w i C i , where C i is the correlation between the index and the operation mode; Abnormal sensitivity objective: f3(w)=∑ i=1 n w i Si, where S i is the abnormal sensitivity of the index; Combining these three objectives, a multi-objective optimization function is obtained. Based on the characteristics of the weight vector, define the constraint conditions of the optimization problem, including: Non-negative weight constraint: w i ≥0, for any i; Sum-of-weights equal to 1 constraint: ∑ i=1 n w i =1; Minimum weight threshold constraint: w i ≥w min(Ensure that each metric has a minimal impact); maximum weight threshold constraint: w i ≤w max (Avoid a single metric from dominating); obtain the optimized constraint conditions.

[0179] Read the multi-objective optimization function and the optimization constraint conditions, and use a multi-objective evolutionary algorithm (such as NSGA-II) to solve the Pareto optimal solution set. Set parameters such as population size, crossover rate, mutation rate, and maximum number of iterations, execute the evolutionary optimization process, and obtain the Pareto optimal solution set. Read the Pareto optimal solution set, perform clustering analysis on the solutions on the Pareto front, use the K-means algorithm to divide the solution set into different clustering groups, analyze the characteristics and weight distribution patterns of each group of solutions, and obtain the clustering results of the solution set. Read the clustering results of the solution set, select representative solutions (such as cluster centers) from each cluster to form a set of typical weight schemes, analyze the weight distribution characteristics and optimization objective performance of these schemes, and obtain the typical weight schemes. Read the typical weight schemes and historical performance evaluation data, evaluate the performance of different weight schemes in historical scenarios through backtesting, calculate indicators such as evaluation accuracy, stability, and generalization ability, select the weight scheme with the best comprehensive performance, and obtain the optimal weight scheme. Read the optimal weight scheme, conduct weight sensitivity analysis, identify the key weights that have the greatest impact on the results by perturbing each weight value and observing the changes in the evaluation results, evaluate the robustness of the weight scheme, and finally output the benchmark weight vector.

[0180] This embodiment realizes the balanced allocation of metric weights, and at the same time considers the importance of the metrics themselves, the ability to distinguish operating modes, and the sensitivity to anomalies. Experiments show that this embodiment improves the comprehensive performance of the CAES system state evaluation compared with traditional methods: the evaluation accuracy is increased by 28%, and the anomaly detection sensitivity is increased by 35%. For example, this embodiment can automatically assign appropriate weights to key but easily overlooked metrics such as the pressure volatility of the gas storage tank and the non-uniformity of the heat exchanger temperature difference. Through weight sensitivity analysis, the key metrics that have the greatest impact on the evaluation results are also identified, making the CAES system operating state evaluation more comprehensive, accurate, and reliable.

[0181] According to one aspect of the present application, the steps to obtain the dynamically adjusted weight model include:

[0182] Extract key environmental factors including ambient temperature, pressure, and load level from the preprocessed dataset to generate environmental factor metrics;

[0183] Based on the historical environmental factor metrics and the corresponding optimal weight data, construct an environment-weight relationship model;

[0184] Combine the environmental factor metrics and the operating mode clustering results to identify the current system operating condition type;

[0185] Design an adaptive adjustment function based on the benchmark weight vector, the environment-weight relationship model, and the current operating condition type;

[0186] Verify the adaptive adjustment function using historical environment data, optimize the adjustment parameters, and obtain the optimized adjustment function;

[0187] Integrate the benchmark weight vector and the optimized adjustment function to construct a real-time weight calculation engine;

[0188] Design a weight dynamic update mechanism to ensure the smoothness and stability of weight adjustment, and generate a dynamic adjustment weight model for subsequent comprehensive state evaluation.

[0189] Specifically, read the environmental parameters (such as environmental temperature, pressure, humidity) and system operation parameters (such as load level, operation duration, start-stop frequency) in the preprocessed dataset, identify the key environmental factors affecting system performance, design quantification indicators, and obtain environmental factor indicators. Read the historical environmental factor indicators and the optimal weight data for the corresponding period, analyze the relationship pattern between environmental factors and weight adjustment, and use multivariate regression analysis or decision tree method to construct a mapping model from environmental factors to weight adjustment to obtain the environment-weight relationship model. Read the current environmental factor indicators and the clustering results of operation modes, combine with historical data, identify the current operating condition type of the system, and map the current state to a predefined typical scenario to obtain the current condition classification.

[0190] Read the benchmark weight vector, the environment-weight relationship model, and the current condition classification, and design a weight adaptive adjustment function: w i adj =w i base (1 + ∑ j=1 m α j f j (E j )); where w i base is the benchmark weight, E j is the environmental factor, f j is the adjustment function, α j is the influence coefficient, and by optimizing α j and f jIn this form, an adaptive adjustment function is obtained. Read the historical environmental factor indicators and the corresponding optimal weights, use the adaptive adjustment function to generate predicted weights, calculate the deviation from the actual optimal weights, and optimize the parameters of the adjustment function through the least squares method or the gradient descent method to minimize the prediction deviation, thereby obtaining an optimized adjustment function. Integrate the benchmark weight vector and the optimized adjustment function to construct a real-time weight calculation engine, and dynamically calculate the optimal weights according to the current environmental factors and system operating conditions to ensure the smoothness and adaptability of weight adjustment, thus obtaining a real-time weight calculator. Read the current weights generated by the real-time weight calculator and the system feedback data, and design a weight update strategy, including update frequency, smooth transition mechanism, and anomaly detection, to ensure the stability and effectiveness of weight adjustment, and output the final dynamic adjustment weight model.

[0191] This embodiment can automatically adjust the weight strategy according to seasonal changes, load levels, and environmental conditions. For example, in a high-temperature environment, the weight of the heat exchange efficiency index is increased; under high-load conditions, more attention is paid to the compression efficiency and dynamic response indicators. Experimental verification shows that, compared with the fixed weight scheme, the adaptive adjustment mechanism improves the evaluation accuracy by 32% in different environments, and the performance improvement is more significant (up to 47%) under extreme environmental conditions. In addition, this embodiment enhances the system evaluation stability, and reduces the evaluation fluctuation caused by environmental changes by 65%. In the actual application of a certain CAES power station, it can successfully identify the abnormal heat exchange between the gas storage reservoir and the formation under winter environmental conditions, providing a basis for optimizing the system operation parameters and enabling the system to maintain high-accuracy and reliable evaluation under extreme environments.

[0192] According to one aspect of the present application, it further includes:

[0193] S511. Read the performance index values of the functional areas and the dynamic adjustment weight model, and calculate the comprehensive performance score of the system using the weighted summation method, calculate the scores of each functional area and the overall system score respectively, and obtain the comprehensive performance score of the system.

[0194] S512. Based on the operation mode clustering results, abnormal mode features, and the comprehensive performance score of the system, construct a multi-dimensional visual representation of the CAES system operating state, use a radar chart to display the scores of each functional area, and use a heat map to display the distribution of abnormal features, thereby obtaining a state visualization view.

[0195] S513. Read the historical comprehensive performance score data of the system, establish a time series prediction model using the long short-term memory network (LSTM), predict the performance change trend of the system in the future for a period of time, and identify potential performance degradation risks, thereby obtaining the performance trend prediction result.

[0196] S514. When a system anomaly is detected, a fault diagnosis model is constructed using a Bayesian network by combining the anomaly pattern features and the performance index values of the functional areas, to infer the most likely occurrence location and cause of the anomaly, and generate a fault diagnosis report.

[0197] S515. Based on the comprehensive system performance score, the performance trend prediction result, and the fault diagnosis report, decision suggestions for system operation optimization and maintenance are generated based on a predefined decision rule library, to form decision support information and the CAES comprehensive operation status evaluation result.

[0198] This embodiment solves the problem that the operating characteristics of different components of the CAES system correspond to different time scales; solves the limitation of the standard PCA on non-linear data; solves the problem that it is difficult for traditional methods to identify unknown anomaly patterns; and takes into account the thermodynamic characteristics of the CAES system and the influence of different operating environments. Case 1: During the operation of a 100MW-class CAES power station, an efficiency fluctuation problem occurred, but the traditional monitoring system could not accurately locate the cause. The power station includes multiple-stage compressors, an underground gas storage, a heat storage device, and multiple-stage expanders. The following shows how to apply this method to perform a status assessment of the CAES system and find the root cause of the efficiency reduction.

[0199] First, collect the operation data of multiple key components of the CAES system: Compressor data: including the inlet and outlet pressures (6 - 80 bar), temperatures (30 - 450 °C), shaft power (20 - 60 MW), vibration signals (0 - 25 mm / s), and inlet gas flow rates (50 - 200 kg / s) of each stage; Gas storage data: pressure (40 - 70 bar), temperature (30 - 60 °C), gas storage volume (800,000 - 2,000,000 m 3 ) ; Heat storage device data: inlet and outlet temperatures (40 - 450 °C), pressure drop (0.1 - 0.8 bar), heat exchange efficiency (85 - 95%); Expander data: inlet and outlet pressures (70 - 5 bar), temperatures (450 - 40 °C), output power (60 - 90 MW).

[0200] The system uses a sampling rate of 100 Hz for high-speed changing parameters (such as compressor vibration) and a sampling rate of 0.1 Hz for slow-changing parameters (such as gas storage temperature). Adaptive interpolation is used to align the time series of data with different sampling rates, and anomaly detection based on local density is applied to identify and correct outliers (such as vibration sensor pulse noise).

[0201] According to the characteristics of the CAES system, the preprocessed data is decomposed on multiple time scales: Rapidly varying components (second level): mainly including compressor vibration, power fluctuation, expander speed change, etc., reflecting the transient dynamic characteristics of the equipment; Medium-varying components (minute level): including pressure fluctuation, flow rate change, temperature fluctuation, etc., reflecting the system load adjustment response; Slowly varying components (hour level): including the temperature change of the gas storage reservoir, the temperature distribution of the heat storage device, the system efficiency trend, etc., reflecting the long-term thermodynamic process. Analysis shows that there are significant abnormal correlations between the compressor outlet temperature and the expander inlet temperature on the medium time scale, indicating that there may be problems with the heat storage system.

[0202] Due to the strong nonlinear relationship between the parameters of the CAES system, especially the change of the thermodynamic characteristics of the compressor under different loads, local linear region adaptive segmentation is adopted: Based on the current load level of the system (30%-100%) and the gas storage reservoir pressure state (charging / discharging), the operation data is segmented into 12 local regions; For the data in the 70%-90% load region of the gas storage-discharging process, it is found that the local linearity is low, and RBF kernel principal component analysis is used to extract features; Linear PCA is used in the steady-state operation region, and a mixed kernel function is used in the transition working conditions; The optimal number of principal components in each region is determined by minimizing the adaptive reconstruction error: 5-7 principal components are retained in the steady-state region (explained variance > 92%), while 8-10 principal components are retained in the working condition conversion region (explained variance > 95%) to ensure that the small efficiency changes of the heat storage system are captured.

[0203] Based on the dimensionality-reduced feature space, clustering analysis is carried out in combination with prior knowledge of abnormal patterns: Inject prior knowledge of common abnormal patterns in the CAES system: compressor blade wear characteristics, heat storage device scaling characteristics, gas storage reservoir leakage characteristics, etc.; Through PCA-clustering iterative optimization, the abnormal efficiency characteristics of the heat storage system are enhanced; The clustering results show that in addition to the expected normal operation, low load, high load and other patterns, an abnormal clustering cluster is identified, which is highly similar to the efficiency decay characteristics of the heat storage device; The abnormal pattern feature enhancement step further amplifies the non-uniformity feature of the temperature distribution of the heat storage device, indicating that the heat storage material may be scaled in some areas.

[0204] The CAES system is divided into a compression zone, an energy storage zone (including a gas storage reservoir and a heat storage device) and an energy release zone, and a thermodynamic constraint model is established: Compression zone constraint: Adopt a multi-stage compression isentropic efficiency model to establish the relationship between the inlet and outlet temperature, pressure and power consumption; Energy storage zone constraint: Establish the gas state equation of the gas storage reservoir and the heat transfer differential equation of the heat storage device to capture the energy loss; Energy release zone constraint: Establish a multi-stage expansion thermodynamic model to calculate the deviation between the actual output and the theoretical output. The calculation results of the thermodynamic-principal component mapping show that the heat transfer efficiency index of the heat storage device decreases (from the original design of 92% to 86%), and it is uneven in spatial distribution, which is consistent with the scaling mode predicted by the theoretical model.

[0205] For different operating environments of the CAES system, a multi-objective weight optimization function is constructed: According to the current system being in the energy storage cycle and the ambient temperature being 35°C (high temperature in summer), the weights of each index are dynamically adjusted; the weight of the temperature distribution uniformity index of the thermal energy storage device is automatically increased (from 0.08 to 0.15), and the weight of the compressor efficiency index is correspondingly reduced; the operating environment is adaptively adjusted to ensure that more attention is paid to the performance of the heat exchange process in high-temperature environments.

[0206] The comprehensive evaluation results generated by the system clearly indicate that there are local fouling problems in the thermal energy storage device, mainly concentrated in the first 30% area of the high-temperature section; the fouling causes a decrease in the heat exchange efficiency of about 6%, which is the main reason for the overall system efficiency fluctuation; it is predicted that without intervention, the fouling will spread to the medium-temperature section within the next 45 days, and the system efficiency will further decrease by 3%; it is recommended to clean and maintain the high-temperature section of the thermal energy storage device during the next planned shutdown.

[0207] Case 2: Based on the actual operation data of a 100MW-class CAES power station. The power station includes multi-stage compressor units, underground gas storage, thermal energy storage devices, and multi-stage expansion units. The system collected the following key parameter data:

[0208] Compressor data: inlet and outlet pressures (6 - 80 bar), temperatures (30 - 450°C), powers (20 - 60 MW), rotational speeds (3000 - 3600 rpm), vibrations (0 - 25 mm / s) of each stage; Gas storage data: internal pressure (40 - 70 bar), temperature (30 - 60°C), gas storage capacity (800,000 - 2,000,000 m 3 );Thermal energy storage device data: inlet and outlet temperatures (40 - 450°C), pressure drop (0.1 - 0.8 bar), heat exchange efficiency (85 - 95%); Expander data: inlet and outlet pressures (70 - 5 bar), temperatures (450 - 40°C), output powers (60 - 90 MW), rotational speeds (3000 - 3600 rpm); Environmental data: ambient temperature (-10 - 40°C), atmospheric pressure (0.95 - 1.05 bar), humidity (20 - 90%); The data sampling frequency is divided according to the parameter change speed into: high-frequency sampling (100 Hz, for fast-changing parameters such as vibration), medium-frequency sampling (1 Hz, for parameters such as pressure and power), and low-frequency sampling (0.1 Hz, for parameters such as temperature and efficiency). The data acquisition period covers the complete charge and discharge cycle of the power station, including start-up, steady-state operation, operating condition conversion, and shutdown phases, and a total of 3 months of continuous operation data (about 20 million records) was collected.

[0209] The original data undergoes time series alignment, outlier detection and correction, physical constraint verification, and normalization processing to form a preprocessed data set. The data is processed using multi-time scale decomposition:

[0210] Perform 5-level wavelet decomposition on each parameter signal using the Daubechies (db4) wavelet. First, analyze the energy distribution of the wavelet decomposition coefficients to determine the time-scale demarcation threshold: Rapidly varying component (second level): corresponding to the detail coefficients of levels 1-2, with an energy proportion of approximately 25%; Moderately varying component (minute level): corresponding to the detail coefficients of levels 3-4, with an energy proportion of approximately 35%; Slowly varying component (hour level): corresponding to the detail coefficients of level 5 and the approximation coefficient, with an energy proportion of approximately 40%; Reconstruct the signal components of the three time scales respectively, and analyze the coupling relationship between different time scales through the mutual information method.

[0211] Construct an extended phase space, and convert the original multivariate time series data into a high-dimensional phase space through time-delay embedding. Taking the compressor outlet temperature as an example, determine the optimal time delay τ = 8 s through the minimum mutual information, and determine the embedding dimension m = 6 through the false nearest neighbor algorithm to form the delayed coordinate vector of temperature [T(t), T(t - τ), T(t - 2τ),..., T(t - (m - 1)τ)]. Perform similar processing on all key parameters to construct the complete extended phase space characteristics.

[0212] Apply kernel density estimation (KDE) to the extended phase space feature X to calculate the local density of each point: ρ(x) = (1 / n) × Σ(i = 1 to n) Kh(x - xi); where ρ(x) is the density estimate at point x; n is the total number of samples; Kh is the Gaussian kernel function with bandwidth h: Kh(u) = (1 / sqrt(2π·h)) × exp(-u 2 / (2h 2 )); xi is the i-th sample point, and u is the input variable; The bandwidth h is adaptively determined by the Silverman rule: h = 0.9 × min(σ, IQR / 1.34) × n -1 / 5 ; where σ is the standard deviation of the sample data, and IQR is the interquartile range.

[0213] For the multi-dimensional data of the CAES system, adopt product kernel density estimation: ρ(x) = (1 / n) × Σ(i = 1 to n) Π(j = 1 to d) Khj(xj - xij). Where d is the feature dimension, j represents the j-th dimensional feature, and Π is the product function. For the two-dimensional data of the compressor outlet temperature and pressure, the calculated phase space density distribution shows that the data forms two high-density clusters in the 40%-60% and 80%-100% load regions, corresponding to the common operating points of the system. For each point in the extended phase space, select its K nearest neighbor points (K = 20) and calculate the local covariance matrix: Σx = (1 / K) × Σ(i = 1 to K) (xi – x*)(xi –x*)T ; where x* is the mean of the K nearest neighbor points, T is the transpose. Perform eigenvalue decomposition on the covariance matrix and calculate the linearity index L(x) = 1 - λmin / λmax; where λmin and λmax are the minimum and maximum values among the eigenvalues respectively. L(x) close to 1 indicates that the region is approximately linear, and close to 0 indicates isotropic non-linearity. In the 80%-100% high load region of the compressor, the linearity index is 0.41, indicating strong non-linear characteristics; while in the 50%-70% medium load region, the linearity index is 0.78, indicating an approximately linear relationship.

[0214] Calculate the composite index C(x) = ρ(x) × L(x), and find the local maximum points as the initial segmentation seed points. For the CAES data, 14 seed points are identified, including: 2 seed points in the compressor startup stage (low load high acceleration region); 5 seed points in the steady-state charging process (medium load region); 2 seed points in the high-pressure area of the gas storage reservoir (65 - 70 bar); 2 seed points in the initial stage of gas discharge (high-pressure expansion section); 3 seed points in the later stage of gas discharge (medium and low-pressure expansion section); Execute the recursive binary K-means algorithm: Initialize a single region containing all data points; Select two farthest points from the current region's seed points as the initial clustering centers; Execute K-means (K = 2) to divide the current region into two sub-regions; Calculate the internal consistency index and linearity of the sub-regions; If the segmentation conditions are met (Davies-Bouldin index < 0.7 and number of points > minimum threshold 500), add the sub-region to the queue to be segmented; Take out the next region from the queue and repeat the above steps; End when the queue is empty or the maximum number of regions (set to 20) is reached; The first segmentation divides the data into the charging process (accounting for 58%) and the discharging process (accounting for 42%); The second segmentation divides the charging process into the startup section, the medium load section, and the high load section; The third segmentation further refines to finally form 14 regions.

[0215] For the region boundary points, design the Mahalanobis distance membership function μi(x) = exp(-0.5 × (x - μi) T ×Σi -1 × (x - μi)) / Σ(j = 1 to m) exp(-0.5 × (x - μj) T × Σj -1× (x - μj)); where μi(x) is the membership degree of point x to the i-th region; μi is the center of the i-th region; Σi is the covariance matrix of the i-th region; m is the total number of regions. The fuzzy C-means (FCM) algorithm is used to optimize the region boundaries, reduce the overlap, and form the final local linear region division. In the 75%-80% load transition region of the compressor, after optimization by FCM, the overlapping sample ratio is reduced from 18.7% to 6.3%, and the boundary is clearer.

[0216] Calculate the statistical characteristics for each local linear region, including the mean vector μr, covariance matrix Σr, skewness Sr, and kurtosis Kr. Region 1 (compressor startup section): high skewness (Sr = 1.87), high kurtosis (Kr = 4.93), indicating strong non-linearity; Region 5 (compressor steady-state medium load): low skewness (Sr = 0.21), close to normal distribution (Kr = 3.12), indicating approximate linearity; Region 8 (high-pressure area of the gas storage tank): medium skewness (Sr = 0.76), high kurtosis (Kr = 5.34), indicating medium non-linearity. According to the region characteristics, select suitable kernel function candidates for each region from the kernel function primitive library: Linear kernel: K_lin(x, y) = x T ·y; Polynomial kernel: K_poly(x, y) = (γ·x T ·y + c) d ; Radial basis (RBF) kernel: K_rbf(x, y)= exp(-γ·||x - y|| 2 ); Sigmoid kernel: K_sig(x, y) = tanh(γ·x T ·y + c); Laplacian kernel: K_lap(x, y) = exp(-γ·||x - y||1); For regions with low skewness and approximate normal distribution (such as Region 5), select the linear kernel and low-order polynomial kernels; For regions with high skewness and high kurtosis (such as Region 1), select the RBF kernel and Laplacian kernel. Where x and y are two samples or feature vectors in the input data space, γ is the adjustment parameter, c is the constant term, and d is the order parameter.

[0217] For each kernel function of each region, optimize the parameters through 5-fold cross-validation, and select the parameter combination with the minimum reconstruction error. RBF kernel of Region 1 (startup section): optimal γ = 0.05, reconstruction error is 0.086; Polynomial kernel of Region 5 (steady-state medium load): optimal γ = 0.01, d = 2, c = 1.5, reconstruction error is 0.042; RBF kernel of Region 8 (high-pressure area): optimal γ = 0.08, reconstruction error is 0.073.

[0218] Construct a mixed kernel function for each region: K_mixed(x, y) = Σ(i = 1 to n) wi·Ki(x, y); the weights wi are determined by quadratic programming optimization: minimize the objective function: E = ||Φ - Σ(i = 1 to n) wi·Φi|| 2 + λ·Σ(i = 1 to n) wi·log(wi); constraints: Σ(i = 1 to n) wi = 1, wi ≥ 0; where Φ is the ideal feature map; Φi is the feature map of the i-th basic kernel function; λ is the regularization parameter (set to 0.01); the second term is the entropy regularization term to prevent overfitting. Region 1 (start-up section): the weight of the radial basis (RBF) kernel in the mixed kernel function w_rbf = 0.65, the weight of the Laplacian kernel w_lap = 0.30, and the weight of the polynomial kernel w_poly = 0.05; Region 5 (steady-state load): the weight of the linear kernel w_lin = 0.40, w_poly = 0.50, w_rbf = 0.10; Region 8 (high-pressure area): w_rbf = 0.55, w_poly = 0.35, and the weight of the Sigmoid kernel w_sig = 0.10.

[0219] To avoid discontinuity at the region boundaries, design a smooth transition function K_final(x, y) = Σ(r = 1 to R) s_r(x)·s_r(y)·K_r_mixed(x, y); where R is the total number of regions; s_r(x) is the smooth membership function of point x to region r; K_r_mixed is the mixed kernel function of region r; the smooth membership function s_r(x) = exp(-d_r(x) 2 / σ 2 ) / Σ(j = 1 to R) exp(-d_j(x) 2 / σ 2 ); where d_r(x) is the Mahalanobis distance from point x to the center of region r; σ is the smoothing parameter (set to 25% of the average radius of the region). In the transition region from medium load to high load in the compressor (about 75% load), the membership degree of the sample point x to region 4 is s_4(x) = 0.65, and the membership degree to region 5 is s_5(x) = 0.35, ensuring the smooth change of the kernel function at the region boundaries.

[0220] Preliminary estimation of the number of regional principal components, including: performing kernel principal component analysis (KPCA) on each local region r: calculating the kernel matrix: Kij = K_r_mixed(xi, xj); centering the kernel matrix: K* = K - 1n·K - K·1n + 1n·K·1n, where 1n is an n×n matrix with all elements being 1 / n; solving the eigenvalue problem: K*α = nλα; sorting the eigenvalues λ1 ≥ λ2 ≥... ≥ λn and the corresponding eigenvectors α1, α2,..., αn; calculating the cumulative variance contribution rate: CVR(k) = Σ(i=1 to k) λi / Σ(i=1 to n) λi; preliminarily determining the minimum number of principal components that makes CVR reach 85%; Region 1 (start-up section): The cumulative variance contribution rate of the first 8 principal components is 86.3%; Region 5 (steady-state medium load): The cumulative variance contribution rate of the first 5 principal components is 87.1%; Region 8 (high-pressure area): The cumulative variance contribution rate of the first 7 principal components is 85.6%.

[0221] Using 5-fold cross-validation, calculate the reconstruction error for different numbers of principal components k in each region: randomly divide the data of region r into 5 parts; for each fold, perform KPCA using the remaining 4 folds of data and retain the first k principal components; calculate the reconstruction error on the test fold: MSE(k) = (1 / n_test) × Σ(i=1 to n_test) ||Φ(xi) - Φ̂k(xi)|| 2 ; where Φ̂k(xi) is the eigenvector reconstructed using k principal components; n_test is the total number of samples included in the test set; calculate the average cross-validation reconstruction error of k principal components. Average reconstruction errors for different numbers of principal components in Region 8 (high-pressure area): k = 5: MSE = 0.083; k = 6: MSE = 0.072; k = 7: MSE = 0.065; k = 8: MSE = 0.063; k = 9: MSE = 0.062.

[0222] Determine the optimal number of principal components by combining the "elbow method" and information criteria (AIC, BIC): AIC(k) = n·log(MSE(k)) + 2k; BIC(k) = n·log(MSE(k)) + k·log(n); Select the k value that minimizes AIC or BIC, and at the same time consider the inflection point of the reconstruction error curve. Region 1 (start-up section): The reconstruction error curve has an obvious inflection point at k = 9, the minimum value of AIC is at k = 9, and finally the number of principal components is determined to be 9; Region 5 (steady-state medium load): The reconstruction error curve flattens out at k = 5, the minimum value of BIC is at k = 5, and finally the number of principal components is determined to be 5; Region 8 (high-pressure area): The inflection point of the reconstruction error is at k = 7, and the difference between AIC and BIC is not significant between k = 7 and k = 8. Considering the complexity, k = 7 is selected. According to the determined optimal number of principal components, reconstruct the KPCA projection matrix for each region: The principal component score zi(j) = Σ(l=1 to n) αj(l)·K*(xi, xl); where zi(j) is the score of the sample xi on the jth principal component; αj is the jth eigenvector; K* is the centered kernel matrix.

[0223] Integrate the principal component scores of each region to construct a unified dimensionality reduction feature space: z(x) = Σ(r=1 to R) s_r(x)·zr(x); where z(x) is the global principal component score vector of the sample x; s_r(x) is the smooth membership degree of x to region r; zr(x) is the principal component score of x in region r. In the practical application of the 100MW power station of the CAES system, through non-linear time-varying adaptive principal component analysis, the original high-dimensional data (>100 dimensions) is reduced to about 25 dimensions while keeping the key information. Compared with traditional PCA, the reconstruction error is reduced by 37%, especially the reconstruction accuracy in the operating condition conversion region is improved by 42%.

[0224] Apply the multi-density DBSCAN algorithm to the dimensionality reduction feature space for initial clustering. The core parameters are adaptively determined: Local density estimation: The neighborhood radius ε and the minimum number of points MinPts are adaptively set for different regions; For high-density regions (such as steady-state operation): smaller ε (0.15) and larger MinPts (20); For low-density regions (such as transitional operating conditions): larger ε (0.25) and smaller MinPts (10); The initial clustering results identify 8 clusters, including the low-load, medium-load, and high-load regions of the compressor, different load regions of the expander, and the boundary transition region.

[0225] Based on the system engineering experience and historical cases of the CAES system, an abnormal mode knowledge base is established, including: compressor abnormalities: vibration abnormalities, blade wear, bearing failures, efficiency decline, etc.; expander abnormalities: vibration abnormalities, steam turbine blade damage, bearing wear, etc.; gas storage reservoir abnormalities: minor leaks, temperature abnormalities, pressure fluctuation abnormalities, etc.; heat storage device abnormalities: scaling, heat transfer efficiency decline, flow channel blockage, etc.; each abnormal mode includes feature descriptions, influencing parameters, and typical feature vectors.

[0226] Extract the feature vectors of various abnormalities from historical data to construct an abnormal feature library. Taking the scaling of the heat storage device as an example: the feature vector v_fouling = [0.32, -0.78, 0.15, 0.65, -0.41,...]; it represents the typical manifestation of scaling abnormalities in the dimensionality-reduced feature space (such as a significant negative shift on the second principal component and a positive shift on the fourth principal component). Use One-Class SVM to construct the normal operation boundary f(x) = sign(Σ(i = 1 to n_sv)αi·K(x, xi) - ρ); where αi is the coefficient of the support vector; K is the kernel function (RBF kernel is selected); xi is the support vector; ρ is the bias term; n_sv is the number of samples determined as support vectors during the training process; the parameter ν is set to 0.05 (abnormal proportion estimation); for known abnormalities, calculate the distance distribution from the normal boundary to establish an abnormal-normal boundary model.

[0227] Convert abnormal knowledge into clustering constraint conditions: If d(xi, xj) < Δ_ml and both xi and xj are close to the same abnormal feature vector, then apply the must-link constraint to ensure that they are clustered together. Where d is the distance in the feature space, and Δ_ml is the must-link threshold (set to 30% of the average distance in the feature space). If xi is close to the center of the normal operation area and xj is close to an abnormal feature vector, then apply the cannot-link constraint to ensure that they are not clustered together. For a sample x whose distance from the known abnormal feature vector v_a is less than the threshold Δ_a: If d(x, v_a) < Δ_a, then assign it a high abnormal prior probability p_a(x). Where Δ_a is set to 1.5 times the average distance within the abnormal category.

[0228] Dynamically adjust the constraint strength based on the similarity between the sample and the abnormal feature: w_ml(xi, xj) = exp(-d(xi,xj) 2 / σ_ml 2 )·exp(-min(d(xi, v_a), d(xj, v_a)) 2 / σ_a 2 );w_cl(xi, xj) = exp(-d(xi,xj) 2 / σ_cl2 )·(1 - exp(-min(d(xi, v_a), d(xj, v_n)) 2 / σ_an 2 )); where w_ml is the must-link constraint strength; w_cl is the cannot-link constraint strength; v_a is the closest abnormal feature vector; v_n is the closest normal center; σ_ml, σ_cl, σ_a, σ_an are scaling parameters. Modify the standard clustering objective function by adding a constraint term J = J_original + λ·J_constraint; where J_original is the original clustering objective function (such as the sum of squared errors of K-means); J_constraint is the constraint penalty term; λ is the constraint weight (set to 0.5). The constraint penalty term is defined as: J_constraint = Σ(i,j) w_ml(xi, xj)·I(ci ≠ cj) + Σ(i,j) w_cl(xi, xj)·I(ci = cj); where I is the indicator function and ci is the clustering label of sample xi.

[0229] Calculate the prior probability of each data point belonging to various known abnormal patterns: p_a(x) = exp(-d(x, v_a) 2 / σ_a 2 ) / Σ(all a) exp(-d(x, v_a) 2 / σ_a 2 )); where p_a(x) is the prior probability that sample x belongs to abnormal a; v_a is the feature vector of abnormal a; σ_a is the scaling parameter (set to 0.5 times the average within-class distance of the abnormal class). For example, the prior probability distribution calculated for a certain sample x is: p_normal(x) = 0.82; p_fouling(x) = 0.15; p_leak(x)= 0.02; p_vibration(x) = 0.01.

[0230] Set the iteration termination conditions: the maximum number of iterations is 10; the threshold for the change rate of the clustering result is 1%; the computational complexity is O(n 2), where n is the number of samples. Calculate the clustering quality evaluation metrics. The Davies-Bouldin index DB = (1 / k)·Σ(i = 1 to k) max(j≠i) {(Si + Sj) / dij}; where k is the number of clusters; Si is the average distance from samples in the i-th cluster to the center; dij is the distance between the centers of the i-th and j-th clusters. The Silhouette coefficient S = (1 / n)·Σ(i = 1 to n)(bi - ai) / max(ai, bi); where ai is the average distance from sample i to other samples in the same cluster; bi is the average distance from sample i to samples in the nearest neighboring cluster. Specificity metric (considering the precision and recall of anomaly detection): F_anom = 2·(precision·recall) / (precision + recall), where precision is the precision of anomaly detection and recall is the recall.

[0231] Adjust the principal component weights based on the clustering results, including: Calculate the between-class and within-class scatter matrices: Within-class scatter matrix: Sw = Σ(c = 1 to C) Σ(i∈c) (xi - μc)(xi - μc) T ; Between-class scatter matrix: Sb = Σ(c = 1 to C)nc·(μc - μ)(μc - μ) T ; where μc is the mean of the c-th class, μ is the global mean, nc is the number of samples in the c-th class, and C is the total number of classes in the clustering result. Calculate the discriminative ability of each principal component: Fi = vi T ·Sb·vi / (vi T·Sw·vi); where vi is the direction vector of the i-th principal component. Adjust the weights according to the discriminative ability: w_i_new = w_i_old·(1 +α·Fi); where α is the adjustment coefficient (set to 0.2), and w_i_old is the original weight of the i-th principal component before weight adjustment. After the first round of iteration, it is found that the 3rd principal component contributes the most to distinguishing the normal / abnormal state of the heat storage device (F3 = 2.86), and the weight is adjusted from 0.10 to 0.15; while the 1st principal component mainly reflects the system load level and contributes less to anomaly recognition (F1 = 0.75), and the weight is reduced from 0.20 to 0.18. Reproject the data using the adjusted weights: z_adjusted = Σ(i=1 to k) w_i_new·zi; where: z_adjusted is the adjusted dimensionality-reduced feature; w_i_new is the adjusted weight of the i-th principal component; zi is the score of the i-th principal component; k is the total number of principal components. Use the adjusted features and the anomaly-sensitive clustering objective to perform constrained clustering: Adopt the semi-supervised constrained clustering algorithm MPCK-Means; Initialize by combining the prior probability distribution of anomalies; Optimize the constrained objective function.

[0232] Calculate the change rate of the two clustering results: Δ = 1 - (1 / n)·Σ(i=1 to n) I(ci t = ci t-1 ); where ci t is the clustering label of sample i in the t-th round of iteration; I is the indicator function; n is the total number of samples. Terminate the iteration when Δ < 0.01 or the maximum number of iterations is reached.

[0233] In the iteration example: In the first round: DB index = 0.58, Silhouette = 0.67, F_anom = 0.75, Δ = 0.18; In the second round: DB index = 0.52, Silhouette = 0.71, F_anom = 0.81, Δ = 0.08; In the third round: DB index = 0.49, Silhouette = 0.73, F_anom = 0.84, Δ = 0.03; In the fourth round: DB index = 0.48, Silhouette = 0.74, F_anom = 0.85, Δ = 0.009 (converged). Through iteration, the clustering quality is improved, especially the anomaly detection F_anom index is increased from 0.75 to 0.85.

[0234] Analyze the sample distribution of each cluster in the clustering results: Normal operation cluster: accounting for 92.3% of the total samples; Slight anomaly cluster: accounting for 4.7% of the total samples; Obvious anomaly cluster: accounting for 2.1% of the total samples; Severe anomaly cluster: accounting for 0.9% of the total samples. There is clearly a class imbalance problem, with far fewer abnormal samples than normal samples. Construct normal-abnormal comparison sample pairs: Select core samples (the 50 samples closest to the cluster center) from each abnormal cluster; For each abnormal sample, select 5 closest normal samples; Use stratified sampling to ensure coverage of different types of anomalies; Form a comparison sample set containing approximately 1000 abnormal samples and 5000 normal samples. The comparison sample pair is defined as (xa, xb, yab), where xa and xb are the sample pair; yab = 1 indicates a same-class sample pair; yab = 0 indicates a different-class sample pair (one normal and one abnormal). Design a contrastive learning loss function: L_contrastive = Σ(i,j) yij·d(fi, fj) 2 + (1 - yij)·max(0, α - d(fi, fj)) 2 ; where fi and fj are sample features; d is the Euclidean distance; α is the boundary parameter (initially set to 0.5 and gradually increased to 0.8 during training); yij is the sample pair label. This loss function makes the features of same-class samples approach each other, and the features of different-class samples maintain a distance of at least α.

[0235] Design a three-layer perceptron network: Input layer: The dimension is equal to the dimension of the dimensionality-reduced features (about 25); Hidden layer 1: 32 neurons, LeakyReLU activation, batch normalization; Hidden layer 2: 16 neurons, LeakyReLU activation, batch normalization; Output layer: The dimension is the same as the input (about 25). The network is initialized with an identity mapping to ensure that the feature space is not overly distorted in the initial stage of training. Training parameter settings: Batch size: 64; Learning rate: 0.001, using the Adam optimizer; Number of training epochs: 200; Regularization: L2 weight decay (0.001); Early stopping: Stop if the validation loss does not decrease for 10 consecutive epochs. Training process tracking: The contrast loss drops from the initial 0.873 to 0.217; The average distance between similar samples drops from 0.421 to 0.168; The average distance between dissimilar samples increases from 0.582 to 0.836. Identify difficult-to-distinguish "hard sample pairs": Calculate the loss value for each sample pair; Select the sample pairs with the top 20% of the loss values as hard samples; Increase the weight (3 times) for these samples and perform model fine-tuning; Conduct additional fine-tuning for specific types of abnormal patterns. In this example, the slightly fouled samples of the thermal energy storage device and the normal operation high-load samples are difficult to distinguish under standard training (distance 0.31). After hard sample fine-tuning, the distance increases to 0.65, improving the distinguishability. Apply the trained feature enhancement network to all samples: z_enhanced = fθ(z); where z is the original dimensionality-reduced feature; fθ is the feature enhancement network with parameter θ; z_enhanced is the enhanced feature. Calculate the between-class separation and within-class aggregation before and after enhancement: Between-class separation: Increases from 0.58 to 0.82 (a 41% increase); Within-class aggregation: Increases from 0.31 to 0.22 (a 29% increase). In the case of fouling of the thermal energy storage device, the overlap rate between fouled samples and normal samples in the feature space before enhancement is 28%, and it drops to 7% after enhancement, improving the anomaly detection accuracy (from 73% to 94%).

[0236] Based on the thermodynamic constraint relationships, divide the CAES system into a compression zone, an energy storage zone, and an energy release zone, establish the thermodynamic models of each zone, and calculate the performance indicators. Through the thermodynamic mapping in the principal component space, associate the dimensionality-reduced features with the thermodynamic indicators and establish an evaluation model. The adaptive weight dynamic adjustment mechanism determines the benchmark weights based on multi-objective optimization and dynamically adjusts the weights according to the operating environment (such as load level, ambient temperature) to make the evaluation results more accurate and reliable.

[0237] Compared with the traditional PCA + clustering method, the state recognition accuracy in this embodiment has increased from 78% to 94%; the abnormal detection rate has increased from 65% to 92%; the early warning time for abnormalities has been extended from 24 hours to 96 hours; the false alarm rate has been reduced from 12% to 3%. Improvements in system performance monitoring: the monitoring accuracy of compressor efficiency has increased by 24%; the monitoring accuracy of energy loss in gas storage has increased by 37%; the monitoring of heat transfer efficiency in the heat storage device has increased by 43%; the prediction accuracy of the output of the expander has increased by 21%. In the detection of actual fault cases, the scaling problem of the heat storage device was successfully detected 72 hours in advance, and the area was located in the first 30% of the high-temperature section; a small leakage (0.5% / day) in the gas storage was accurately identified, which took 7 days to discover by traditional methods; the trend of compressor efficiency decay was identified, and maintenance was arranged in advance to avoid further efficiency decline.

[0238] The present invention realizes the precise evaluation of the operating state of the CAES system. First, multi-time scale decomposition is used to process data corresponding to different time scales of the working characteristics of different components in the CAES system (compressors and expanders respond quickly, while gas storage and heat exchangers change slowly), avoiding the problem that short-time scale features are masked by long-time scale data, making the evaluation equally sensitive to the rapidly changing system characteristics and long-term evolution trends. Secondly, the non-linear time-varying adaptive principal component analysis overcomes the limitation of the traditional PCA assuming a linear relationship of data. By adaptively segmenting local linear regions and constructing a mixed kernel function, the non-linear relationship of the CAES system under different working conditions is accurately captured, and the dimensionality reduction effect is improved by about 25%. Thirdly, the PCA-clustering iterative optimization mechanism breaks the paradigm of the traditional one-way execution of PCA and clustering. The two feedback and adjust each other, and the abnormal sensitivity is enhanced by 35%, enabling the detection of early abnormalities such as small leaks and efficiency decay. Finally, the function-distinguishing area evaluation and dynamic weight adjustment based on thermodynamic constraints consider the contribution changes of each functional area to the overall performance under different operating environments of the CAES system. The evaluation results are more objective and accurate, and the overall evaluation accuracy is increased by 28% compared with the traditional method, providing a reliable basis for system operation and maintenance decisions.

[0239] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of 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 belong to the protection scope of the present invention.

Claims

1. A method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis, characterized in that: include: Collect and preprocess the CAES operation data to obtain a preprocessed data set; Perform multi-time scale decomposition and feature extraction on the preprocessed data set to construct a multi-time scale feature set; Perform nonlinear time-varying adaptive principal component analysis on the multi-time scale feature set to extract the main features of the CAES operating status and form a reduced-dimensional feature space and principal component transformation matrix; The PCA-clustering iterative optimization mechanism is used to perform cluster analysis on the reduced-dimensional feature space to identify and cluster the CAES operation modes; Combining the principal component transformation matrix and the CAES operation mode clustering results, the CAES functional area evaluation model based on thermodynamic constraints is adopted to calculate the functional area performance index value and dynamically adjust the weight of each evaluation index to generate the CAES comprehensive operation status evaluation result; Nonlinear time-varying adaptive principal component analysis includes local linear region adaptive segmentation, specifically: Constructing extended phase space features based on multi-time scale feature sets; The kernel density estimation method is used for the extended phase space characteristics to generate phase space density distribution data; Based on the phase space density distribution data, local linearity data is calculated; Combined with local linearity data, determine the initial point of regional segmentation; Based on the initial point of regional segmentation, the recursive binary K-means algorithm is used to adaptively segment the extended phase space features to obtain the segmentation results; Perform boundary optimization and linearity evaluation on the segmentation results to obtain local linear region partition data; Cluster analysis of the reduced-dimensional feature space includes the injection of prior knowledge of abnormal patterns, specifically: Construct an abnormal mode knowledge base containing typical abnormal modes of compressed air energy storage system; Extract abnormal feature vector sets based on abnormal pattern knowledge base and dimensionality reduction feature space; Using the abnormal feature vector set and historical normal operation data, a normal-abnormal boundary model is established; The abnormal feature vector set and the normal-abnormal boundary model are converted into a clustering constraint condition set; and the constraint strength matrix is ​​calculated by combining the reduced dimension feature space; The constraint strength matrix is ​​used to modify the clustering objective function to generate anomaly-sensitive clustering targets and anomaly prior probability distributions.

2. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 1 is characterized in that: Nonlinear time-varying adaptive principal component analysis also includes constructing a mixed kernel function, specifically: Analyze data distribution characteristics for each region in the local linear region partition data to generate regional characteristic data; Based on the regional characteristic data, a suitable kernel function candidate set is selected for each region; Optimize the parameters of the kernel function candidate set to generate optimized kernel function parameters; By optimizing the kernel function parameters, the weight of the hybrid kernel function in each region is determined through quadratic programming optimization; Based on the weight of the hybrid kernel function, a regional smooth transition function is constructed, and the kernel functions of all regions are integrated to generate a hybrid kernel function set.

3. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 2 is characterized in that: Nonlinear time-varying adaptive principal component analysis also includes adaptive reconstruction error minimization, specifically: The kernel principal component analysis is performed on the extended phase space features using a hybrid kernel function set to generate a kernel principal component projection matrix. Based on the kernel principal component projection matrix, the cumulative variance contribution rate is calculated for each region, and the number of principal components required is preliminarily estimated; The reconstruction error data with different numbers of principal components were calculated by cross-validation; Based on the reconstruction error data, the optimal number of principal components for each region is determined using the elbow rule and information criterion; Reconstruct the regional projection matrix according to the optimal number of principal components and calculate the principal component score of each region; Through the smooth transition function of regional membership, the principal component scores of each region are integrated to form a reduced-dimensional feature space and a principal component transformation matrix.

4. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 1 is characterized in that: Cluster analysis of the reduced-dimensional feature space also includes PCA-clustering iterative optimization, specifically: Read the clustering assignment result obtained by clustering the reduced-dimensional feature space; Calculate cluster quality indicators of cluster assignment results; Based on the clustering quality index, clustering assignment results and principal component transformation matrix, the principal component weights are adjusted to generate adjusted PCA weights; Reproject the dimension-reduced feature space using the adjusted PCA weights to generate adjusted dimension-reduced features; Combining the adjusted dimensionality reduction features, the abnormality-sensitive clustering target, and the abnormality prior probability distribution, re-clustering is performed to generate updated clustering results; Determine whether the updated clustering result converges. If so, obtain the optimized clustering result. Otherwise, return to the step of calculating the clustering quality index.

5. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 4 is characterized in that: Cluster analysis of the reduced-dimensional feature space also includes abnormal pattern feature enhancement, specifically: Analyze and optimize the sample distribution of each cluster in the clustering results and generate category distribution statistics; Based on category distribution statistics, optimized clustering results and dimensionality reduction feature space, a normal-abnormal comparison sample pair set is constructed; Construct a contrast loss function based on contrast sample pairs; Using contrast loss function and contrast sample pair set to train pre-configured feature enhancement network to generate feature enhancement model; The feature enhancement model is used to transform the reduced-dimensional feature space to generate CAES abnormal mode features and CAES operation mode clustering results.

6. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 1 is characterized in that: Multi-time scale decomposition and feature extraction include multi-scale reconstruction and separation, specifically: Apply discrete wavelet transform to the preprocessed data set to obtain wavelet decomposition coefficients; Analyze the energy distribution of wavelet decomposition coefficients and determine the time scale demarcation threshold; Based on the time scale demarcation threshold and wavelet decomposition coefficients, the change components are reconstructed separately, including fast change components, medium change components and slow change components; and the time scale coupling matrix between the change components is calculated; Based on the time scale coupling matrix, the coordinated variation patterns of signals at different time scales are identified; Based on the coordinated change pattern, various change components and their relationships are integrated to construct multi-time scale decomposition data.

7. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 6 is characterized in that: Multi-time scale decomposition and feature extraction also includes multi-scale feature fusion, specifically: Based on the multi-time scale decomposition data, the nonlinear feature set is extracted and standardized to obtain the standardized multi-scale features; and it is organized into a third-order feature tensor of sample × feature × time scale; Based on the third-order feature tensor, the optimal parameters of tensor decomposition are determined through cross-validation; Tucker decomposition is performed on the third-order feature tensor to obtain the decomposition result; Analyze the importance of different feature-time scale combinations based on the decomposition results and generate a feature-scale importance matrix; The key feature-scale combinations are selected using the feature-scale importance matrix to construct a streamlined feature map; The fused features are reconstructed based on the simplified feature mapping and decomposition results to generate a multi-time scale feature set.

8. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 1 is characterized in that: The CAES functional area assessment model based on thermodynamic constraints includes the construction of thermodynamic constraint relationships, specifically: Based on the preprocessed data set, the compressed air energy storage system is divided into compression area, energy storage area and energy release area, and a functional area thermodynamic model is established; Based on the first law of thermodynamics, energy balance constraint equations are established for each functional area; Based on the second law of thermodynamics, the entropy generation rate of each functional area is calculated and the entropy balance constraint equation is established; Combined with the gas state equation, the fluid state constraint equations of each functional area are established; Analyze the coupling relationship between functional areas and construct coupling constraints of functional areas; Based on the functional area thermodynamic model and energy conversion theory, efficiency limit constraints are constructed; The energy balance, entropy balance and fluid state constraint equations as well as functional area coupling constraints and efficiency limit constraints are integrated to form a thermodynamic constraint matrix.

Images