Compressed air energy storage system operation state evaluation method based on principal component clustering analysis

Abstract

The invention discloses a compressed air energy storage system operation state evaluation method based on principal component clustering analysis. The method comprises the steps that system operation data are collected and preprocessed; performing multi-time scale decomposition and feature extraction on the preprocessed data set; performing nonlinear time-varying adaptive principal component analysis on the multi-time scale feature set, and extracting main features of a system operation state; carrying out clustering analysis on the dimension reduction feature space by adopting a PCA-clustering iterative optimization mechanism; establishing an evaluation model of each functional area of the system based on thermodynamic constraints, and calculating performance index values of the functional areas in combination with the principal component transformation matrix and an operation mode clustering result; dynamically adjusting the weight of each evaluation index; and generating a CAES comprehensive operation state evaluation result. According to the method, the problems of poor PCA adaptability, difficulty in multi-time scale data fusion, low abnormal mode detection sensitivity and the like under the nonlinear time-varying characteristic of the CAES system are solved.

Description

Technical Field

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

[0002] Compressed air energy storage system (CAES) is an efficient large-scale energy storage technology that plays an important role in renewable energy grid connection, power peak load regulation and energy optimization configuration. Accurately evaluating the operating status of the CAES system is of great significance for ensuring the safe and efficient operation of the system, extending the life of the equipment and reducing the operation and maintenance costs. The CAES system integrates multiple equipment such as compressors, gas storage, heat exchangers and expanders. The operation process involves the conversion of multiple energy forms such as mechanical energy, electrical energy, thermal energy and potential energy. The system has many operating parameters, high data dimensions, complex relationships between parameters, and the operating characteristics of each component correspond to different time scales - the state of the compressor and expander changes quickly (seconds), while the thermodynamic characteristics of the gas storage and heat accumulator change slowly (minutes to hours). Such high-dimensional, multi-time-scale and strong nonlinear characteristics make the evaluation of the operating status of the CAES system face severe challenges, and advanced data analysis and evaluation methods are needed to accurately characterize the system status and promptly identify abnormal conditions.

[0003] At present, the evaluation of the operation status of CAES systems mainly adopts physical model-based methods and data-driven methods. The physical model-based method relies on accurate thermodynamic models to construct energy balance equations, fluid state equations, etc. Although it can reflect the physical nature of the system, the model is highly complex, parameter adjustment is difficult, and it is difficult to adapt to the dynamic characteristics of the system under multiple working conditions. Data-driven methods such as principal component analysis (PCA) and cluster analysis are widely used in industrial process monitoring. Traditional PCA projects high-dimensional data into low-dimensional space through linear transformation for anomaly detection and state classification; while cluster analysis performs unsupervised grouping of 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 CAES system status evaluation models, or use gray correlation analysis, hierarchical analysis method, etc. to determine the weights of evaluation indicators. However, these methods usually perform PCA and clustering as independent steps linearly, lack interactive optimization mechanisms, and have obvious limitations when dealing with nonlinear, time-varying, and multi-time scale data unique to CAES systems.

[0004] In summary, the existing PCA methods have the following main problems in the evaluation of the operating status of CAES systems: First, the standard PCA assumes that the data distribution satisfies a linear relationship, while the CAES system exhibits obvious nonlinear time-varying characteristics under different operating conditions and loads. In particular, when the compressor and expander are in the start-up and shutdown process and the load changes rapidly, the relationship between the parameters is strongly nonlinear, resulting in poor dimensionality reduction effect of the standard PCA during the operating condition conversion and inaccurate principal component extraction. Second, the existing PCA-clustering process uses a unified time window and weight for all data, which cannot effectively handle the multi-time scale data fusion problem unique to the CAES system, resulting in rapidly changing characteristics being submerged by long-time scale data or long-term trends being obscured by short-term fluctuations. Third, PCA dimensionality reduction and clustering analysis are usually performed in a linear sequence, lacking an iterative optimization mechanism. When the system has unknown operating abnormal modes (such as small leaks, heat exchange efficiency attenuation, etc.), these abnormal features may be filtered in the PCA dimensionality reduction or misclassified in the clustering, reducing the sensitivity of abnormality detection. Finally, the existing evaluation method fails to fully consider the thermodynamic constraints and mutual influence of each functional area of ​​the CAES system. It uses a general entropy weight method to determine the indicator weights, ignoring the dynamic contribution of each functional area's performance 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 purpose of the invention is to provide a method for evaluating the operating status of a compressed air energy storage system based on principal component clustering analysis, in order to solve one of the above-mentioned problems existing in the prior art.

[0006] Technical solution, a method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis, comprising the following steps:

[0007] Collect and preprocess the CAES operation data 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 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;

[0010] 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;

[0011] Combining the principal component transformation matrix and the CAES operation mode clustering results, a 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.

[0012] Beneficial effects: The present invention effectively solves the problems of poor PCA adaptability, difficulty in multi-time-scale data fusion and low sensitivity of abnormal detection under the nonlinear time-varying characteristics of CAES system through technologies such as multi-time-scale data decomposition, nonlinear time-varying adaptive principal component analysis, PCA-clustering iterative optimization and functional partition evaluation under thermodynamic constraints. It can accurately capture the operating characteristics of CAES system under different working conditions, especially the state changes in the process of working condition conversion such as compressor start-stop and gas storage filling and deflation, improve the abnormal detection rate by 43%, and improve the prediction accuracy by 28%, providing reliable technical support for the safe operation and preventive maintenance of the system. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 It is a flowchart of the steps of a method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis provided by an embodiment of the present invention.

[0014] Figure 2 It is a flowchart of the steps of adaptive segmentation of local linear regions provided by an embodiment of the present invention.

[0015] Figure 3 It is a flowchart of the steps of constructing a hybrid kernel function provided by an embodiment of the present invention.

[0016] Figure 4 It is a flowchart of the steps of minimizing the adaptive reconstruction error provided by an embodiment of the present invention.

[0017] Figure 5 It is a flowchart of the steps of injecting prior knowledge of abnormal patterns provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0018] In order to enable those skilled in the art to better understand the scheme of the present invention, the technical scheme 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 part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of the present invention.

[0019] It should be noted that in order to clearly show the steps of this application, serial numbers are marked for each step in the specification. These serial numbers are only used for the convenience of explanation and do not limit the order of execution of the steps. In actual operation, according to the technical requirements of the specific implementation scenario, the steps can be executed in a different order than that shown in the specification, and in some cases, parallel processing between steps can be achieved.

[0020] like Figure 1As shown, a method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis comprises the following steps:

[0021] S1. Collecting and preprocessing the operating data of various components of the compressed air energy storage system 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 thermal storage device data. The compressor operating data can be inlet and outlet pressure, temperature, power, speed and vibration; the expander operating data can be inlet and outlet pressure, temperature, power and speed; the gas storage data can be the pressure, temperature and gas storage volume in the reservoir; the thermal storage device data can be inlet and outlet temperature and heat exchange efficiency.

[0023] S2, decompose the preprocessed data set into multiple time scales, extract features of different time scales, and construct a multi-time scale feature set;

[0024] Specifically, the preprocessed data set can be broken down into subsequences at three different levels: short-term, medium-term, and long-term, to separate out possible patterns or fluctuations in 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 subsequence at each time scale, the corresponding features are extracted. For example, for short-term subsequences, local fluctuations, peaks, mutations, and other characteristics are focused on; while for medium- and long-term subsequences, overall trends, cyclical changes, or cumulative effects may be focused on. Capture information at different time levels.

[0025] S3, performing nonlinear time-varying adaptive principal component analysis on the multi-time scale feature set, extracting the main features of the system operation state, and forming a reduced-dimensional feature space and a principal component transformation matrix;

[0026] Specifically, a principal component analysis (PCA) method that can capture nonlinear and dynamic changes is used to "simplify" complex data into a few "main clues" (i.e., principal components) that can accurately reflect the current operating status of the system. In this way, the feature space after dimensionality reduction is used to reduce the computational complexity and retain the most representative key information of the system status.

[0027] S4. Use PCA-clustering iterative optimization mechanism to perform cluster analysis on the reduced dimension feature space, identify the system operation mode, and obtain the operation mode clustering results and abnormal mode characteristics;

[0028] Specifically, clustering is performed in low-dimensional space to identify each typical operating state; the clustering results are then 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 abnormal data that is far away from the main groups, thereby helping to discover abnormal operating states.

[0029] S5. Based on thermodynamic constraints, an evaluation model for each functional area of ​​the system is established. Combined with the principal component transformation matrix and the operation mode clustering results, the performance index values ​​of the functional areas are calculated and the weights of each evaluation index are dynamically adjusted to generate CAES comprehensive operation status evaluation results and decision support information.

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

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

[0032] According to one aspect of the present application, the step of pre-processing comprises:

[0033] S11. Read the compressor operation data, expander operation data, gas storage data and heat storage device data, adopt a stratified sampling strategy, use high-frequency sampling for fast-changing parameters (such as pressure and power), and use low-frequency sampling for slow-changing parameters (such as temperature and heat exchange efficiency) to obtain the original multi-source heterogeneous data set.

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

[0035] S13, read the time series synchronization data set, adopt the anomaly 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, correct the marked outliers through high-order spline interpolation, and obtain the anomaly 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 conservation of energy and conservation of mass), build a physical constraint relationship matrix, calculate the constraint deviation, adjust the data that does not meet the physical constraints, and obtain a physically consistent data set.

[0037] S15. Read the physical consistency data set, and according to the physical characteristics and variation range of different parameters, use an adaptive normalization method to calculate normalization parameters (mean, standard deviation, maximum and minimum values) for different types of parameters, perform normalization processing, and 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 and constructing a multi-time scale feature set include:

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

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

[0041] S23. Based on the time scale characteristic distribution and wavelet decomposition coefficients, the decomposed coefficients are divided into fast-changing components (corresponding to 1-2 level detail coefficients), medium-changing components (corresponding to 3-4 level detail coefficients) and slow-changing components (corresponding to 5 level and above detail coefficients and approximate coefficients), and wavelet reconstruction is performed on each component to obtain multi-time scale decomposition data.

[0042] S24. Nonlinear feature extraction methods (including sample entropy, approximate entropy, and Lyapunov exponent) are used to extract indicators characterizing the nonlinear dynamic characteristics of each time scale for the fast-changing components, medium-changing components, and slow-changing components in the multi-time scale decomposition data, respectively, to obtain a nonlinear feature set.

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

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

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

[0046] Analyze the energy distribution of wavelet decomposition coefficients and determine the time scale demarcation threshold;

[0047] Based on the time scale demarcation threshold and wavelet decomposition coefficient, fast-changing components, medium-changing components and slow-changing components are reconstructed respectively;

[0048] Calculate the time scale coupling matrix between fast-varying components, medium-varying components, and slow-varying components;

[0049] Based on the time scale coupling matrix, the coordinated variation patterns of signals at different time scales are identified;

[0050] Based on the coordinated change pattern, fast-changing components, medium-changing components, slow-changing components and their relationships are integrated to construct multi-time scale decomposition data for subsequent nonlinear feature extraction.

[0051] Specifically, read the time scale characteristic distribution and wavelet decomposition coefficients, calculate the cumulative energy distribution function of the wavelet coefficients at each decomposition level, set the energy threshold (usually 85%, 95% and 99%), determine the time scale boundaries of fast changes, medium changes and slow changes, and obtain the time scale demarcation threshold. Read the wavelet decomposition coefficients and the time scale demarcation threshold, select the wavelet coefficients corresponding to fast changes (usually 1-2 level wavelet detail coefficients), use the adaptive threshold method based on Shannon entropy to filter the noise, retain the main characteristics 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 demarcation threshold, select the wavelet coefficients corresponding to medium changes (usually 3-4 level wavelet detail coefficients), 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 time scale demarcation thresholds, select the wavelet coefficients corresponding to slow changes (5-level and above detail coefficients and approximate coefficients), use the wavelet soft threshold method to reduce the fluctuations in the long-term trend, reconstruct the long-term trend signal through the inverse discrete wavelet transform (IDWT), and obtain the slow-changing component. Read the fast-changing component, medium-changing component, and slow-changing 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-changing component, medium-changing component, slow-changing component, and time scale coupling matrix, use the dynamic time warping (DTW) algorithm to analyze the coordinated change pattern between signals of different time scales, identify the response sequence and coupling relationship of different time scales in the system state transition process, and obtain the coordinated change pattern map. Read the fast-changing component, medium-changing component, slow-changing component, time scale coupling matrix, and coordinated change pattern map, organize the signals of the three time scales and their mutual relationships into a unified data structure, construct a multi-scale signal representation, and form multi-time scale decomposition data.

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

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

[0055] Based on the multi-time scale decomposition data, the nonlinear feature set is extracted and standardized to obtain the standardized multi-scale features;

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

[0057] 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 results, including the core tensor and three factor matrices;

[0058] Analyze the importance of different feature-time scale combinations based on the decomposition results and generate a feature-scale importance matrix;

[0059] The key feature-scale combinations are selected using the feature-scale importance matrix to construct a streamlined feature map;

[0060] Based on the simplified feature mapping and decomposition results, the fusion features are reconstructed to generate a multi-time scale feature set for subsequent nonlinear time-varying adaptive principal component analysis.

[0061] Specifically, read the multi-time scale decomposition data and nonlinear feature set, use the Z-score standardization method for the features of different time scales, calculate the mean and standard deviation of the features of each time scale, and perform standardization processing to ensure that the features of different scales are comparable in subsequent fusion, and obtain standardized multi-scale features. Read the standardized 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, ensure that the tensor structure retains the relationship between time scales, and obtain a 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 the objective function based on reconstruction error and information retention rate, iteratively optimize the decomposition parameters, and 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 the core tensor and three factor matrices, where the three factor matrices correspond to the sample mode, feature mode and time scale mode respectively, and obtain the Tucker decomposition result, including the core tensor and the three factor matrices. Read the Tucker decomposition result, calculate the energy distribution of each element in the core tensor, analyze the contribution of different mode combinations, construct the 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 Tucker decomposition result, set the importance threshold (usually 90% cumulative contribution rate), select the feature-scale combination with importance exceeding the threshold, construct a simplified feature set, and retain the mapping relationship between the original features and the selected features to obtain a simplified feature map. Read the simplified feature map and Tucker decomposition result, reconstruct the fused feature representation based on the selected core tensor elements and the corresponding factor matrix columns, ensure that the reconstructed features retain the key information of the original multi-scale features, and reduce the redundancy, and generate the final multi-time scale feature set.

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

[0064] According to one aspect of the present application, the steps of extracting the main features of the system operation state and forming a dimensionality reduction 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 through the mutual information method, determine the optimal embedding dimension through the false nearest neighbor method, and obtain the extended phase space feature.

[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, a hybrid kernel function adapted to the nonlinear characteristics of each local region is constructed. The radial basis kernel function (RBF) and the polynomial kernel function are combined, and the optimal kernel parameters are determined through cross-validation to obtain a hybrid kernel function set.

[0068] S34. Combining the extended phase space features and the mixed kernel function set, performing kernel principal component analysis on each local region, calculating the kernel matrix, solving the eigenvalue problem, and obtaining the kernel principal component projection matrix of each region.

[0069] S35. Based on the kernel principal component projection matrix, the optimal number of principal components for each local area is determined by the pre-optimization method. The number of principal components is adaptively adjusted using the principle of minimizing the reconstruction error. The principal component scores of each area are calculated and integrated to form a reduced-dimensional feature space and a principal component transformation matrix.

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

[0071] Constructing extended phase space features based on multi-time scale feature sets;

[0072] The kernel density estimation method is used for the extended phase space characteristics to generate phase space density distribution data;

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

[0074] Combining phase space density distribution data and local linearity data, determine the initial point of regional segmentation;

[0075] 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;

[0076] The segmentation results are subjected to boundary optimization and linearity evaluation to obtain local linear region partition data for subsequent hybrid kernel function construction.

[0077] Specifically, the extended phase space features are read, the kernel density estimation (KDE) method is used, and the adaptive bandwidth Gaussian kernel function is used to calculate the local density of each point in the phase space, construct a density distribution map, identify high-density areas and low-density areas in the phase space, and obtain the phase space density map. The extended phase space features and phase space density map are read, and each data point is centered and its K nearest neighbor points are selected (K is determined by cross-validation, usually 10-30), the local covariance matrix of these points is calculated, and eigenvalue decomposition is performed. The linearity of the local area is evaluated by the eigenvalue ratio, and a local linearity map is constructed. The phase space density map and the local linearity map are read, the composite index of density × linearity is calculated, the local extreme value points of the composite index are found, and the peak detection algorithm is applied to determine the seed points for the initial segmentation. These points usually correspond to the centers of the areas with high linearity and dense samples in the phase space, and the initial segmentation seed points are obtained.

[0078] Read the extended phase space features and initial segmentation seed points, and perform a recursive binary K-means clustering process: initialize a region containing all data points; select the two farthest points from the initial segmentation seed points of the current region as the initial clustering 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 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 to be segmented; take the next region from the queue to be segmented, and repeat the above steps until the queue to be segmented is empty or the preset maximum number of regions is reached, and the initial region division result is obtained.

[0079] Read the initial region division results, design a membership function based on Mahalanobis distance for the region boundary points, calculate the membership 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 multivariate linear regression) on the data points in each region, calculate the root mean square error (RMSE) of the fitting residual as an accurate evaluation of the linearity of the region, and obtain the regional linearity score. Read the optimized region division and regional linearity score, analyze the linearity similarity and boundary smoothness of adjacent regions, design the region merging criteria, merge adjacent regions with similar linear characteristics, optimize the overall segmentation results, and output the final local linear region division.

[0080] This embodiment determines the boundary based on the intrinsic structure of the data, and performs more detailed divisions in the CAES system operating condition conversion area (such as the transition from low to high compressor load, and the turning point from inflation to deflation of the gas storage). Experiments have shown that compared with the standard PCA, this embodiment reduces the dimensionality reduction reconstruction error by 37%, especially during periods of rapid changes in system operating conditions, the reconstruction accuracy is improved by 42%. For example, when the gas storage is full (65-70bar high pressure area), the gas compressibility coefficient changes more sharply, and the PVT relationship becomes more nonlinear. This embodiment automatically identifies and finely divides this key area, thereby improving the monitoring accuracy of the gas storage status.

[0081] like Figure 3 As shown, according to one aspect of the present application, the step of obtaining a hybrid kernel function set includes:

[0082] Analyze data distribution characteristics for each region in the local linear region partition data to generate regional characteristic data;

[0083] Based on the regional characteristic data, a suitable kernel function candidate set is selected for each region;

[0084] The kernel function candidate set is optimized by cross-validation method to generate optimized kernel function parameters;

[0085] By optimizing the kernel function parameters, the weight of the hybrid kernel function in each region is determined through quadratic programming optimization;

[0086] Based on the weights 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 for subsequent kernel principal component analysis.

[0087] Specifically, read the local linear region division and extended phase space characteristics, calculate the statistical characteristics of each linear region, including 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 regional characteristic index. Predefine a variety of kernel function primitives, including linear kernel, polynomial kernel (different orders), radial basis kernel (different bandwidth), sigmoid kernel and Laplace kernel, etc., build a kernel function primitive library, set the parameter range for each kernel function, and form a kernel function primitive set. Read the regional characteristic index and kernel function primitive set, and select the kernel function type suitable for the regional characteristics from the kernel function primitive set according to the characteristic index of each region (such as linearity and distribution form), generate a kernel function candidate set for each region, and obtain a regional kernel function candidate set.

[0088] Read the regional kernel function candidate set and the extended phase space features. For each candidate kernel function in each region, use grid search and cross-validation methods to optimize the kernel function parameters (such as the order of the polynomial kernel and the bandwidth σ of the RBF kernel) to minimize the reconstruction error of the data in the region and obtain the optimized kernel function parameter set. Read the regional 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 between reconstruction error and complexity, determine the mixed weights of each basic kernel function through the quadratic programming optimization method, and obtain the mixed kernel weight vector. Read the regional kernel function candidate set, optimize the kernel function parameter set and the mixed kernel weight vector, and construct a specific mixed kernel function for each region based on the optimized parameters and weights: Kmixed(x, y)=∑ i=1 n wi·Ki(x, y); where Ki is the basic kernel function and wi is the corresponding weight, satisfying ∑ i=1 n wi=1, wi≥0, get the optimized hybrid 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 hybrid kernel function of each region, design a smooth transition function between regions, ensure that the kernel function changes smoothly at the region boundary, avoid discontinuity at the region boundary, and integrate to get a complete set of hybrid kernel functions.

[0089] This embodiment enables the kernel principal component analysis (KPCA) to more accurately extract the nonlinear characteristics of the CAES system. For example, for the strong nonlinear region in the startup phase of the compressor, the system automatically increases the RBF kernel weight; for the near-linear region of steady-state operation, the linear kernel weight is increased. Experiments show that compared with KPCA with 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 hybrid kernel function achieves continuous transformation between regions through a smooth transition function, avoiding the discontinuity of the feature representation of the CAES system when switching operating conditions (such as inflation to deflation), and providing a coherent and consistent feature basis for the overall evaluation of the system status.

[0090] like Figure 4 As shown, according to one aspect of the present application, the steps of forming a reduced-dimensional feature space and a principal component transformation matrix include:

[0091] 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.

[0092] 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;

[0093] The reconstruction error data with different numbers of principal components were calculated by cross-validation;

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

[0095] Reconstruct the regional projection matrix according to the optimal number of principal components and calculate the principal component score of each region;

[0096] 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.

[0097] Specifically, read the kernel principal component projection matrix, use the cumulative variance contribution method to calculate the cumulative contribution rate curve of the eigenvalue for each local area, set the initial threshold (usually 85%), determine the minimum number of principal components that meet the threshold, and obtain the initial principal component number estimate. Read the extended phase space features and the initial principal component number estimate, construct a candidate set of the number of principal components for each local area, ranging from [initial estimate -3, initial estimate +5] (make sure it is 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 perform the following operations for each candidate number of principal components in each local area: on the training set, use the first k principal components to project the data to obtain the reduced-dimensional data; use the same k principal components to reconstruct the reduced-dimensional data back to the original space; 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, 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 optimal number of regional principal components. Read the kernel principal component projection matrix and the optimal number of regional principal components, reconstruct the dimension reduction projection matrix for each region based on the determined optimal number of principal components, retain only the important principal components, and obtain the optimized regional projection matrix. Read the extended phase space features and regional projection matrix, transform the data points in each local region using the corresponding optimized projection matrix, calculate the principal component score after dimension reduction, and obtain the regional principal component score of each region. Read the local linear region division, regional principal component score and regional projection matrix, design a smooth transition function based on regional membership, integrate the principal component scores of each region, construct a unified global feature representation, save the projection matrix of each region and its applicable range, and output the dimension reduction feature space and principal component transformation matrix.

[0099] This embodiment configures different operating areas of the CAES system differently - retaining more principal components in areas with complex operating condition changes (such as start-stop processes) and selecting fewer principal components in stable operating areas. Experiments show that this embodiment reduces the reconstruction error by an average of 25%, while increasing the average dimensionality reduction rate by 15%. For example, during the gas storage filling process, the system automatically determines that 7 principal components are required to capture complex thermodynamic changes; for steady-state operation, only 5 principal components are required. Through the reconstruction of the regional projection matrix and the smooth integration of the regional principal component scores, the problem of regional boundary discontinuity in the traditional regionalized PCA method is solved, providing a high-quality low-dimensional feature representation for CAES system state assessment.

[0100] According to one aspect of the present application, the step of performing cluster analysis to obtain operation mode clustering results includes:

[0101] S41. Read the reduced-dimensional feature space, use the multi-density DBSCAN algorithm to perform initial clustering on the reduced-dimensional features, and automatically determine the optimal density parameters of each region through the density adaptability method to obtain the initial clustering results.

[0102] S42, analyzing the boundary points and noise points in the initial clustering results, calculating the Mahalanobis distances between these points and the centers of each cluster, constructing a fuzzy membership function, and performing enhanced processing on the boundary points to obtain a boundary-enhanced clustering result.

[0103] S43. Based on the prior knowledge of typical abnormal patterns of CAES systems, an abnormal pattern feature library is constructed, and this prior knowledge is injected into the clustering process through constraints, and the clustering objective function is adjusted to obtain abnormality-sensitive clustering targets.

[0104] S44, read the boundary enhanced clustering results and abnormality sensitive clustering targets, build an iterative optimization mechanism for PCA and clustering, in each iteration, adjust the PCA weight matrix according to the current clustering results, and then optimize the clustering boundaries according to the new PCA results. After multiple iterations, the optimized clustering results are obtained.

[0105] S45. Based on the optimized clustering results, the feature vector of the abnormal cluster is extracted, and the contrastive learning method is used to enhance the separation between the abnormal mode and the normal mode. An abnormal-normal contrast loss function is constructed, and the enhanced abnormal features are obtained by minimizing the loss function.

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

[0107] like Figure 5 As shown, according to one aspect of the present application, the step of obtaining an abnormally sensitive clustering target includes:

[0108] Construct an abnormal mode knowledge base containing typical abnormal modes of compressed air energy storage system;

[0109] Extract abnormal feature vector sets based on abnormal pattern knowledge base and dimensionality reduction feature space;

[0110] Using the abnormal feature vector set and historical normal operation data, a normal-abnormal boundary model is established;

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

[0112] Based on the reduced dimension feature space and the clustering constraint condition set, the constraint strength matrix is ​​calculated;

[0113] The constraint strength matrix is ​​used to modify the clustering objective function to generate anomaly-sensitive clustering targets and anomaly prior probability distributions for subsequent iterative clustering analysis.

[0114] Specifically, based on the engineering experience and historical operation data of the CAES system, typical abnormal operation modes are identified, including compressor abnormalities (such as abnormal vibration and reduced efficiency), expander abnormalities (such as blade damage and bearing wear), gas storage abnormalities (such as leakage and temperature abnormalities), and heat exchanger abnormalities (such as scaling and reduced heat exchange efficiency). Feature descriptions are established for each abnormal mode, and an abnormal mode knowledge base is constructed. The abnormal mode knowledge base and historical reduced-dimensional feature space data are read, and for each known abnormal mode, its typical feature vector in the reduced-dimensional feature space is extracted, and the feature statistical distribution (mean, variance, skewness, etc.) is calculated to form a mathematical expression of the abnormal mode and obtain an abnormal feature vector set. The abnormal feature vector set and the reduced-dimensional feature space data of historical normal operation are read, and the one-class support vector machine (One-Class SVM) or isolation forest (IsolationForest) algorithm is applied to learn the boundaries of the normal operation area, while considering the distribution of known abnormal modes, and constructing a normal-abnormal boundary model.

[0115] Read the abnormal feature vector set and the normal-abnormal boundary model, and transform the abnormal pattern knowledge into the constraints of the clustering process, including: Must-Link: samples with the same abnormal pattern should be clustered together; Cannot-Link: samples with different abnormal patterns should not be clustered together; Abnormal sample identification constraints: samples close to known abnormal feature vectors should be marked as potential anomalies; Normal area constraints: samples located in the normal area tend to be clustered into normal classes; Combining these constraints, form an abnormal constraint set. Read the abnormal constraint set and the current dimensionality reduction feature space, analyze the similarity between the current data and the historical abnormal pattern, design the distance metric function, calculate the distance between the current data point and each abnormal feature vector, adaptively adjust the constraint strength based on the distance, and generate a constraint strength matrix. Read the abnormal constraint set and the constraint strength matrix, modify the objective function of the standard clustering algorithm, and add the constraint term: J=J original +λJ constraint Among them, 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, abnormality-sensitive clustering targets and principal component transformation matrix, set the iteration termination conditions (such as the maximum number of iterations, the clustering result change rate threshold), initialize the iteration counter and the intermediate result storage space, and establish the iterative optimization framework. Read the current clustering assignment results (initially the boundary-enhanced clustering results), calculate the clustering quality evaluation indicators, including the Davies-Bouldin index, the silhouette coefficient and the specificity index (considering the precision and recall rate of abnormality detection), and obtain the clustering quality index. Read the current clustering assignment results, clustering quality indicators and principal component transformation matrix, analyze the separation between different clusters, identify clusters with blurred or overlapping boundaries, calculate the inter-class and intra-class scatter matrices based on the Fisher discriminant analysis principle, adjust the weights of each principal component in the PCA transformation matrix, enhance the principal component weights that contribute to distinguishing different clusters, and obtain the adjusted PCA weights.

[0125] Read the extended phase space features and adjusted PCA weights, apply the adjusted weights to perform PCA transformation, reproject the data into the feature space, and generate adjusted dimensionality reduction features. Read the adjusted dimensionality reduction features, abnormality-sensitive clustering targets, and abnormality prior probability distributions, apply a semi-supervised constrained clustering algorithm (such as COP-KMeans or MPCK-Means), perform re-clustering, and obtain updated clustering results. Read the current updated clustering results and the clustering assignment results of the previous round, calculate the rate of change of the two clustering results (such as Rand index, mutual information), and check whether the convergence conditions are met. If not converged and the maximum number of iterations has not been reached, the updated clustering results are used as the new clustering assignment results, and the calculation of the clustering quality index is returned to continue iterating. If converged or the maximum number of iterations has been reached, the final optimized clustering results are output. Read all updated clustering results in multiple rounds of iterations, apply ensemble learning methods (such as Majority Voting or label propagation), integrate the clustering results of multiple rounds of iterations, improve clustering stability and reliability, and generate the final optimized clustering results.

[0126] This embodiment improves the accuracy of state classification, especially in the operating condition conversion area and near the abnormal boundary, and the clustering accuracy is improved by 31%. For example, in the process of the compressor switching from medium load to high load, this embodiment can accurately identify whether it is a normal parameter fluctuation caused by load change or a parameter offset caused by equipment abnormality. Iterative convergence is usually completed within 3-5 rounds, and the computational overhead is controllable. In collaboration with the injection of prior knowledge of abnormal patterns, it not only improves the accuracy of CAES system state classification, but also enhances the detection sensitivity of unknown anomalies, providing more accurate operating 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 abnormal mode characteristics include:

[0128] Analyze and optimize the sample distribution of each cluster in the clustering results and generate category distribution statistics;

[0129] Based on category distribution statistics, optimized clustering results and dimensionality reduction feature space, a normal-abnormal comparison sample pair set is constructed;

[0130] Construct a contrast loss function based on contrast sample pairs;

[0131] Using contrast loss function and contrast sample pair set to train pre-configured feature enhancement network to generate feature enhancement model;

[0132] Identify the hard-to-distinguish samples in the comparison sample pair set, fine-tune the feature enhancement model, and obtain the optimized feature enhancement model;

[0133] The reduced-dimensional feature space is transformed using an optimized feature enhancement model to generate abnormal pattern features and operating mode clustering results for subsequent functional zoning evaluation under thermodynamic constraints.

[0134] Specifically, read the optimized clustering results, count the number of samples in each cluster, calculate the sample ratio of normal clusters to abnormal clusters, detect the degree of category imbalance, provide a reference for subsequent feature enhancement, and obtain category distribution statistics. Read the optimized clustering results, dimension reduction feature space, and category distribution statistics, select samples from normal clusters and abnormal clusters to construct contrast sample pairs, adopt a stratified sampling strategy to ensure that the sample pairs cover different types of abnormal patterns and normal states, and generate a contrast sample pair set. Based on the contrast sample pair set, design the loss function of 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 ; Among them, f i and f j is the sample feature, y ij is the sample pair label (1 for the same class, 0 for different classes), d is the distance function, α is the boundary parameter, and the initial α value is set (usually half of the average distance between classes) to obtain the contrast loss function.

[0135] Read the dimension of the reduced feature space, design the feature enhancement network architecture, use the multi-layer perceptron (MLP) structure, including 2-3 hidden layers, each layer uses batch normalization and LeakyReLU activation function, the output layer dimension is the same as the input feature, initialize the network parameters, and build the feature enhancement network. Read the contrast sample pair set, contrast loss function and feature enhancement network, set the training parameters (learning rate, batch size, number of training rounds), use the small batch stochastic gradient descent optimization algorithm, train the feature enhancement network, minimize the contrast loss, and obtain the trained feature enhancement model. Read the trained feature enhancement model and the contrast sample pair set, identify the "hard sample pairs" (sample pairs with high loss values) that are difficult to distinguish, increase the weight 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 reduced feature space, transform all samples through the feature enhancement model, generate enhanced feature representations, calculate the inter-class separation and intra-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 unbalanced data characteristics of anomaly detection in the CAES system. Experiments have shown that after feature enhancement, the separation between abnormal samples and normal samples in the feature space is increased by an average of 56%, and the purity of abnormal clustering is increased by 38%. For example, for anomalies with weak signals such as small leaks in gas storage reservoirs, 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 reduced heat transfer efficiency of the heat storage device, discovered the problem 3 weeks earlier than traditional methods, and avoided further deterioration of system efficiency. By identifying and focusing on "difficult samples", this embodiment further enhances the sensitivity to boundary areas and early anomalies, and improves the preventive maintenance capabilities of the CAES system.

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

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

[0139] S52. Combining the functional zone thermodynamic model with the first and second laws of thermodynamics, the constraint relationship between the functional zones in the system is constructed, 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 of 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 a functional area evaluation indicator system.

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

[0142] S55. Combine the operation mode clustering results, thermodynamic-principal component mapping model and functional area evaluation index system to calculate the performance index value of each functional area under different operation modes, analyze the change law of each functional area performance with the operation mode, and obtain the functional area performance index value.

[0143] According to one aspect of the present application, the step of forming a thermodynamic constraint matrix comprises:

[0144] 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;

[0145] Based on the first law of thermodynamics, energy balance constraint equations are established for each functional area;

[0146] 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;

[0147] Combined with the gas state equation, the fluid state constraint equations of each functional area are established;

[0148] Analyze the coupling relationship between functional areas and construct coupling constraints of functional areas;

[0149] Based on the functional area thermodynamic model and energy conversion theory, efficiency limit constraints are constructed;

[0150] The energy balance, entropy balance and fluid state constraint equations as well as the functional area coupling constraints and efficiency limit constraints are integrated to form a thermodynamic constraint matrix for subsequent functional area performance indicator definition and evaluation.

[0151] Specifically, the functional area thermodynamic model is read, and based on the first law of thermodynamics, energy balance equations are established for the compression area, energy storage area, and energy release area of ​​the CAES system, respectively. Considering energy input, output, storage, and loss, the energy flow relationship is calculated for each functional area, including the conversion of mechanical energy, electrical energy, thermal energy, and pressure energy, and the energy balance constraint equation is obtained. The functional area thermodynamic model is read, and based on the second law of thermodynamics, the entropy generation rate of each functional area of ​​the CAES system is calculated, and the entropy increase caused by irreversible processes (such as friction, heat conduction, mixing, etc.) is analyzed, and the entropy balance equation is established. Considering the entropy flow and entropy generation in the system, the entropy balance constraint equation is obtained. The pressure, temperature, volume and other parameters in the functional area thermodynamic model are read, and based on the ideal gas state equation and the actual gas correction equation (such as the van der Waals equation), the non-ideal behavior of air under high pressure conditions is considered, and the pressure-volume-temperature (PVT) relationship equation of air in the CAES system is established to obtain the fluid state constraint equation. The energy balance constraint equation, entropy balance constraint equation and fluid state constraint equation are read, and the coupling relationship between the functional areas of the CAES system is analyzed, including the mass flow and energy flow transfer between the compression area and the energy storage area, and the energy release process between the energy storage area and the energy release area. The boundary condition equations connecting different functional areas are established to obtain the functional area coupling constraints.

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

[0153] The constraint framework based on physical principles in this embodiment ensures that the evaluation results comply with the basic thermodynamic laws of the CAES system. Experiments have shown that after the introduction of thermodynamic constraints, the physical rationality of the functional area performance evaluation is improved by 65%, and the evaluation reliability under extreme conditions is improved by 47%. For example, this embodiment can identify abnormal conditions where the surface data appears normal but violates the laws of thermodynamics. For example, in an actual 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. This type of anomaly is often ignored in traditional methods, but it is successfully captured through thermodynamic constraints, and the anomaly detection rate is increased by 42%. Not only does it improve the accuracy of CAES system evaluation, but it also enhances the interpretability of the results.

[0154] According to one aspect of the present application, the step of obtaining a thermodynamic-principal component mapping model comprises:

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

[0156] Decomposing the thermodynamic evaluation index into the functional expression of basic physical quantities, and obtaining the index decomposition expression;

[0157] Based on the physical-feature mapping relationship and the index decomposition expression, a candidate set of nonlinear regression models from principal components to thermodynamic indexes is constructed;

[0158] Analyze the data distribution in the dimension-reduced feature space and determine the segmented modeling scheme;

[0159] Using the historical functional area performance index values, the parameters of the candidate set of nonlinear regression models are optimized to obtain the optimized regression model; and its performance on the test data is verified to generate a verified regression model;

[0160] The validated regression model, physical-feature mapping relationship and segmented modeling scheme were integrated to construct a thermodynamic-principal component mapping model for calculating the performance index values ​​of functional areas.

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

[0162] Read the regression model candidate set and the dimensionality reduction feature space, analyze the data distribution in the principal component space, determine the boundaries of the area that needs segmented modeling, automatically divide the modeling area through clustering or decision tree methods, select suitable regression models for different areas, and obtain segmented modeling schemes. Read the regression model candidate set, segmented modeling schemes, and historical functional area performance index values, use cross-validation methods and grid search strategies to optimize the parameters of each regression model, minimize the prediction error, determine the optimal model parameters, and obtain the optimized regression model. Read the optimized regression model and test data set, evaluate the prediction performance of the model on unseen data, and calculate the prediction error and determination coefficient (R 2) values ​​and mean absolute percentage error (MAPE) are used to test the generalization ability of the model, and the models with poor performance are adjusted to obtain the verified regression model. The verified regression model, physical-feature mapping relationship and segmented modeling scheme are read, all regression models are integrated into a unified mapping framework, and a complete mapping relationship from the principal component space to the thermodynamic indicator space is established, including mapping functions, applicable conditions and error estimates, and the thermodynamic-principal component mapping model is output.

[0163] This embodiment can directly calculate the thermodynamic performance indicators of the CAES system from the principal component space, such as compression isentropic efficiency, heat exchanger efficiency, etc., avoiding the complex process of first reconstructing the original physical quantity and then calculating the indicator, and improving the calculation efficiency by 65% ​​and the stability of the indicator calculation by 41%. Experimental verification shows that the thermodynamic indicator values ​​obtained through this mapping are 92% consistent with the results calculated directly based on the physical formula, but the noise resistance is improved by 3.5 times. The segmented modeling strategy takes into account the changes in the nonlinear characteristics of the CAES system under different operating conditions, such as using a more complex nonlinear model under the high-pressure state of the gas storage reservoir, so that the indicator calculation maintains high accuracy in the entire operating range, and the maximum error is controlled within 3.8%, providing accurate performance measurement for system operation optimization.

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

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

[0166] S57, read the operation mode clustering results and the functional area performance index values, 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-indicator correlation matrix.

[0167] S58. Read the abnormal pattern characteristics and the functional area performance index value, analyze the sensitivity of each performance index to the abnormal pattern, use the information gain method to calculate the contribution of each index to abnormal identification, and obtain the index abnormal sensitivity.

[0168] S59. Combining the initial importance of indicators, the mode-indicator correlation matrix and the abnormal sensitivity of indicators, a multi-objective optimization problem is constructed. Considering the importance of indicators, the ability to distinguish operating modes and the sensitivity of abnormal detection, a genetic algorithm is used to solve the optimal weight combination and 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, a dynamic weight adjustment mechanism is established, an adaptive adjustment function is designed, and the reference weight vector is adjusted in real time to obtain a dynamic adjustment weight model.

[0170] According to one aspect of the present application, the step of obtaining a reference weight vector includes:

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

[0172] Combining the operation mode clustering results and the functional area performance index values, a mode-index correlation matrix is ​​constructed;

[0173] Calculate the sensitivity of each evaluation indicator to the abnormal pattern and generate the indicator abnormal sensitivity;

[0174] Based on the initial importance of indicators, the pattern-indicator correlation matrix and the abnormal sensitivity of indicators, a multi-objective optimization function including importance objectives, pattern distinction objectives and abnormal sensitivity objectives is constructed;

[0175] Define weight optimization constraints, including non-negative weight constraints, weight sum equals 1 constraints, and weight threshold constraints;

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

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

[0178] Specifically, read the initial importance of the indicator, the pattern-indicator correlation matrix and the abnormal sensitivity of the indicator, and construct three objective functions: Importance objective: f1(w)=∑ i=1 n w i I i , where I i is the indicator importance; mode differentiation target: f2(w)=∑ i=1 n w i C i , where C i is the correlation between the indicator and the operating mode; abnormally sensitive target: f3(w)=∑ i=1 n w i Si, where S i is the abnormal sensitivity of the indicator; combining these three objectives, we get the multi-objective optimization function. Based on the characteristics of the weight vector, we define the constraints of the optimization problem, including: non-negative weight constraint: w i ≥0, any i; weight sum is 1 constraint: ∑ i=1 n w i =1; minimum weight threshold constraint: w i ≥w min(Ensure that each indicator has a minimum impact); Maximum weight threshold constraint: w i ≤w max (Avoid single indicator dominance); obtain optimization constraints.

[0179] Read the multi-objective optimization function and optimization constraints, 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 cluster analysis on the solutions on the Pareto frontier, use the K-means algorithm to divide the solution set into different cluster groups, analyze the characteristics and weight distribution pattern of each group of solutions, and obtain the solution set clustering results. Read the solution set clustering results, select representative solutions (such as cluster centers) from each cluster to form a typical weight scheme set, analyze the weight distribution characteristics and optimization target performance of these schemes, and obtain a typical weight scheme. Read the typical weight scheme and historical performance evaluation data, evaluate the performance of different weight schemes in historical scenarios through the backtesting method, calculate indicators such as evaluation accuracy, stability, and generalization ability, select the weight scheme with the best overall performance, and obtain the optimal weight scheme. Read the optimal weight scheme, perform 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 achieves a balanced configuration of indicator weights, while taking into account the importance of the indicators themselves, the ability to distinguish operating modes, and the sensitivity to anomalies. Experiments show that this embodiment improves the comprehensive performance of CAES system status assessment compared to traditional methods: the assessment accuracy is improved by 28%, and the sensitivity of anomaly detection is improved by 35%. For example, this embodiment can automatically assign appropriate weights to key but easily overlooked indicators such as the pressure fluctuation rate of the gas storage reservoir and the temperature difference non-uniformity of the heat exchanger. Through weight sensitivity analysis, the key indicators that have the greatest impact on the assessment results are also identified, making the CAES system operating status assessment more comprehensive, accurate, and reliable.

[0181] According to one aspect of the present application, the step of obtaining a dynamically adjusted weight model includes:

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

[0183] Based on historical environmental factor indicators and corresponding optimal weight data, an environment-weight relationship model is constructed;

[0184] Combine environmental factor indicators and operation mode clustering results to identify the current system operation 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] Using historical environmental data to verify the adaptive adjustment function, optimize the adjustment parameters, and obtain the optimized adjustment function;

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

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

[0189] Specifically, read the environmental parameters (such as ambient temperature, pressure, humidity) and system operating parameters (such as load level, operating time, start-stop frequency) in the preprocessed data set, identify the key environmental factors that affect system performance, design quantitative indicators, and obtain environmental factor indicators. Read the historical environmental factor indicators and the optimal weight data of the corresponding period, analyze the relationship between environmental factors and weight adjustment, use multivariate regression analysis or decision tree methods to build a mapping model from environmental factors to weight adjustment, and obtain an environment-weight relationship model. Read the current environmental factor indicators and operating mode clustering results, combine historical data, identify the current operating condition type of the system, map the current state to a predefined typical scenario, and obtain the current operating condition classification.

[0190] Read the baseline weight vector, environment-weight relationship model and current working condition classification, and design the 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 base weight, E j is an environmental factor, j is the adjustment function, α j is the influence coefficient, by optimizing α j and f j, and obtain the adaptive adjustment function. Read the historical environmental factor indicators and the corresponding optimal weights, use the adaptive adjustment function to generate the predicted weights, calculate the deviation from the actual optimal weights, optimize the adjustment function parameters by the least squares method or gradient descent method, minimize the predicted deviation, and obtain the optimized adjustment function. Integrate the benchmark weight vector and the optimized adjustment function, build a real-time weight calculation engine, dynamically calculate the optimal weights according to the current environmental factors and system conditions, ensure the smoothness and adaptability of the weight adjustment, and obtain the real-time weight calculator. Read the current weights and system feedback data generated by the real-time weight calculator, design the weight update strategy, including update frequency, smooth transition mechanism, and anomaly detection, to ensure the stability and effectiveness of the 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 indicator is increased; under high load conditions, more attention is paid to 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 stability of system evaluation, and reduces the evaluation fluctuations caused by environmental changes by 65%. In the actual application of a CAES power station, it can successfully identify the abnormal heat exchange between the gas storage reservoir and the formation under winter environmental conditions, which provides a basis for the optimization of system operating parameters, so that the system can still maintain high accuracy and reliability evaluation in extreme environments.

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

[0193] S511, read the functional area performance index value and dynamically adjust the weight model, use the weighted summation method to calculate the comprehensive performance score of the system, calculate the score 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 characteristics and system comprehensive performance scores, a multi-dimensional visualization representation of the CAES system operation status is constructed. A radar chart is used to display the scores of each functional area, and a heat map is used to display the distribution of abnormal characteristics to obtain a state visualization view.

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

[0196] S514. When a system abnormality is detected, a Bayesian network is used to construct a fault diagnosis model based on the abnormal pattern characteristics and the functional area performance index values, to infer the most likely location and cause of the abnormality, and to generate a fault diagnosis report.

[0197] S515, based on the comprehensive system performance score, performance trend prediction results and fault diagnosis report, generates decision recommendations for system operation optimization and maintenance based on the predefined decision rule library, and forms decision support information and CAES comprehensive operation status evaluation results.

[0198] This embodiment solves the problem that the operating characteristics of different components of the CAES system correspond to different time scales; solves the limitations of standard PCA on nonlinear data; solves the problem that traditional methods are difficult to identify unknown abnormal patterns; and takes into account the thermodynamic characteristics of the CAES system and the impact of different operating environments. Case 1: A 100MW CAES power station experienced efficiency fluctuations during operation, but the traditional monitoring system was unable to accurately locate the cause. The power station includes multi-stage compressors, underground gas storage, heat storage devices, and multi-stage expanders. The following shows how to apply this method to perform a status assessment on the CAES system and find the root cause of the reduced efficiency.

[0199] First, the operating data of multiple key components of the CAES system are collected: compressor data: including inlet and outlet pressures (6-80 bar), temperature (30-450°C), shaft power (20-60 MW), vibration signal (0-25 mm / s), and intake flow rate (50-200 kg / s); gas storage data: pressure (40-70 bar), temperature (30-60°C), gas storage capacity (800,000-2,000,000 m 3 ); heat storage device data: inlet and outlet temperature (40-450℃), pressure drop (0.1-0.8bar), heat exchange efficiency (85-95%); expander data: inlet and outlet pressures of each level (70-5bar), temperature (450-40℃), output power (60-90MW).

[0200] The system uses a 100Hz sampling rate for high-speed changing parameters (such as compressor vibration) and a 0.1Hz sampling rate 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 used to identify and correct outliers (such as vibration sensor pulse noise).

[0201] According to the characteristics of CAES system, the pre-processed data is decomposed into multiple time scales: fast-changing components (second level): mainly including compressor vibration, power fluctuation, expander speed change, etc., reflecting the transient dynamic characteristics of the equipment; medium-changing components (minute level): including pressure fluctuation, flow change, temperature fluctuation, etc., reflecting the load adjustment response of the system; slow-changing components (hour level): including gas storage temperature change, thermal storage device temperature distribution, system efficiency trend, etc., reflecting the long-term thermodynamic process. The analysis found that there was a significant correlation anomaly between the compressor outlet temperature and the expander inlet temperature on the medium time scale, indicating that there may be problems with the thermal storage system.

[0202] Due to the strong nonlinear relationship between CAES system parameters, especially the changes in 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 pressure state of the gas storage reservoir (inflation / deflation), the operating data is divided into 12 local regions; for the 70%-90% load area data of the gas storage-deflation 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 area, and a mixed kernel function is used for the transition condition; 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 area (explained variance>92%), and 8-10 principal components are retained in the operating condition conversion area (explained variance>95%), to ensure that the slight efficiency changes of the thermal storage system are captured.

[0203] Based on the reduced-dimensional feature space, cluster analysis is performed in combination with prior knowledge of abnormal patterns: prior knowledge of common abnormal patterns of CAES systems is injected: characteristics of compressor blade wear, scaling characteristics of thermal storage devices, leakage characteristics of gas storage reservoirs, etc.; abnormal characteristics of thermal storage system efficiency are enhanced through PCA-clustering iterative optimization; clustering results show that in addition to expected normal operation, low load, high load and other modes, an abnormal cluster is identified, which is highly similar to the efficiency attenuation characteristics of thermal storage devices; the abnormal pattern feature enhancement step further amplifies the uneven temperature distribution characteristics of the thermal storage device, indicating that scaling may have occurred in some areas of the thermal storage material.

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

[0205] A multi-objective weighted optimization function is constructed for different operating environments of the CAES system: the weights of various indicators are dynamically adjusted according to the current system being in the energy storage cycle and the ambient temperature being 35°C (high temperature in summer); the weight of the temperature distribution uniformity indicator of the heat storage device is automatically increased (from 0.08 to 0.15), and the weight of the compressor efficiency indicator 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 a high temperature environment.

[0206] The comprehensive evaluation results generated by the system clearly pointed out that the heat storage device had local scaling problems, mainly concentrated in the first 30% of the high-temperature section; scaling caused the heat exchange efficiency to drop by about 6%, which was the main reason for the overall efficiency fluctuation of the system; it was predicted that without intervention, scaling would extend to the medium-temperature section within the next 45 days, and the system efficiency would further drop by 3%; it was recommended to clean and maintain the high-temperature section of the heat storage device during the next planned shutdown.

[0207] Case 2: Based on the actual operation data of a 100MW CAES power station. The power station includes a multi-stage compressor unit, an underground gas storage, a heat storage device and a multi-stage expansion unit. The system collects the following key parameter data:

[0208] Compressor data: inlet and outlet pressures (6-80 bar), temperature (30-450 ° C), power (20-60 MW), speed (3000-3600 rpm), vibration (0-25 mm / s); gas storage data: internal pressure (40-70 bar), temperature (30-60 ° C), gas storage capacity (800,000-2,000,000 m 3 ); thermal storage device data: inlet and outlet temperature (40-450℃), pressure drop (0.1-0.8bar), heat exchange efficiency (85-95%); expander data: inlet and outlet pressures (70-5bar), temperature (450-40℃), output power (60-90MW), speed (3000-3600rpm); environmental data: ambient temperature (-10-40℃), atmospheric pressure (0.95-1.05bar), humidity (20-90%); data sampling frequency is divided into high-frequency sampling (100Hz, used for fast-changing parameters such as vibration), medium-frequency sampling (1Hz, used for pressure, power, etc.) and low-frequency sampling (0.1Hz, used for temperature, efficiency, etc.) according to the speed of parameter changes. The data collection cycle covers the complete charging and discharging cycle of the power station, including startup, steady-state operation, operating condition conversion and shutdown stages, and a total of 3 months of continuous operation data (about 20 million records) were collected.

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

[0210] Daubechies wavelet (db4) is used to perform 5-level wavelet decomposition on each parameter signal. First, the energy distribution of the wavelet decomposition coefficients is analyzed to determine the time scale demarcation threshold: fast-changing component (second level): corresponding to 1-2 level detail coefficients, energy accounting for about 25%; medium-changing component (minute level): corresponding to 3-4 level detail coefficients, energy accounting for about 35%; slow-changing component (hour level): corresponding to 5 level detail coefficients and approximate coefficients, energy accounting for about 40%; reconstruct the signal components of the three time scales respectively, and analyze the coupling relationship between different time scales by the mutual information method.

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

[0212] Apply kernel density estimation (KDE) to the extended phase space feature X to calculate the local density at 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 a Gaussian kernel function with bandwidth h: Kh(u) = (1 / sqrt(2π·h)) × exp(-u 2 / (2h 2 )); xi is the ith sample point, 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 multidimensional data of the CAES system, the product kernel density estimation is used: ρ(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 dimension 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 neighbors (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 neighbors, T is transposed. 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 ​​in the eigenvalues, respectively. L(x) close to 1 indicates that the region is approximately linear, and close to 0 indicates isotropic nonlinearity. In the high load region of 80%-100% of the compressor, the linearity index is 0.41, indicating a strong nonlinear characteristic; while in the medium load region of 50%-70%, 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 point as the initial segmentation seed point. For CAES data, 14 seed points are identified, including: 2 seed points in the compressor startup phase (low load and high acceleration area); 5 seed points in the steady-state inflation process (medium load area); 2 seed points in the high pressure area of ​​the gas storage reservoir (65-70bar); 2 seed points in the initial deflation period (high pressure expansion section); 3 seed points in the late deflation period (medium and low pressure expansion section); perform recursive binary K-means algorithm: initialize a single area containing all data points; select the two farthest points from the seed points in the current area as the initial clustering centers; perform K-means (K=2) to divide the current area into The sub-regions are divided into two sub-regions; the internal consistency index and linearity of the sub-regions are calculated; if the segmentation condition is met (Davies-Bouldin index < 0.7 and the number of points > minimum threshold 500), the sub-region is added to the queue to be segmented; the next region is taken out from the queue and the above steps are repeated; the queue ends when it is empty or the maximum number of regions (set to 20) is reached; the first segmentation divides the data into the inflation process (accounting for 58%) and the deflation process (accounting for 42%); the second level of segmentation divides the inflation process into the startup segment, the medium load segment and the high load segment; the third level of segmentation is further refined, and finally 14 regions are formed.

[0215] For the regional boundary points, 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 of point x to the ith region; μi is the center of the ith region; Σi is the covariance matrix of the ith region; and m is the total number of regions. The fuzzy C-means (FCM) algorithm is used to optimize the region boundaries and reduce overlaps to form the final local linear region division. In the 75%-80% load transition region of the compressor, after FCM optimization, the proportion of overlapping samples dropped from 18.7% to 6.3%, and the boundaries became clearer.

[0216] Statistical characteristics are calculated 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 nonlinearity; Region 5 (compressor steady-state load): low skewness (Sr = 0.21), close to normal distribution (Kr = 3.12), indicating approximate linearity; Region 8 (gas storage high pressure area): moderate skewness (Sr = 0.76), high kurtosis (Kr = 5.34), indicating moderate nonlinearity. According to the regional characteristics, suitable kernel function candidates are selected 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 function (RBF) kernel: K_rbf(x, y) = exp(-γ·||xy|| 2 ); Sigmoid kernel: K_sig(x, y) = tanh(γ·x T ·y + c); Laplace kernel: K_lap(x, y) = exp(-γ·||xy||1); For areas with low skewness and approximately normal distribution (such as area 5), ​​select linear kernel and low-order polynomial kernel; for areas with high skewness and high kurtosis (such as area 1), select RBF kernel and Laplace 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 in each region, the parameters are optimized through 5-fold cross validation, and the parameter combination with the smallest reconstruction error is selected. The RBF kernel in region 1 (startup segment): the optimal γ=0.05, the reconstruction error is 0.086; the polynomial kernel in region 5 (steady-state load): the optimal γ=0.01, d=2, c=1.5, the reconstruction error is 0.042; the RBF kernel in region 8 (high-voltage region): the optimal γ=0.08, the 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 weight wi is 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 (startup segment): The weight of the radial basis (RBF) kernel in the mixed kernel function is w_rbf = 0.65, the weight of the Laplacian kernel is w_lap = 0.30, and the weight of the polynomial kernel is w_poly = 0.05; Region 5 (steady-state load): The weight of the linear kernel is w_lin = 0.40, w_poly = 0.50, and w_rbf = 0.10; Region 8 (high-pressure area): w_rbf = 0.55, w_poly = 0.35, and the weight of the Sigmoid kernel is w_sig = 0.10.

[0219] To avoid discontinuity at the region boundary, a smooth transition function K_final(x, y) = Σ(r=1 to R) s_r(x)·s_r(y)·K_r_mixed(x, y) is designed; 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 load to high load transition region of the compressor (about 75% load), the membership of sample point x to region 4 is s_4(x)=0.65, and the membership to region 5 is s_5(x)=0.35, ensuring the smooth change of the kernel function at the region boundary.

[0220] Preliminary estimation of the number of regional principal components, including: performing kernel principal component analysis (KPCA) for each local region r: calculating the kernel matrix: Kij = K_r_mixed(xi, xj); centralizing the kernel matrix: K* = K - 1n·K - K·1n + 1n·K·1n, where 1n is an n×n matrix with all elements equal to 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 determine the minimum number of principal components that make CVR reach 85%; Region 1 (start-up phase): 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-voltage area): the cumulative variance contribution rate of the first 7 principal components is 85.6%.

[0221] Using 5-fold cross validation, the reconstruction error is calculated for different numbers of principal components k in each region: the data of region r is randomly divided into 5 parts; for each fold, KPCA is performed using the remaining 4 folds of data, retaining the first k principal components; the reconstruction error on the test fold is calculated: MSE(k) = (1 / n_test) × Σ(i=1 to n_test) ||Φ(xi) - Φ̂k(xi)|| 2 ; where Φ̂k(xi) is the feature vector reconstructed using k principal components; n_test is the total number of samples included in the test set; the average cross-validation reconstruction error of k principal components is calculated. The average reconstruction error of 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] Combine the "elbow rule" and information criteria (AIC, BIC) to determine the optimal number of principal components: 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 consider the inflection point of the reconstruction error curve. Region 1 (startup section): The reconstruction error curve has an obvious inflection point at k=9, the minimum AIC is at k=9, and the number of principal components is finally determined to be 9; Region 5 (steady-state medium load): The reconstruction error curve tends to be flat at k=5, the minimum BIC is at k=5, and the number of principal components is finally determined to be 5; Region 8 (high-voltage area): The inflection point of the reconstruction error is at k=7, and the AIC and BIC are not much different between k=7 and k=8. Weighing the complexity, k=7 is selected. According to the determined optimal number of principal components, the KPCA projection matrix is ​​reconstructed for each region: principal component score zi(j) = Σ(l=1 to n) αj(l)·K*(xi, xl); where zi(j) is the score of sample xi on the jth principal component; αj is the jth eigenvector; K* is the centralized kernel matrix.

[0223] Integrate the principal component scores of each region and 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 sample x; s_r(x) is the smoothed membership of x to region r; zr(x) is the principal component score of x in region r. In the actual application of the CAES system 100MW power station, the original high-dimensional data (>100 dimensions) was reduced to about 25 dimensions through nonlinear time-varying adaptive principal component analysis, while maintaining key information. Compared with traditional PCA, the reconstruction error was reduced by 37%, especially in the operating condition conversion area, the reconstruction accuracy was improved by 42%.

[0224] The multi-density DBSCAN algorithm is applied to the reduced-dimensional 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 transition conditions): larger ε (0.25) and smaller MinPts (10); the initial clustering results identify 8 clusters, including low-load, medium-load, and high-load regions of the compressor, different load regions of the expander, and boundary transition regions.

[0225] Based on CAES system engineering experience and historical cases, an abnormal pattern knowledge base is established, including: compressor abnormalities: abnormal vibration, blade wear, bearing failure, efficiency reduction, etc.; expander abnormalities: abnormal vibration, turbine blade damage, bearing wear, etc.; gas storage abnormalities: small leaks, temperature abnormalities, pressure fluctuation abnormalities, etc.; heat storage device abnormalities: scaling, reduced heat exchange efficiency, flow channel blockage, etc.; each abnormal pattern contains feature descriptions, influencing parameters and typical feature vectors.

[0226] Extract the feature vectors of various anomalies from historical data and construct an anomaly feature library. Take the fouling of thermal storage device as an example: the feature vector v_fouling = [0.32, -0.78, 0.15, 0.65, -0.41, ...]; it represents the typical performance of fouling anomaly in the reduced dimension feature space (such as significant negative offset on the second principal component and positive offset on the fourth principal component). Use one-class support vector machine (One-Class SVM) to construct the normal operating 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 training; parameter ν is set to 0.05 (abnormal proportion estimation); for known anomalies, calculate the distance distribution from the normal boundary and establish the abnormal-normal boundary model.

[0227] Convert abnormal knowledge into clustering constraints: If d(xi, xj) < Δ_ml and xi, xj are close to the same abnormal feature vector, then impose a must-connect constraint to ensure that they are clustered together. Where d is the feature space distance, and Δ_ml is the must-connect threshold (set to 30% of the average distance in the feature space). If xi is close to the center of the normal operating area and xj is close to a certain abnormal feature vector, then impose a no-connect constraint to ensure that they are not clustered together. For samples 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 constraint strength based on the similarity between sample and abnormal features: 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-connect constraint strength; w_cl is the do-not-connect 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 and add the 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 cluster label of sample xi.

[0229] Calculate the prior probability of belonging to various known abnormal modes for each data point: 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 anomaly a; v_a is the feature vector of anomaly a; σ_a is the scaling parameter (set to 0.5 times the average distance within the anomaly class). For example, the prior probability distribution calculated for a 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 clustering result change rate threshold is 1%; the computational complexity is O(n 2), where n is the number of samples. Calculate the clustering quality evaluation index, Davies-Bouldin index DB = (1 / k)·Σ(i=1to k) max(j≠i) {(Si + Sj) / dij}; where k is the number of clusters; Si is the average sample-to-center distance of the i-th cluster; dij is the distance between the i-th and j-th cluster centers. Silhouette coefficient S = (1 / n)·Σ(i=1 to n)(bi - ai) / max(ai, bi); where ai is the average distance between sample i and other samples in the same cluster; bi is the average distance between sample i and the nearest neighbor cluster sample. Specificity index (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 inter-class and intra-class scatter matrices: Intra-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 cth class, μ is the global mean, nc is the number of samples in the cth class, and C is the total number of classes in the clustering result. Calculate the discriminant ability of each principal component: Fi = vi T ·Sb·vi / (vi T·Sw·vi); where vi is the direction vector of the ith principal component. Adjust the weight 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 ith principal component before the weight adjustment. After the first round of iterations, it was found that the third principal component contributed the most to distinguishing the normal / abnormal state of the thermal storage device (F3=2.86), and the weight was adjusted from 0.10 to 0.15; while the first principal component mainly reflects the system load level and contributes less to abnormality identification (F1=0.75), and the weight was reduced from 0.20 to 0.18. Apply the adjusted weights to reproject the data: z_adjusted = Σ(i=1 to k) w_i_new·zi; where: z_adjusted is the adjusted dimensionality reduction feature; w_i_new is the adjusted weight of the ith principal component; zi is the score of the ith principal component; k is the total number of principal components. Using the adjusted features and anomaly-sensitive clustering targets, constrained clustering is performed: the semi-supervised constrained clustering algorithm MPCK-Means is used; initialization is performed in combination with the anomaly prior probability distribution; and the constrained objective function is optimized.

[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 cluster label of sample i in the tth iteration; I is the indicator function; n is the total number of samples. The iteration is terminated when Δ < 0.01 or the maximum number of iterations is reached.

[0233] In the iterative example: Round 1: DB index = 0.58, Silhouette = 0.67, F_anom = 0.75, Δ = 0.18; Round 2: DB index = 0.52, Silhouette = 0.71, F_anom = 0.81, Δ = 0.08; Round 3: DB index = 0.49, Silhouette = 0.73, F_anom = 0.84, Δ = 0.03; Round 4: DB index = 0.48, Silhouette = 0.74, F_anom = 0.85, Δ = 0.009 (convergence). 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: 92.3% of the total samples; slightly abnormal cluster: 4.7% of the total samples; obviously abnormal cluster: 2.1% of the total samples; seriously abnormal cluster: 0.9% of the total samples. There is obviously a class imbalance problem, and the abnormal samples are far less than the normal samples. Construct normal-abnormal contrast sample pairs: select core samples (50 samples closest to the cluster center) from each abnormal cluster; for each abnormal sample, select the 5 closest normal samples; use stratified sampling to ensure that different types of abnormalities are covered; form a contrast sample set containing about 1000 abnormal samples and 5000 normal samples. The contrast sample pair is defined as (xa, xb, yab), where xa, xb are sample pairs; yab=1 indicates the same type of sample pair; yab=0 indicates a different type of sample pair (one normal, one abnormal). Design contrast learning loss function: L_contrastive =Σ(i, j) yij·d(fi, fj) 2 + (1-yij)·max(0, α - d(fi, fj)) 2 ; where fi, fj are sample features; d is the Euclidean distance; α is the boundary parameter (initialized to 0.5, gradually increased to 0.8 with training); yij is the sample pair label. This loss function makes the features of similar samples approach each other, and the features of heterogeneous samples keep at least α away.

[0235] Design a 3-layer perceptron network: input layer: dimension equal to the dimension of the reduced feature (about 25); hidden layer 1: 32 neurons, LeakyReLU activation, batch normalization; hidden layer 2: 16 neurons, LeakyReLU activation, batch normalization; output layer: dimension the same as the input (about 25). The network is initialized as an identity mapping to ensure that the feature space is not overly distorted in the early stage of training. Training parameter settings: batch size: 64; learning rate: 0.001, using Adam optimizer; number of training rounds: 200; regularization: L2 weight decay (0.001); early stopping: stop if the validation loss does not decrease for 10 consecutive rounds. Training process tracking: contrast loss dropped from the initial 0.873 to 0.217; the average distance of similar samples dropped from 0.421 to 0.168; the average distance of heterogeneous samples increased from 0.582 to 0.836. Identify "hard sample pairs" that are difficult to distinguish: Calculate the loss value for each sample pair; Select the sample pairs with the top 20% loss value as hard samples; Increase the weight of these samples (3 times) and fine-tune the model; Perform additional fine-tuning for specific types of abnormal patterns. In this example, the sample of slight fouling of the thermal storage device and the sample of normal operation with high load are difficult to distinguish under standard training (distance 0.31). After fine-tuning with hard samples, the distance increases to 0.65, which improves the distinguishability. Apply the trained feature enhancement network to all samples: z_enhanced = fθ(z); where z is the original dimension reduction feature; fθ is the feature enhancement network with parameter θ; z_enhanced is the enhanced feature. Calculate the inter-class separation and intra-class aggregation before and after enhancement: Inter-class separation: increased from 0.58 to 0.82 (an increase of 41%); Intra-class aggregation: increased from 0.31 to 0.22 (an increase of 29%). In the case of scaling of the heat storage device, the overlap rate between scaling samples and normal samples in the feature space was 28% before enhancement, and dropped to 7% after enhancement, which improved the accuracy of anomaly detection (from 73% to 94%).

[0236] Based on the thermodynamic constraints, the CAES system is divided into compression zone, energy storage zone and energy release zone, and the thermodynamic model of each zone is established and the performance index is calculated. Through the thermodynamic mapping in the principal component space, the dimension reduction features are associated with the thermodynamic indexes to establish an evaluation model. The adaptive weight dynamic adjustment mechanism determines the benchmark weight based on multi-objective optimization, and dynamically adjusts the weight according to the operating environment (such as load level, ambient temperature), making the evaluation results more accurate and reliable.

[0237] Compared with the traditional PCA+ clustering method, the state recognition accuracy of this embodiment is improved from 78% to 94%; the anomaly detection rate is improved from 65% to 92%; the anomaly warning lead time is extended from 24 hours to 96 hours; and the false alarm rate is reduced from 12% to 3%. Improvements in system performance monitoring: The accuracy of compressor efficiency monitoring is improved by 24%; the accuracy of gas storage energy loss monitoring is improved by 37%; the heat transfer efficiency monitoring of the heat storage device is improved by 43%; and the accuracy of expander output prediction is improved by 21%. In the actual fault case detection, the scaling problem of the heat storage device was successfully detected 72 hours in advance and located in the first 30% area of ​​the high-temperature section; the small leak of the gas storage reservoir (0.5% / day) was accurately identified, which would take 7 days to detect with the traditional method; the compressor efficiency attenuation trend was identified, and maintenance was arranged in advance to avoid further efficiency decline.

[0238] The present invention realizes accurate evaluation of the operating status of the CAES system. First, the multi-time scale decomposition processes the data of different time scales corresponding to the working characteristics of different components in the CAES system (the compressor and expander respond quickly, and the gas storage and heat exchanger change slowly), avoiding the problem that the short-time scale characteristics are masked by the long-time scale data, making the evaluation equally sensitive to the rapidly changing system characteristics and long-term evolution trends. Secondly, the nonlinear time-varying adaptive principal component analysis overcomes the limitation of the traditional PCA assumption of linear relationship of data, and accurately captures the nonlinear relationship of the CAES system under different working conditions through local linear region adaptive segmentation and hybrid kernel function construction, and the dimension 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 are mutually fed back and adjusted, and the abnormal sensitivity is enhanced by 35%, which can detect early abnormalities such as small leaks and efficiency attenuation. Finally, the functional area zoning evaluation and dynamic weight adjustment based on thermodynamic constraints take into account 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 improved 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 are described in detail above; however, the present invention is not limited to the specific details in the above embodiments. Within the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all 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, a 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.

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 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; The segmentation results are subjected to boundary optimization and linearity evaluation to obtain local linear region partition data.

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 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.

4. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 3 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.

5. 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 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.

6. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 5 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.

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: 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.

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: 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.

9. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 8, 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.

10. The method for evaluating the operating status of a compressed air energy storage system based on principal component cluster analysis according to claim 1, 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