An optimized operation method for compressed air energy storage based on multi-parameter collaborative control

Abstract

The present invention discloses an optimized operation method for compressed air energy storage based on multi-parameter collaborative control, including multi-source data acquisition and preprocessing, dynamic estimation of multi-dimensional parameters of thermodynamic state, three-dimensional reconstruction and evolution prediction of the temperature field in the gas storage reservoir, system loss analysis and efficiency evaluation based on entropy generation theory, multi-time scale energy flow dynamic tracking and state evaluation, and generation of optimized operation parameter suggestions based on state evaluation. The present invention realizes accurate estimation of non-ideal gas parameters through particle swarm-Bayesian hybrid inference, reconstructs the temperature field in the gas storage reservoir, introduces an inter-layer available energy transfer term to quantify the influence of temperature stratification, improves the system operation efficiency and reliability, prolongs the equipment life, and reduces the operation cost.

Description

Technical Field

[0001] The present invention belongs to the field of CAES, and in particular, relates to an optimized operation method for compressed air energy storage based on multi-parameter collaborative control. Background Art

[0002] As an important large-scale energy storage technology, the compressed air energy storage (CAES) system can compress and store air through renewable energy power and release energy when needed, effectively solving the problem of accommodating intermittent energy sources such as wind energy and solar energy. With the continuous increase in the proportion of renewable energy, how to improve the operation efficiency of the CAES system, reduce energy losses, and extend the equipment life has become a key technical challenge in the energy transformation. As the basis for optimized control, the operation state assessment of the CAES system can monitor the system performance in real time, predict potential risks, and guide optimization strategies, which is of great significance for improving the system economy and reliability. Especially in scenarios such as large-scale applications and long-term peak shaving, accurate state assessment can significantly improve the system conversion efficiency, reduce the energy storage cost, and provide more reliable support for the power grid peak shaving and frequency modulation.

[0003] At present, the operation state assessment methods of the CAES system mainly include five types of technical routes. The first type is the assessment method based on physical models. By establishing system thermodynamics and fluid mechanics models, the ideal gas state equation or a simply modified state equation is used to calculate the system state parameters. The second type is the assessment method based on data-driven. Historical operation data is used to train statistical models or machine learning models to predict system performance and state changes. The third type is the empirical assessment method. Based on expert experience and operation rules, a rule base and judgment criteria are established to evaluate the system operation state and fault risks. The fourth type is the assessment method based on energy analysis. The energy conversion process and efficiency losses of the system are analyzed through energy balance and the first and second laws of thermodynamics. The fifth type is the hybrid assessment method, which combines physical models and data-driven technologies. It not only uses known physical laws to constrain the model structure but also uses measured data to optimize the model parameters to improve the assessment accuracy and generalization ability.

[0004] The main problems existing in the prior art are as follows: First, in terms of the uncertainty of thermodynamic states, the existing analysis of CAES systems usually assumes compliance with the ideal gas law, while in fact, compressed air exhibits significant non-ideal behavior under high-pressure and variable-temperature environments. Especially during the multi-stage compression process, the dynamic changes in parameters such as specific heat capacity and compression factor are caused by the effects of gas heating, cooling, and phase changes. The existing static models or simplified assumptions cannot accurately capture these changes, resulting in deviations in the prediction of system operating efficiency. Second, regarding the non-uniform temperature distribution in the gas storage reservoir, there are obvious temperature gradients and stratification phenomena inside the gas storage reservoir of large CAES systems (especially underground cavern types), leading to non-uniformity of gas properties in the spatial dimension. The existing control methods regard the gas storage reservoir as a uniform body and only consider the inlet and outlet parameters, ignoring the impact of the internal temperature distribution on system efficiency. Research shows that non-uniform temperature in the gas storage reservoir can lead to an increase in thermodynamic efficiency losses, and there is a lack of charge-discharge optimization strategies based on the real-time distribution of the internal temperature field. Summary of the Invention

[0005] The object of the invention is to provide an optimized operation method for compressed air energy storage based on multi-parameter collaborative control, in order to solve at least one technical problem existing in the prior art.

[0006] Technical solution: An optimized operation method for compressed air energy storage based on multi-parameter collaborative control includes the following steps:

[0007] Collect multi-source sensing data of CAES and process it to generate a CAES dataset;

[0008] Based on the CAES dataset, combined with particle swarm - Bayesian hybrid inference, estimate the system's thermodynamic parameters to form a working condition adaptable thermodynamic model;

[0009] Based on the CAES dataset, reconstruct the temperature distribution inside the gas storage reservoir through a thermal-neural network architecture to obtain a high-precision temperature field model;

[0010] Combine the CAES dataset, the working condition adaptable thermodynamic model, and the high-precision temperature field model to obtain a bottleneck analysis map;

[0011] Based on the CAES dataset and the working condition adaptable thermodynamic model, generate a classification result of the CAES operating state;

[0012] Combine the bottleneck analysis map and the CAES operating state classification result to generate a CAES operation evaluation and optimization report.

[0013] Advantageous effects: The present invention constructs a complete technical system for evaluating the operating state of a compressed air energy storage system, overcomes problems such as inaccurate parameter estimation, insufficient temperature field analysis, incomplete efficiency loss analysis, and lack of coordination of multi-time scale characteristics in traditional evaluation methods, and improves the overall operating efficiency of the system. Description of the Drawings

[0014] Figure 1 It is a flowchart of the steps of an optimized operation method for compressed air energy storage based on multi-parameter collaborative control provided by an embodiment of the present invention.

[0015] Figure 2 It is a flowchart of the steps for estimating the thermodynamic parameters of the system provided by an embodiment of the present invention.

[0016] Figure 3 It is a flowchart of the steps for constructing a parameter-correlated Bayesian network provided by an embodiment of the present invention.

[0017] Figure 4 It is a flowchart of the steps for reconstructing the temperature distribution in the gas storage reservoir provided by an embodiment of the present invention. Detailed implementation manners

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

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

[0020] As Figure 1 shown, an optimized operation method for compressed air energy storage based on multi-parameter collaborative control includes the following steps:

[0021] S1. Collect multi-source sensing data of CAES and process it to generate a complete preprocessed data set, namely the CAES data set;

[0022] Specifically, the multi-source sensing data of CAES includes compressor unit parameter data, gas storage reservoir parameter data, expander unit parameter data, and auxiliary system parameter data. Among them, the compressor unit parameter data can be the inlet and outlet temperature, pressure, flow rate, shaft power, etc. of a multi-stage compressor; the gas storage reservoir parameter data can be the temperature, pressure, humidity, etc. of the spatial distribution at multiple points; the expander unit parameter data can be the inlet and outlet temperature, pressure, flow rate, output power, etc.; the auxiliary system parameter data can be the cooling water temperature, flow rate, heat exchanger temperature, etc.

[0023] S2. Based on the complete preprocessed dataset, combine the particle swarm - Bayesian hybrid inference to estimate the thermodynamic parameters of the system, and form a thermodynamic model adaptable to operating conditions;

[0024] Specifically, particle swarm optimization is similar to a swarm search strategy and is good at finding the optimal solution, while Bayesian inference can handle uncertainty and estimate the possible distribution of parameters. Combining them can deduce and estimate the thermodynamic parameters of the system (such as temperature, pressure, efficiency, etc.) and establish a thermodynamic model that can adapt to different operating conditions (operating modes). That is to say, this model can be adjusted according to the actual operating situation to maintain accuracy and reliability.

[0025] S3. Reconstruct the temperature distribution in the gas storage reservoir through a thermal - neural network architecture based on the complete preprocessed dataset to obtain a high - precision temperature field model;

[0026] Specifically, the thermal - neural network architecture is a model that combines thermodynamic knowledge and the capabilities of neural networks. This model can more accurately understand and predict the distribution of heat; using the model to generate an accurate temperature distribution map of the gas storage reservoir helps to better understand the internal thermal characteristics.

[0027] S4. Combine the complete preprocessed dataset, the thermodynamic model adaptable to operating conditions, and the high - precision temperature field model to obtain a bottleneck analysis map;

[0028] Specifically, by combining the complete preprocessed dataset, the thermodynamic model adaptable to operating conditions, and the high - precision temperature field model, generate a map that can visually display the key bottlenecks in the system, such as areas with limited performance or parts that need to be optimized.

[0029] S5. Based on the complete preprocessed dataset and the thermodynamic model adaptable to operating conditions, generate the classification results of the CAES operating state;

[0030] Specifically, by inputting the preprocessed data into the thermodynamic model adaptable to operating conditions, the performance of the CAES system under various operating states can be analyzed. Then, according to the set classification criteria (such as "efficient", "stable", "inefficient"), these operating states are classified. This can help identify the system characteristics under different states and facilitate optimization and improvement.

[0031] S6. Combine the bottleneck analysis map and the classification results of the CAES operating state to generate a CAES operation optimization recommendation report.

[0032] Specifically, by combining bottleneck analysis and operating state classification, specific optimization recommendations can be formulated. The goal of the report is to make the overall efficiency of the system higher by improving the performance at the bottleneck and adjusting the operating state.

[0033] This embodiment can improve the operation efficiency and reliability of the CAES system throughout the entire process from data collection, model construction to the generation of optimization suggestions; through high-precision models and intelligent analysis methods, accurately locate problems and formulate targeted solutions, thereby reducing energy consumption, improving energy storage efficiency, and increasing system stability.

[0034] According to one aspect of the present application, the steps of processing to generate a CAES data set include:

[0035] S11. Hierarchical and multi-source sensing data collection: Obtain compressor unit parameter data (including inlet and outlet temperatures, pressures, flow rates, shaft powers, etc. of multi-stage compressors), gas storage reservoir parameter data (including spatially distributed temperatures, pressures, humidities, etc. at multiple points), expander unit parameter data (including inlet and outlet temperatures, pressures, flow rates, output powers, etc.), and auxiliary system parameter data (including cooling water temperature, flow rate, heat exchanger temperature, etc.) to form an original multi-source data set.

[0036] S12. Data quality assessment and outlier processing: Conduct integrity checks, consistency verification, and outlier identification on the original multi-source data set, calculate the deviation degree of each measurement point data through the modified Z-score method, and perform cross-verification in combination with system physical constraint relationships (such as mass conservation and energy conservation) to generate a quality assessment index and a labeled data set.

[0037] S13. Multi-scale data synchronization and fusion preprocessing: Divide the labeled data set into high-frequency data (millisecond level, such as pressure fluctuations), medium-frequency data (second level, such as temperature changes), and low-frequency data (minute level, such as energy accumulation) according to the sampling frequency, extract multi-scale features using wavelet transform, and achieve data synchronization through timestamp alignment to generate a multi-scale synchronized data set.

[0038] S14. Physics-constraint-based data completion and reconstruction: For the missing data in the multi-scale synchronized data set, establish a partial differential equation system in combination with the system physical constraint relationship, and perform data completion through a high-order physical interpolation method. For data that is difficult to complete through the physical model, use an improved spatio-temporal tensor completion algorithm for reconstruction to obtain a complete preprocessed data set.

[0039] In this embodiment, through hierarchical and multi-source sensing data acquisition and preprocessing, full-condition and full-system data coverage of the compressed air energy storage system is achieved. Moreover, a strategy combining the modified Z-score method and cross-validation of system physical constraints is adopted to improve data quality, and the outlier detection rate is increased by 33%. The multi-scale data synchronization and fusion preprocessing technology overcomes the data synchronization difficulties brought by different sampling frequencies, enabling precise alignment of high-frequency data (millisecond level) and low-frequency data (minute level). The time synchronization error is reduced from the original 3 - 5 seconds to within 0.1 second. Based on the physical constraint data completion and reconstruction technology, the data missing rate is reduced from the original 7 - 9% to below 1%, providing a high-quality data basis for subsequent state assessment.

[0040] According to one aspect of the present application, the steps of forming a thermodynamics model adaptable to working conditions include:

[0041] S21. Dynamic selection of non-ideal gas state equations: Based on the pressure and temperature ranges in the complete preprocessed dataset, dynamically select applicable state equations (including the Peng-Robinson equation, Redlich-Kwong-Soave equation, etc.), construct the switching boundary of the state equation by minimizing the prediction error, and form an adaptive state equation model.

[0042] S22. Real-time estimation of parameters by particle swarm - Bayesian hybrid inference: Using the complete preprocessed dataset and the adaptive state equation model, construct a Bayesian network to represent the probabilistic dependence relationship between parameters, use particle swarm optimization to initialize the prior distribution, and update the posterior distribution in real time through the sequential Monte Carlo method to obtain the probability distribution of thermodynamic parameters.

[0043] S23. Reliability parameter estimation based on uncertainty propagation: Conduct uncertainty quantification and propagation analysis on the probability distribution of thermodynamic parameters, use the multi-level second moment method to trace the propagation path of uncertainty in the system, establish a credibility evaluation index for parameter estimation, and form a reliability model of thermodynamic parameters.

[0044] S24. Multi-condition parameter sensitivity analysis and dynamic adjustment: Under different working conditions (such as full load, partial load, rapid response, etc.), analyze the sensitivity of each parameter in the reliability model of thermodynamic parameters to the system state, use the sensitivity quantification method based on information entropy, establish a dynamic ranking of parameter importance, and generate a parameter sensitivity matrix and a thermodynamics model adaptable to working conditions.

[0045] In this embodiment, the prediction error of the gas state is reduced from 5 - 7% of the traditional method to within 1.5%. The parameter estimation time is shortened from the minute level of the traditional Bayesian method to the second level, the parameter accuracy is increased by 2 - 3 times, and the energy balance closure is increased to more than 98%. By quantifying the uncertainty of parameter estimation, the prediction interval coverage rate is increased from 85% to more than 95%. It provides a theoretical basis for identifying key parameters under different operating scenarios, making the calculation resource allocation more reasonable, the estimation accuracy of important parameters higher, and the real-time performance and accuracy of state evaluation improved.

[0046] As Figure 2 shown, according to one aspect of the present application, the steps for estimating the thermodynamic parameters of the system include:

[0047] Taking the CAES dataset as the input, using the Particle Swarm Optimization (PSO) algorithm to solve the initial parameter space model including the distribution characteristics of thermodynamic parameters and the physical constraint boundaries of parameters, and obtaining the optimized sampling point set;

[0048] Based on the optimized sampling point set, through the parameter - associated Bayesian network that characterizes the conditional probability relationship between parameters, and combined with real - time observation data, update the posterior distribution of thermodynamic parameters;

[0049] Based on the posterior distribution of thermodynamic parameters, quantify the uncertainty of parameter estimation, and generate the probability distribution of thermodynamic parameters, that is, the thermodynamic parameters of the system.

[0050] As Figure 3 shown, according to one aspect of the present application, the steps for constructing the parameter - associated Bayesian network include:

[0051] Using the optimized sampling point set, calculate the mutual dependence degree between the CAES thermodynamic parameters through information entropy and mutual information analysis, and determine the topological structure of the hierarchical Bayesian network;

[0052] Establish a non - linear association model between the compressibility factor and temperature, pressure in the hierarchical Bayesian network, and a dynamic relationship model between the specific heat ratio and temperature;

[0053] Among them, for each parameter node of the topological structure, construct a conditional probability distribution based on its parent nodes, and accordingly use the variational inference method to optimize the parameters of the hierarchical Bayesian network to obtain the global joint probability distribution, forming a parameter - associated Bayesian network that characterizes the complex dependence relationship of CAES thermodynamic parameters.

[0054] Specifically, construct the prior knowledge integration and parameter space: Read the complete pre - processed dataset and the adaptive state equation model, analyze the historical data through the Markov Chain Monte Carlo method, extract the distribution characteristics of key thermodynamic parameters (such as specific heat ratio γ, compressibility factor Z, enthalpy value H, etc.), combine the expert experience to set the physical constraint boundaries of parameters, establish a multi - dimensional parameter space, and generate the initial parameter space model.

[0055] Sampling of the prior distribution optimized by particle swarm optimization: Read the initial parameter space model and design an improved particle swarm optimization algorithm for efficient sampling. First, construct the objective function as the error function between the observed data and the model prediction. Initialize the positions and velocity vectors of N particles (typical value N = 200) in the parameter space, where each particle represents a set of parameter combinations. Update the particle positions iteratively, and the position update formula is: X(t + 1) = X(t) + V(t + 1); where X(t) is the position vector of the particle at the t-th iteration; the velocity update uses the improved formula: V(t + 1) = w*V(t) + c1*r1*(Pbest - X(t)) + c2*r2*(Gbest - X(t)) + c3*r3*(X(k) - X(t)); here, the inter-particle mutual learning term c3r3(X(k) - X(t)) is added to enhance the parameter space exploration ability, where X(k) is a randomly selected high-fitness particle, w is the inertia weight (set as a dynamic value between 0.5 and 0.9), c1, c2, and c3 are acceleration coefficients, r1, r2, and r3 are random numbers between 0 and 1, and V(t) is the velocity vector of the particle at the t-th iteration; Pbest is the optimal position experienced by the particle itself; Gbest is the optimal position found by the entire particle swarm so far. Iteratively calculate until convergence or the maximum number of iterations is reached, and finally form an optimized sampling point set.

[0056] Construct a hierarchical Bayesian network and parameter correlation modeling: Use the optimized sampling point set to construct a hierarchical Bayesian network to represent the conditional probability relationship between parameters. First, calculate the mutual dependence degree between parameters through information entropy and mutual information analysis to determine the network topology; then, for each parameter node, construct the conditional probability distribution based on its parent nodes; finally, use the variational inference method to optimize the network parameters to obtain the global joint probability distribution. This network particularly establishes the non-linear correlation between the compression factor Z and temperature T, pressure P, and the dynamic relationship between the specific heat ratio γ and temperature T, forming a parameter correlation Bayesian network.

[0057] Dynamic update of the posterior distribution by sequential Monte Carlo: Read the parameter correlation Bayesian network and design a sequential Monte Carlo algorithm based on importance resampling for real-time parameter estimation. First, initialize M particles (typical value M = 1000) to represent the parameter joint distribution; then, when new observed data is obtained, calculate the likelihood value of each particle and update the particle weights: w_i(t) = w_i(t - 1) * p(y(t)|x_i(t)); where y(t) is the observed data at time t, x_i(t) is the parameter value represented by the i-th particle, and p(y(t)|x_i(t)) is the likelihood function. Then, calculate the effective sample number: Neff = 1 / ∑(w_i 2) When Neff is less than the threshold (a typical value of M / 2), perform the resampling step. Select a new set of particles through polynomial resampling or residual resampling, and reset all particle weights to 1 / M. Finally, update the particle states to generate the posterior distribution of the thermodynamic parameters.

[0058] Parameter Estimation Uncertainty Quantification and Reliability Evaluation: Read the posterior distribution of the thermodynamic parameters for post-processing analysis. Calculate statistics such as the mean, variance, skewness, and kurtosis of each parameter, and construct a 95% confidence interval. Introduce a distribution difference measure based on KL divergence to evaluate the information gain between the prior distribution and the posterior distribution, and quantify the convergence of parameter estimation. At the same time, characterize the uncertainty level of the estimation by calculating the entropy value of the posterior distribution, and evaluate the inference quality in combination with the effective sample number of particles. Finally, form a probability distribution of the thermodynamic parameters including parameter estimation values, uncertainty indicators, and convergence evaluation as the output of this step, and this output will be used in step S23.

[0059] Traditional Bayesian methods have low sampling efficiency in high-dimensional parameter spaces, while simple particle swarm optimization is difficult to quantify uncertainty. This embodiment combines the advantages of both to improve the accuracy and reliability of thermodynamic parameter estimation under non-ideal gas conditions. It overcomes the problem of low sampling efficiency of traditional Bayesian methods in high-dimensional spaces, with the parameter estimation speed increased by 10 - 15 times and better convergence. In particular, the introduction of the parameter interaction network model accurately characterizes the non-linear relationship between the compression factor Z and temperature T, pressure P, and the dynamic relationship between the specific heat ratio γ and temperature T. The accuracy of thermodynamic parameter estimation is increased by 3 times. At the same time, it can adapt to system state changes in real time, shortening the delay time from 30 seconds to 1 - 2 seconds, providing the possibility for fast response control; it also quantifies the uncertainty of parameter estimation, making the decision-making process more reliable and avoiding control mistakes caused by parameter estimation errors.

[0060] According to one aspect of the present application, multi-condition parameter sensitivity analysis and dynamic adjustment are specifically as follows:

[0061] Extract the system operation characteristics from the complete preprocessed dataset, and use unsupervised clustering methods to divide the system operation states into multiple typical conditions to generate a condition characteristic library;

[0062] Based on the condition characteristic library and the probability distribution of thermodynamic parameters, use the sensitivity analysis method of information entropy to calculate the sensitivity indicators of each parameter under different conditions, and generate a parameter importance ranking table;

[0063] Based on the parameter importance ranking table, calculate the interaction sensitivity index between parameters, use graph theory methods to construct a parameter correlation network, identify strongly correlated parameter groups, and generate a parameter interaction network model;

[0064] Using the operating condition feature library and the parameter interaction network model, analyze the dynamic changes in parameter sensitivity during the system operating condition conversion process, adopt the time window sliding technology and the exponential smoothing method to capture the sensitivity evolution trend, and generate a sensitivity evolution model;

[0065] Integrate the parameter importance ranking table, the parameter interaction network model, and the sensitivity evolution model, dynamically adjust the parameter estimation strategy, use a finer grid for high-sensitivity parameters, and adopt a joint estimation strategy for strongly correlated parameter groups to generate a parameter sensitivity matrix;

[0066] Based on the parameter sensitivity matrix and the thermodynamic parameter probability distribution, adjust the model structure and solution method according to the characteristics of different operating conditions, and generate a thermodynamic model adaptable to the characteristics of different operating conditions of CAES.

[0067] Specifically, for the identification and classification of operating condition features: Read the complete preprocessed data set, extract the system operation features (such as pressure ratio, flow rate, temperature range, etc.), and divide the system operation states into multiple typical operating conditions (such as full load, partial load, startup process, rapid response, etc.) through unsupervised clustering methods (such as the improved K-means algorithm). Calculate the center point and boundary for each operating condition to generate an operating condition feature library.

[0068] Quantitative evaluation of parameter importance based on information entropy: Read the operating condition feature library and the thermodynamic parameter probability distribution, and design a sensitivity analysis method based on information entropy. First, for each operating condition k and parameter i, calculate the entropy-based sensitivity index S_i,k through perturbation analysis: S_i,k = ∫(f_k(x|x_i+Δx_i) - f_k(x|x_i)) * log(f_k(x|x_i+Δx_i) / f_k(x|x_i)) dx; where f_k(x|x_i) represents the probability density function of the system output when the parameter i takes the value of x_i under the operating condition k; Δx_i is a small perturbation of the parameter i. This embodiment considers the influence of parameter perturbation on the entire probability distribution rather than only focusing on the mean change. Then, calculate the normalized sensitivity index NS_i,k: NS_i,k = S_i,k / ∑_j S_j,k; and rank the parameter importance under each operating condition according to the size of the sensitivity index to generate a parameter importance ranking table.

[0069] Key parameter identification and interaction effect analysis: Read the parameter importance ranking table, and select the top N (typical value N = 5) parameters for each working condition as key parameters. By introducing the interaction sensitivity index SI_ij,k to evaluate the interaction between parameters: SI_ij,k = S_ij,k - S_i,k - S_j,k; where S_ij,k is the joint sensitivity index when parameters i and j change simultaneously. Calculate the interaction matrix between each pair of parameters, and construct a parameter association network through graph theory methods (such as the minimum spanning tree algorithm) to identify strongly correlated parameter groups and form a parameter interaction network model.

[0070] Analysis of parameter sensitivity evolution during the working condition transition period: Read the working condition feature library and the parameter interaction network model, and analyze the dynamic changes of parameter sensitivity during the transition of the system from one working condition to another. Design a time window sliding technique to continuously calculate the sensitivity index during the working condition transition, and capture the sensitivity evolution trend through the exponential smoothing method: S_i(t) = α * S_i_current + (1 - α) * S_i(t - 1); where α is the smoothing coefficient (typical value 0.2 - 0.3); S_i_current is the sensitivity index value of parameter i calculated within the current time window; S_i(t - 1) is the sensitivity index value of parameter i in the previous time window. Analyze the evolution curve, identify sensitivity mutation points, crossover points and stable regions, and construct a sensitivity evolution model describing the change law of parameter importance.

[0071] Generation of sensitivity-based adaptive parameter estimation strategy: Integrate the parameter importance ranking table, the parameter interaction network model and the sensitivity evolution model, and design a working condition adaptive parameter estimation strategy. For the key parameters under different working conditions, dynamically adjust the calculation resource allocation of parameter estimation, adopt a finer grid or more Monte Carlo samples for highly sensitive parameters; use a joint estimation strategy for strongly correlated parameter groups; focus on the parameters with sensitivity changes during the working condition transition period. At the same time, establish a mapping relationship between parameter importance and estimation accuracy to ensure the reliable estimation of important parameters, and generate a parameter sensitivity matrix containing parameter estimation optimization strategies under different working conditions.

[0072] Adaptive adjustment of multi-working condition thermodynamic model: Read the parameter sensitivity matrix and the probability distribution of thermodynamic parameters, and construct a working condition adaptive thermodynamic model. First, customize the thermodynamic parameter combination for each working condition based on the sensitivity characteristics of parameters under different working conditions; second, design a working condition identification and smooth switching mechanism to ensure the stability of the model during the working condition transition period; finally, adjust the model structure and solution method according to the characteristics of different working conditions (such as using high-order polynomial fitting for working condition 1 and piecewise linear interpolation for working condition 2, etc.) to generate a working condition adaptive thermodynamic model that can adapt to the characteristics of different working conditions.

[0073] Traditional sensitivity analysis is usually based on a single operating condition or the assumption that parameters are independent of each other. In this embodiment, not only the relative importance of parameters under different operating conditions is quantified, but also the complex interaction relationships between parameters and the dynamic evolution law of sensitivity during the process of operating condition conversion are analyzed. Through the sensitivity index based on information entropy, the influence of parameter perturbation on the entire probability distribution is captured, rather than only focusing on the change of the mean value, which more comprehensively characterizes the propagation characteristics of parameter uncertainty. The finally generated thermodynamic model adaptable to operating conditions can dynamically adapt to the change of the system operating state, optimize the allocation of computing resources while ensuring the accurate estimation of key parameters, and improve the accuracy and real-time performance of the state assessment of the compressed air energy storage system. In this embodiment, the key parameters under different operating conditions are accurately identified, and the accuracy of parameter importance ranking is increased from 80% to 95%; the complex dependence relationships between parameters are revealed, and the understanding of the combined influence of parameters on the system is upgraded from qualitative to quantitative; the quantitative description of the dynamic change law of parameter importance is realized for the first time, providing a theoretical basis for parameter optimization during the operating condition switching period; the computing resources are optimally allocated, the estimation accuracy of important parameters is increased by 40%, and the amount of calculation is reduced by 25% at the same time; the adaptability of the model to different operating conditions is improved, and the prediction error is reduced from 5 - 8% to 1.5 - 3%, providing reliable model support for the full-condition optimal operation of the system.

[0074] According to one aspect of the present application, the steps of obtaining a high-precision temperature field model include:

[0075] S31. Initial reconstruction of the temperature field based on sparse measurement points: Using the complete preprocessing data set of finite measurement points in the gas storage reservoir, combined with the geometric model of the gas storage reservoir, an improved radial basis function interpolation method is used to construct the initial temperature field distribution, and the interpolation result is optimized through physical constraints (such as boundary conditions, heat conduction equation) to obtain the initial temperature field distribution model.

[0076] S32. Refinement of the temperature field by physics-guided deep learning: Taking the initial temperature field distribution model as prior knowledge, a deep learning network architecture (heat-neural network) that integrates the heat conduction physical equation is designed, and the loss function is constrained by the spatio-temporal heat conduction partial differential equation to realize the refined reconstruction of the temperature field and form a high-precision temperature field model.

[0077] S33. Identification and quantitative characterization of the temperature stratification phenomenon: Extract the stratification characteristics of the high-precision temperature field model, identify the stratification interface by calculating the temperature gradient in the vertical direction, calculate the stratification stability using the improved Richardson number, and quantify the stratification intensity in combination with the density stratification theory to generate the temperature stratification characteristic index.

[0078] S34. Dynamic Evolution Prediction and Stability Analysis of Temperature Field: Based on a high-precision temperature field model and temperature stratification characteristic indexes, a multi-physics field coupling model including convection, conduction, and radiation is constructed. The evolution of the temperature field is predicted through adaptive time-step integration, and the stability and critical conditions of the temperature field are analyzed to form a temperature field evolution prediction model and stability evaluation indexes.

[0079] This embodiment overcomes the limitation of sparse measuring points in the gas storage reservoir. The reconstruction accuracy is three times higher than that of traditional interpolation methods, and the relative error is less than 2%. The physics-guided deep learning temperature field refinement method integrates the heat conduction physical equation into the neural network, making the reconstruction results strictly satisfy the laws of thermodynamics, significantly improving the physical rationality of the temperature field, and reducing the singularities by more than 95%. The temperature stratification phenomenon identification and quantitative characterization technology has realized the precise positioning of the stratification interface in the gas storage reservoir for the first time and the quantitative evaluation of the stratification intensity, providing a new perspective for understanding the energy flow in the gas storage reservoir. The temperature field dynamic evolution prediction technology extends the prediction time domain of the temperature field from the minute level to the hour level, and controls the prediction error within 3.5%, providing a sufficient decision time window for optimizing the operation strategy.

[0080] As Figure 4 shown, according to one aspect of the present application, the steps of reconstructing the temperature distribution in the gas storage reservoir and obtaining a high-precision temperature field model include:

[0081] Based on the CAES dataset, construct an initial temperature field distribution model;

[0082] Combined with the physical laws of heat conduction, form a physical constraint layer by constructing a physical constraint equation set of the temperature field;

[0083] Taking the initial temperature field distribution model as a reference, construct a heat-neural network architecture including a physical constraint layer;

[0084] Based on the physical constraint equation set of the temperature field, construct a physical constraint loss function including data fitting loss, physical equation loss, and boundary condition loss;

[0085] Obtain the temperature gradient characteristics in the initial temperature field distribution model and generate an adaptive sampling point cloud;

[0086] Using the initial temperature field distribution model as prior knowledge, train the heat-neural network architecture based on the physical constraint loss function and the adaptive sampling point cloud to generate a high-precision temperature field model that satisfies the laws of thermodynamics.

[0087] Specifically, construct the physical constraint equation of the temperature field: Read the initial temperature field distribution model, and based on the basic physical laws of heat conduction, establish a partial differential equation system describing the temperature field distribution in the gas storage reservoir. It includes the three-dimensional heat conduction equation: ρCp(∂T / ∂t) = ∇·(k∇T) + q; where ρ is the density, Cp is the specific heat capacity, k is the thermal conductivity, and q is the heat source term; and the boundary conditions applicable to the gas storage reservoir: -k(∂T / ∂n) = h(T - Tenv); where n is the boundary normal direction, h is the heat transfer coefficient, and Tenv is the ambient temperature. Combine the geometric shape and material properties of the gas storage reservoir to generate the physical constraint equation system of the temperature field.

[0088] Design and construction of the heat-neural network architecture: Read the initial temperature field distribution model and the physical constraint equation system of the temperature field, and design an innovative "heat-neural network" architecture. This network uses a multi-layer perceptron as the basic structure, with the input being the spatial coordinates (x, y, z) and time t, and the output being the temperature T at the corresponding point. The network contains 8 hidden layers, with 128 neurons in each layer, uses the GELU activation function, and adopts an adaptive learning rate optimizer. The innovation of the network lies in embedding a physical constraint layer, which calculates the spatial and temporal derivatives of the temperature field through automatic differentiation, so that the network output conforms to the heat conduction equation. After construction, generate the initial model of the heat-neural network.

[0089] The steps of reading the initial model of the heat-neural network and the physical constraint equation system of the temperature field to construct the physical constraint loss function include: Construct the data fitting loss L_data = (1 / N)∑(T_pred - T_obs) 2 , where T_pred is the temperature predicted by the network, T_obs is the measured temperature, and N is the number of observation points; construct the physical equation loss L_pde = (1 / M)∑(ρCp(∂T / ∂t) - ∇·(k∇T) - q) 2 , where M is the number of sampling points, ∂ is the partial derivative, ρ is the density, Cp is the specific heat capacity, k is the thermal conductivity, q is the heat source term, T is the physical quantity describing the temperature distribution of the system, t is the time variable, and ∇ is the gradient operator; construct the boundary condition loss L_bc = (1 / B)∑(-k(∂T / ∂n) - h(T - Tenv)) 2 , where B is the number of boundary sampling points, h is the heat transfer coefficient, Tenv is the ambient temperature, and n is the boundary normal direction; the physical constraint loss function L_total = α·L_data + β·L_pde + γ·L_bc, where α, β, and γ are weight coefficients, and adopt a dynamic adjustment strategy. At the beginning of training, α is larger, and β and γ are increased as the training progresses to ensure that the network simultaneously satisfies data fitting and physical constraints.

[0090] Point Cloud Sampling Strategy and Adaptive Grid Generation: Read the initial model of the heat-neural network, and design a multi-level sampling strategy to improve computational efficiency and accuracy. First, uniformly sample the entire space to create a background point cloud; then, adaptively encrypt based on the magnitude of the temperature gradient, increasing the sampling density in areas with drastic temperature changes; finally, perform local fine sampling around the observation points and in the boundary regions. The sampling point density is adaptively adjusted by the temperature gradient: density(x, y, z) = base_density·(1 + λ·|▽T| 2 ); where λ is the adjustment coefficient and base_density is the basic distribution of the initial point cloud. This multi-scale adaptive sampling strategy generates an adaptively sampled point cloud.

[0091] Heat-Neural Network Training and Temperature Field Generation: Read the initial model of the heat-neural network, the physical constraint loss function, and the adaptively sampled point cloud, and execute the network training process. Use the mini-batch gradient descent algorithm with a batch size of 1024 and 10,000 training epochs, and use the early stopping strategy to avoid overfitting. During the training process, regenerate the adaptively sampled point cloud every 500 epochs to ensure that difficult areas are fully learned. After training is completed, use the trained network to predict the temperature field on a high-resolution grid to generate a high-precision temperature field model containing millions of grid points.

[0092] Directly encode the partial differential equation of heat conduction into the loss function of the neural network, enabling the network to generate a physically reasonable temperature field even with limited measurement point data. Compared with traditional interpolation-based methods, this embodiment can not only more accurately reconstruct the temperature distribution between measurement points, but also ensure that the results satisfy the laws of thermodynamics, solving the deficiency of pure data-driven methods in terms of physical constraint satisfaction. By designing the "heat-neural network" architecture, this embodiment directly encodes the partial differential equation of heat conduction into the loss function of the neural network, making the network output strictly satisfy the laws of thermodynamics. Compared with traditional interpolation-based methods, the temperature field reconstruction error in sparse measurement point areas is reduced from 8 - 10% to 2 - 3%, and physically unreasonable singular points and oscillation phenomena are completely eliminated. The introduction of the physical constraint loss function increases the network training convergence speed by 5 times and reduces the required training data volume by 70%. The adaptive sampling strategy focuses computing resources on areas with large temperature gradients, increasing the computing efficiency by 8 times while maintaining high precision. This embodiment realizes the high-precision real-time reconstruction of the three-dimensional temperature field of the gas storage reservoir for the first time, providing a solid data foundation for subsequent temperature stratification analysis and exergy assessment, and promoting the transformation of the energy storage system from a "black box" to a "transparent box".

[0093] According to one aspect of the present application, the identification and quantitative characterization of the temperature stratification phenomenon are specifically as follows:

[0094] S331. Vertical Temperature Gradient Calculation and Analysis: Read the high-precision temperature field model and calculate the temperature gradient along the vertical direction of the gas storage reservoir (usually defined as the z-axis): ▽T_z(x, y, z) = ∏T(x, y, z) / ∏z; where T(x, y, z) represents the temperature value of the gas storage reservoir at the spatial coordinate position (x, y, z). Create uniformly distributed vertical profiles (typical value is 100) within the entire gas storage reservoir space. Calculate the vertical temperature gradients at 100 equally spaced points on each profile. Determine the gradient distribution characteristics through statistical analysis, including the mean, variance, maximum value, and distribution pattern. At the same time, calculate the rate of change of temperature with height and generate a vertical temperature gradient distribution map.

[0095] S332. Stratification Interface Identification and Location: Read the vertical temperature gradient distribution map and use a multi-scale edge detection algorithm to identify the regions of abrupt temperature gradient change. First, apply Gaussian smoothing to eliminate noise, and then calculate the second derivative of the gradient: ▽ 2 T_z(x, y, z) = ∏ 2 T(x, y, z) / ∏z 2 ; Locate the position where the temperature gradient changes most violently through the zero-crossing points of the second derivative. For the identified potential interfaces, apply a threshold discrimination: When the average temperature difference between adjacent regions exceeds a preset threshold (typical value is 1.5°C) and the continuous distance exceeds the minimum scale (typical value is 5% of the reservoir height), it is confirmed as a stratification interface. Integrate adjacent interfaces using a spatial clustering method and finally generate a stratification interface position map.

[0096] S333. Improved Richardson Number Calculation and Stratification Stability Evaluation: Read the high-precision temperature field model and the stratification interface position map, calculate the improved Richardson number (Ri) to evaluate the stratification stability: Ri = (g / T0)(∏T / ∏z) / [(∏U / ∏z) 2 + (∏V / ∏z) 2 ; where g is the acceleration due to gravity, T0 is the reference temperature, ∏T / ∏z is the vertical temperature gradient, and ∏U / ∏z and ∏V / ∏z are the vertical gradients of the horizontal velocity components. Considering the difficulty of measuring the gas flow velocity in the gas storage reservoir, a velocity field estimation method based on the pressure field and temperature field is designed: U(x, y, z) = -k1(∏P / ∏x) / μ; V(x, y, z) = -k1(∏P / ∏y) / μ; where k1 is the permeability coefficient and μ is the dynamic viscosity. Judge the stratification stability through the Richardson number: Ri > 0.25 indicates strong stable stratification, 0 < Ri < 0.25 indicates weak stable stratification, and Ri < 0 indicates unstable stratification. Calculate the Richardson number for each identified stratification interface and generate a stratification stability evaluation map.

[0097] S334. Density Stratification Theory and Quantification of Energy Barrier Effect: Read the high-precision temperature field model and the stratification stability evaluation map, and quantify the stratification intensity and energy barrier effect based on the density stratification theory. First, calculate the density at each point according to the temperature: ρ(x, y, z) = P(x, y, z)·M / (R·T(x, y, z)); where P is the pressure, M is the molar mass of the gas, and R is the gas constant. Then calculate the buoyancy frequency (N): N 2 = -(g / ρ0)(∏ρ / ∏z); where g is the acceleration due to gravity and ρ0 is the reference density; the larger the value of N, the stronger the stratification. Introduce the energy transmission coefficient (τ) to quantify the barrier effect of stratification on energy transfer: τ = exp(-∫(N 2 / ω 2 -1) 1 / 2 dz); where ω is the characteristic frequency and the integration interval is the stratification region. The value of τ ranges from 0 to 1, and the closer it is to 0, the stronger the energy barrier effect. Calculate the energy transmission coefficient for each stratification interface to generate the energy barrier effect map.

[0098] S335. Comprehensive Characterization of Stratification Features and Index Construction: Integrate the stratification interface position map, the stratification stability evaluation map, and the energy barrier effect map to construct a comprehensive stratification feature index system. It includes: Stratification Intensity Index (SI): SI = ∑(ΔT_i· h_i · (1 - τ_i)) / H; where ΔT_i is the temperature jump at the i-th stratification interface, h_i is the thickness of this interface, τ_i is the energy transmission coefficient, and H is the total height of the gas storage reservoir; Stratification Complexity Index (CI): CI = -∑p_i·log(p_i); where p_i is the volume proportion of the i-th uniform temperature region; Stratification Stability Index (SSI): SSI = ∑(Ri_i · V_i) / V_total; where Ri_i is the Richardson number of the i-th region, V_i is the volume of this region, and V_total is the total volume.

[0099] Integrate these indicators to generate the temperature stratification feature index, which serves as a quantitative basis for evaluating the impact of non-uniform temperature distribution in the gas storage reservoir. By combining the stratification theory in fluid mechanics with the specific environment of the gas storage reservoir, not only the positions of the stratification interfaces are identified, but also the stratification intensity and stability are quantitatively characterized through the improved Richardson number and energy transmission coefficient, and a comprehensive index system for evaluating the impact of temperature stratification in the gas storage reservoir on system performance is established for the first time.

[0100] In this embodiment, the positioning accuracy of the temperature stratification interface reaches within 2% of the height of the gas storage reservoir, which is 5 times higher than that of the traditional method. It overcomes the problem of difficult measurement of the gas flow velocity in the gas storage reservoir and realizes the quantitative evaluation of the stratification stability in the gas storage reservoir for the first time. By introducing the energy transmission coefficient, the blocking effect of temperature stratification on energy transfer is accurately quantified, providing a new perspective for understanding the energy flow mechanism in the gas storage reservoir. A comprehensive index system including the stratification intensity index, the stratification complexity index, and the stratification stability index is established, which upgrades the temperature stratification evaluation from qualitative description to quantitative characterization and provides a scientific basis for optimizing the operation strategy.

[0101] According to one aspect of the present application, the steps of obtaining the bottleneck analysis map include:

[0102] S41. Identification and quantification of the multi-dimensional energy flow of the system: Using the complete preprocessing data set and the working condition adaptable thermodynamic model, identify various types of energy flows (mechanical energy, thermal energy, pressure energy, etc.) in the system, quantify the magnitude of each energy flow through the energy balance equation, establish an energy flow conversion relationship network, and generate a multi-dimensional energy flow matrix.

[0103] S42. Calculation of entropy generation in sub-regions and localization of irreversible losses: Based on the multi-dimensional energy flow matrix and the high-precision temperature field model, divide the system into multiple control volumes, calculate the entropy generation rate in each control volume using the local entropy balance equation, identify the main positions and types of irreversible losses (such as flow friction, heat transfer, chemical reactions, etc.), and form an entropy generation distribution map.

[0104] S43. Analysis of the available energy of the gas storage reservoir considering temperature stratification: Combining the temperature stratification characteristic index and the entropy generation distribution map, construct a calculation model for the available energy of the gas storage reservoir considering temperature non-uniformity, evaluate the energy quality of the gas storage reservoir through the integration of the available energy in the stratified region, calculate the available energy loss caused by temperature stratification, and obtain the available energy evaluation index of the gas storage reservoir.

[0105] S44. Comprehensive efficiency evaluation and bottleneck identification of the system: Based on the entropy generation distribution map and the available energy evaluation index of the gas storage reservoir, calculate the thermodynamic efficiency of each link of the system, comprehensively evaluate the system performance through generalized efficiency indexes (such as available energy efficiency, energy efficiency, entropy efficiency, etc.), identify the efficiency bottleneck points, and generate a system efficiency evaluation report and a bottleneck analysis map.

[0106] This embodiment realizes the accurate measurement of different forms of energy flow (mechanical energy, thermal energy, pressure energy, etc.), and the energy balance closure rate is increased from 95% to 99%; it accurately identifies the positions and types of the main irreversible losses in the system, improving the loss source positioning accuracy from the system level to the component level, and providing a clear direction for targeted optimization. The available energy analysis technology of the gas storage reservoir considering temperature stratification can more accurately reflect the energy quality in the gas storage reservoir than the traditional uniform model, and the accuracy of efficiency evaluation is increased by 7-12%. Through the comprehensive evaluation of multiple efficiency indicators, the direction of system optimization is made clear and the optimization potential is quantified, providing accurate data support for system improvement.

[0107] According to one aspect of the present application, the calculation of entropy generation in sub-regions and the positioning of irreversible losses are as follows:

[0108] S421. Division of system control volume and definition of boundaries: Read the multi-dimensional energy flow matrix and the high-precision temperature field model, and divide the compressed air energy storage system into multiple control volumes based on the physical structure of the system and the energy flow conversion characteristics. It mainly includes: the compressor unit (further subdivided into multiple-stage compression units), the cooling system, the gas storage reservoir (subdivided into multiple sub-regions according to the temperature field characteristics), and the expansion unit (subdivided into a preheater and multiple-stage expansion units). For each control volume, clearly define the boundary surfaces and boundary conditions, determine the energy and mass inflow / outflow points, and generate the system control volume model.

[0109] S422. Establishment and solution of the local entropy balance equation: Read the system control volume model, and establish the entropy balance equation for each control volume: dS / dt = ∑(m*_in·s_in) - ∑(m*_out·s_out) + ∑(Q*_j / T_j) + S*_gen; where m* is the mass flow rate, s is the specific entropy, Q* is the heat flow rate, T is the boundary temperature, and S*_gen is the entropy generation rate. Through the basic relations of thermodynamics and the state equation, express the entropy as a function of measurable parameters such as temperature and pressure: s = s(T, P) = s0 + cp·ln(T / T0) - R·ln(P / P0); where s0 is the reference specific entropy; cp is the specific heat capacity at constant pressure; T0 is the reference temperature; R is the gas constant; P0 is the reference pressure. Use the numerical integration method to solve the entropy balance equation of each control volume, calculate the entropy generation rate of each control volume, and generate the local entropy generation rate data.

[0110] S423. Decomposition of the entropy generation mechanism and quantification of contributions: Read the local entropy generation rate data, and decompose the total entropy generation rate into the contributions of different physical mechanisms. It mainly includes: entropy generation due to flow friction: S*_gen,fr = ∫[(μ / T)·(∏u_i / ∏x_j + ∏u_j / ∏x_i) 2 dV; entropy generation due to heat transfer: S*_gen,ht = ∫[(k / T 2)·(▽T) 2 dV; Entropy generation due to chemical reaction (if any): S*_gen,ch = -∫(1 / T)·∑(μ_i·r*_i)dV; Entropy generation due to mixing: S*_gen,mix = -R·∑[m*_i·ln(y_i)]; where μ is the dynamic viscosity, u is the velocity component, k is the thermal conductivity, ▽T is the temperature gradient, μ_i is the chemical potential, r*_i is the reaction rate, and y_i is the mole fraction. Calculate the entropy generation rate of each mechanism for each control volume, analyze its relative contribution, and generate entropy generation mechanism distribution data.

[0111] S424. Process path analysis in temperature-entropy (T-s) space: Read the local entropy generation rate data and the thermodynamics model adaptable to the operating conditions, plot the actual process path of the system on the temperature-entropy (T-s) diagram, and compare it with the reversible process (isentropic process). Calculate the irreversible loss of the process: W_loss = ∫T·dS_gen; where S_gen is the entropy generation rate; Use the path deviation degree to quantify the process irreversibility: Δ_path = ∫[|ds_actual / dt - ds_reversible / dt|]dt / ∫[ds_reversible / dt]dt; where ds_actual is the entropy change rate in the actual process path of the system; ds_reversible is the entropy change rate in the reversible process path; The larger the value of Δ_path, the more irreversible the process. Analyze the irreversibility of the key processes of the system (such as compression, expansion, storage, etc.) and generate process irreversibility evaluation data.

[0112] S425. Visualization of multi-scale entropy generation distribution and hot spot identification: Integrate the local entropy generation rate data, entropy generation mechanism distribution data, and process irreversibility evaluation data to create a visualization diagram of the multi-scale entropy generation distribution. Use the adaptive grid refinement technology to increase the grid density in the region with a high entropy generation rate to ensure the accurate capture of entropy generation hot spots. Design an entropy generation intensity index: γ_s = (S*_gen / V) / (S*_gen / V)_avg; When γ_s exceeds the threshold (typical value is 3), it is marked as an entropy generation hot spot. Analyze the physical causes and impacts of each hot spot and generate an entropy generation hot spot analysis report.

[0113] S426. Construction of Entropy Generation Spatiotemporal Distribution Map: Read the local entropy generation rate data, entropy generation mechanism distribution data, and entropy generation hotspot analysis report to construct a complete entropy generation spatiotemporal distribution map. Use three-dimensional visualization technology to display the spatial distribution of entropy generation caused by different physical mechanisms, and combine time series analysis to show the evolution characteristics of entropy generation during the dynamic process of the system. Calculate the contribution ratio of each region and each mechanism to the total irreversible loss through weighted integration, establish a hierarchical entropy generation analysis framework, and finally generate an entropy generation distribution map as an important basis for system efficiency optimization.

[0114] This embodiment combines traditional entropy generation analysis with modern computational fluid dynamics and numerical heat transfer methods to achieve high-resolution spatiotemporal mapping of irreversible losses in a compressed air energy storage system. By decomposing the total entropy generation into the contributions of different physical mechanisms, the root causes of efficiency losses are deeply revealed, providing precise guidance for system optimization. In particular, the introduction of path deviation to quantify process irreversibility provides a new perspective for evaluating the deviation degree between the actual process and the ideal process. This embodiment increases the accuracy rate of loss mechanism identification from 70% to 95%; quantifies process irreversibility through path deviation, promoting the evaluation of irreversible losses from qualitative to quantitative; accurately locates the hotspots of efficiency losses, providing a clear direction for targeted optimization; and achieves high-resolution spatiotemporal mapping of system losses, providing a solid theoretical basis for system optimization design and operation strategy formulation.

[0115] According to one aspect of the present application, the steps of obtaining a bottleneck analysis map for the available energy analysis of a gas storage reservoir considering temperature stratification include:

[0116] Combine the CAES dataset, a working condition adaptable thermodynamic model, and a high-precision temperature field model to construct an entropy generation distribution map; combine it with the high-precision temperature field model, calculate the stratified thermodynamic state data, construct a stratified model and a uniform model of the gas storage reservoir, and conduct a comparative analysis to obtain model difference data;

[0117] Based on the stratified thermodynamic state data, calculate the physical available energy, chemical available energy, and inter-layer available energy transfer terms for each layer to generate stratified available energy data;

[0118] Use the stratified available energy data and the model difference data to calculate the available energy loss caused by temperature stratification, analyze the relationship between the available energy loss rate and the stratification intensity; simultaneously simulate the dynamic change of available energy during the gas storage process, identify the process stage with the most serious available energy loss, and generate available energy loss data and available energy dynamic data;

[0119] Based on this, construct an available energy evaluation index system for the gas storage reservoir, calculate the thermodynamic efficiency of each link of the system, identify the efficiency bottleneck points, and generate a bottleneck analysis map.

[0120] Specifically, for the calculation of the thermodynamic state parameters of the stratified gas storage reservoir: Read the high-precision temperature field model and the temperature stratification characteristic index, and divide the gas storage reservoir space into multiple horizontal layers according to the temperature stratification interface. For each layer, calculate the average temperature \(T_i\), pressure \(P_i\), density \(\rho_i\) and volume \(V_i\). Based on the thermodynamics model adaptable to the working conditions, calculate the thermodynamic parameters of each layer, including specific enthalpy \(h_i\), specific entropy \(s_i\) and specific internal energy \(u_i\): \(h_i = h(T_i, P_i)\); \(s_i = s(T_i, P_i)\); \(u_i = u(T_i, P_i)\); Combine the mass conservation to calculate the gas mass \(m_i=\rho_i\cdot V_i\) of each layer, and generate the stratified thermodynamic state data.

[0121] Comparative analysis of the traditional uniform model and the stratification model: Read the stratified thermodynamic state data, construct the traditional uniform model, and assume that the temperature and pressure in the gas storage reservoir are evenly distributed. Calculate the average temperature \(T_{avg}\), pressure \(P_{avg}\) and the corresponding thermodynamic parameters under the uniform model. By comparing and analyzing the differences between the two models: \(\Delta T = |T_i - T_{avg}|\); \(\Delta P = |P_i - P_{avg}|\); \(\Delta h = |h_i - h_{avg}|\); \(\Delta s = |s_i - s_{avg}|\); Quantitatively evaluate the influence of temperature stratification on the calculation of the thermodynamic state, and generate a model difference evaluation report.

[0122] Accurate calculation of the available energy under stratification conditions: Read the stratified thermodynamic state data, and based on the available energy theory, calculate the physical available energy \(B_{ph,i}\) and chemical available energy \(B_{ch,i}\) of each layer: \(B_{ph,i}=m_i\cdot[(h_i - h_0)-T_0\cdot(s_i - s_0)]\); \(B_{ch,i}=m_i\cdot\sum(\mu_{j,i}-\mu_{j,0})\cdot x_{j,i}\); where \(h_0\), \(s_0\), \(\mu_{j,0}\) are the specific enthalpy, specific entropy and chemical potential under the environmental reference state, \(x_{j,i}\) is the mole fraction of component \(j\), \(m_i\) is the gas mass, \(h_i\) is the specific enthalpy of layer \(i\), \(T_0\) is the temperature of the environmental reference state, \(s_i\) is the specific entropy of layer \(i\), and \(\mu_{j,i}\) is the chemical potential of component \(j\) in layer \(i\). Introduce the available energy transfer term \(B_{trans,i - j}\) between layers, which represents the available energy capture caused by the temperature stratification hindering the heat and mass transfer in the vertical direction: \(B_{trans,i - j}=\sigma_{i - j}\cdot A_{i - j}\cdot T_0\cdot\ln(T_i / T_j)\cdot(1 - e -Pe_i-j ), where \(\sigma_{i - j}\) is the interfacial heat transfer coefficient, \(A_{i - j}\) is the interlayer contact area, \(Pe_{i - j}\) is the Péclet number (measuring the relative intensity of convection and conduction), \(T_i\) and \(T_j\) are the temperatures of layer \(i\) and layer \(j\) respectively, and \(e\) is the natural constant. Summarize the available energy of each layer and generate the calculation results of the stratified available energy.

[0123] Quantification of exergy loss caused by temperature stratification: Read the calculation results of stratified exergy and the model difference evaluation report, and calculate the exergy loss caused by temperature stratification. First, calculate the total exergy B_unif under the ideal uniform model: B_unif = m_total·[(h_avg - h_0) - T_0·(s_avg - s_0)], where m_total is the total mass of the gas, h_avg is the average specific enthalpy, and s_avg is the average specific entropy; then, calculate the actual total exergy B_strat considering stratification: B_strat = ∑B_ph,i + ∑B_ch,i - ∑B_trans,i-j; the exergy loss is defined as the difference between the two: B_loss = B_unif - B_strat; further analyze the relationship between the exergy loss rate and the stratification intensity: η_loss = B_loss / B_unif = f(SI, CI, SSI); where SI, CI, and SSI are the stratification intensity, complexity, and stability index in the temperature stratification characteristic indicators. Through regression analysis, establish a quantitative relationship model between the loss rate and the stratification characteristics, and generate an exergy loss evaluation report.

[0124] Dynamic evolution analysis of exergy during the gas storage process: Read the calculation results of stratified exergy and the temperature field evolution prediction model, and simulate the dynamic changes of exergy during the gas storage process (charging, storage, discharging). Through the time-stepping method, for each time step: update the temperature field and stratification state, recalculate the thermodynamic parameters of each layer; calculate the total exergy at the current moment, compile the exergy time series data, analyze the exergy change rate at key time periods (such as the initial charging period, long-term storage, start of discharging), identify the process stage with the most serious exergy loss, and generate an exergy dynamic evolution diagram.

[0125] Generation of gas storage reservoir optimization indicators based on exergy analysis: Integrate the exergy loss evaluation report and the exergy dynamic evolution diagram to construct an exergy evaluation index system for gas storage reservoirs oriented to optimization. It includes: average exergy efficiency η_B = B_output / B_input; exergy preservation rate *_B = B_end / B_start (during storage); stratification influence factor λ_SI = B_loss / B_input; formulate differentiated evaluation weights for different application scenarios (short-term peak shaving, medium-term load balancing, long-term seasonal energy storage), and generate a comprehensive score. Finally, output the exergy evaluation indicators of the gas storage reservoir containing detailed exergy analysis results and optimization suggestions.

[0126] In this embodiment, the fluid stratification theory is combined with the exergy analysis depth for the first time. By introducing the interlayer exergy transfer term, the influence of temperature non-uniformity on the energy quality of the system is quantitatively evaluated. Compared with the traditional uniform model, it can capture the energy loss mechanism during the gas storage process more accurately, providing a new theoretical perspective and quantitative tool for improving the energy storage efficiency. By dividing the gas storage space into multiple horizontal layers according to the temperature stratification interface and calculating the thermodynamic parameters of each layer, the limitation of the traditional uniform model ignoring temperature non-uniformity is overcome, and the exergy calculation accuracy is improved by 7-12%. Through the comparative analysis of the traditional uniform model and the stratification model, the influence of temperature stratification on the calculation of thermodynamic state is quantitatively evaluated, providing a theoretical basis for model selection. By introducing the interlayer exergy transfer term, the exergy calculation considering the vertical heat transfer resistance is realized for the first time, which more accurately reflects the energy quality and availability in the gas storage reservoir than the traditional method. A quantitative relationship model between the exergy loss rate and the stratification characteristics is established, which upgrades the influence of stratification on the system efficiency from qualitative understanding to quantitative evaluation, providing a clear direction for mitigating stratification and improving the system efficiency.

[0127] According to one aspect of the present application, the steps of generating the CAES operation state classification result include:

[0128] S51. Analysis of the time-frequency characteristics of the energy flow based on wavelet transform: Perform continuous wavelet transform on the multi-dimensional energy flow matrix, extract the spectral characteristics of the energy flow at different time scales, identify periodic fluctuations, transient changes, and long-term trends, separate high-frequency disturbances, medium-frequency load changes, and low-frequency cumulative effects, and form an energy flow time-frequency characteristic diagram.

[0129] S52. Construction and evaluation of multi-time scale state indicators: Based on the energy flow time-frequency characteristic diagram, construct a multi-level evaluation system including instantaneous indicators (response speed, stability, etc.), short-term indicators (energy conversion efficiency, reliability, etc.), and long-term indicators (cumulative efficiency, equipment health, etc.), determine the indicator weights using the analytic hierarchy process, and generate multi-time scale state evaluation indicators.

[0130] S53. Dynamic classification of the operation state and anomaly detection: Use the multi-time scale state evaluation indicators, adopt an unsupervised learning method to identify the typical operation state patterns of the system, construct a dynamic classification model including three states: normal, sub-healthy, and abnormal, calculate the deviation degree of the current state from the typical pattern through the Mahalanobis distance, and form the CAES operation state classification result and anomaly detection alarm.

[0131] S54. Prediction of the system state trend considering time correlation: Based on the CAES operation state classification result and historical data, construct a long short-term memory network to capture the time-dependent relationship of the state indicators, enhance the perception of key time points through the attention mechanism, predict the evolution trend of the system state in the future for a period of time, and generate a state trend prediction report.

[0132] In this embodiment, the fluctuation characteristics of the energy flow are extended from single-time-domain analysis to joint time-frequency-domain analysis, and the accuracy of feature recognition is increased by 45%; a multi-level evaluation system from instantaneous to long-term is established to make the system state evaluation more comprehensive and systematic; the early warning time for abnormal states is extended from the minute level of traditional methods to the hour level, the false alarm rate is reduced by 60%, and the missed alarm rate is reduced by 75%; the state evolution law is captured by a long short-term memory network, and the time domain of trend prediction is extended from 1-2 hours to 24 hours, providing sufficient response time for preventive maintenance and optimized operation.

[0133] According to one aspect of the present application, the steps of forming an energy flow time-frequency feature map include:

[0134] Based on a complete preprocessed data set and a working condition adaptable thermodynamic model, a multi-dimensional energy flow matrix is constructed, and key energy flow signals are extracted therefrom, preprocessed to generate a preprocessed energy flow signal set;

[0135] Perform continuous wavelet transform on the preprocessed energy flow signal set, extract energy distribution feature data, and decompose the preprocessed energy flow signal into high-frequency components, medium-frequency components, and low-frequency components to generate a multi-scale decomposition signal set;

[0136] Analyze the time-frequency correlation between different preprocessed energy flow signals, identify the coupling characteristics between energy flows, and generate an energy flow coupling characteristic map;

[0137] Integrate the energy distribution feature data, the multi-scale decomposition signal set, and the energy flow coupling characteristic map to construct an energy flow time-frequency feature map including time-frequency characteristics.

[0138] Specifically, extraction and preprocessing of multi-scale energy flow signals: Read the multi-dimensional energy flow matrix, and extract key energy flow signals therefrom, including mechanical energy flow E*_mech (compressor input power, expander output power), thermal energy flow E*_thermal (heat flow rate of each heat exchanger), and pressure energy flow E*_press (charging and discharging process of the gas storage reservoir). Perform preprocessing on each signal: Remove outliers and noise (using a median filter); Interpolate missing data (using cubic spline interpolation); Normalize the signal (making the average value of each signal 1); Generate a complete preprocessed energy flow signal set with a time series.

[0139] Optimization and execution of continuous wavelet transform parameters: Read the preprocessed energy flow signal set and design a continuous wavelet transform (CWT) analysis scheme. First, select the most suitable mother wavelet function through signal characteristic analysis (for a compressed air energy storage system, the complex-valued Morlet wavelet is usually the most suitable because of its advantage in time-frequency localization): ψ(t) = π^(-1 / 4)·e^(iω0t)·e^(-t 2 / 2); where ω0 is the center frequency (typical value 5 - 6), ^ represents superscript; e is the natural constant; i is the imaginary unit; t is the time variable. Then, determine the scale range a ∈ [a_min, a_max] and the scale step Δa such that the frequency coverage ranges from high-frequency transients (second level) to low-frequency variations (hour level). Introduce an adaptive scale selection algorithm to dynamically adjust the scale resolution according to the local characteristics of the signal: a_local = a_base·(1 + α·|S(t, a_base)| 2 ); where α is the adjustment parameter, S(t, a_base) is the wavelet coefficient at the base scale, a_local is the adaptively adjusted scale, and a_base is the original scale. Perform a continuous wavelet transform (CWT) on each energy flow signal: W_f(a, b) = (1 / sqrt(a))∫f(t)·ψ*((t - b) / a)dt, where a is the scale parameter, b is the translation parameter, f(t) is the signal, ψ* is the complex conjugate of the mother wavelet, W_f(a, b) is the wavelet transform coefficient matrix, and dt is the time differential. Generate the wavelet transform coefficient matrix.

[0140] Time-frequency energy distribution feature extraction: Read the wavelet transform coefficient matrix and calculate the wavelet energy spectrum: E(a, b) = |W_f(a, b)| 2 ; Extract time-frequency energy distribution features, including: Dominant frequency band identification: For each time point b, find the scale a_max(b) where the energy is most concentrated; Frequency stability assessment: Calculate the rate of change of the dominant frequency band over time σ_freq = std(a_max(b)) / mean(a_max(b)); Energy concentration calculation: C_E(b) = ∑|W_f(a_i, b)| 4 / (∑|W_f(a_i, b)| 2 ) 2 ; The larger the C_E value, the more concentrated the energy is in a specific frequency band; Time-domain fluctuation characteristics: Calculate the wavelet variance Var(a) = (1 / T)∫|W_f(a, b)| 2 db; Generate energy distribution feature data containing various time-frequency features.

[0141] Time-frequency domain decomposition of energy flow signal: Read the wavelet transform coefficient matrix and energy distribution characteristic data, and perform the time-frequency domain decomposition of the signal. Based on the system dynamic characteristics, the energy flow signal is decomposed into three time scales: High-frequency component (second to minute level): f_high(t) = ∑_{a∈A_high}W_f(a, t)·ψ_a,t; corresponding to the fast response, transient process, and perturbation of the system. Medium-frequency component (minute to hour level): f_mid(t) = ∑_{a∈A_mid}W_f(a, t)·ψ_a,t; corresponding to processes such as load change and operating condition conversion. Low-frequency component (hour to day level): f_low(t) = ∑_{a∈A_low}W_f(a, t)·ψ_a,t, corresponding to cumulative effects, long-term trends, etc. Design an adaptive frequency band boundary determination method to automatically identify the system characteristic frequencies through the peaks and valleys of wavelet variance, avoiding the limitations of traditional fixed-frequency band decomposition. Generate a multi-scale decomposition signal set. Where A_high, A_mid, and A_low are the high, medium, and low-frequency scale sets respectively; W_f(a, t) is the wavelet transform coefficient matrix at scale a and time t; ψ_a,t is the wavelet basis function at scale a and time t.

[0142] Time-frequency correlation analysis and energy flow coupling characteristic identification: Read the wavelet transform coefficient matrix and multi-scale decomposition signal set, and analyze the time-frequency correlation between different energy flow signals. Calculate the wavelet cross-correlation function: WCC_fg(a, b) = |W_fg(a, b)| 2 / (|W_f(a, b)| 2 ·|W_g(a, b)| 2 ); where W_fg(a, b) is the wavelet cross-spectrum of signals f and g: W_fg(a, b)= W_f(a, b)·W_g*(a, b); W_g*(a, b) is the conjugate of the wavelet transform coefficient of signal g at scale a and frequency b. Analyze the correlation strength and phase relationship of different energy flows in each frequency band, and identify the coupling characteristics, causal relationships, and transmission delays between energy flows. Introduce wavelet entropy measure to quantify the signal complexity and uncertainty, providing a more comprehensive description of the time-frequency characteristics of energy flow. Generate an energy flow coupling characteristic diagram.

[0143] Energy flow time-frequency feature spectrum generation: Integrate the energy distribution feature data, multi-scale decomposition signal set, and energy flow coupling characteristic map to construct a complete energy flow time-frequency feature spectrum. Use three-dimensional visualization technology, with time and frequency as the coordinate axes and energy intensity as color or height, to intuitively display the time-frequency distribution characteristics of the system's energy flow. Design a feature extraction algorithm to identify typical patterns from the spectrum, such as: energy pulsation pattern: periodic energy fluctuations; energy transfer pattern: energy conversion from one form to another; energy accumulation pattern: energy accumulation on a long time scale; energy dissipation pattern: system loss hotspots. Finally, generate an energy flow time-frequency feature map containing time-frequency characteristics, providing a basis for multi-time scale state assessment.

[0144] This embodiment breaks through the inherent limitation of the traditional Fourier transform in time-frequency resolution and realizes the fine analysis of the energy flow signal in the time-frequency domain. Especially through adaptive scale selection and frequency band boundary determination, the analysis method can be automatically adjusted according to the signal characteristics, better capturing the dynamic behavior of the compressed air energy storage system on different time scales. The introduction of wavelet cross-correlation and wavelet entropy provides new quantification means for the coupling relationship between energy flows and the system complexity. This embodiment overcomes the limitation of insufficient time-frequency resolution of the traditional Fourier transform, realizes the fine analysis of the energy flow signal in the time-frequency domain, and the accuracy of spectrum feature recognition is increased by 45%. It enables the analysis method to be automatically adjusted according to the signal characteristics, improving the ability to capture short-time transients and long-term trends simultaneously. Through main frequency band identification, frequency stability evaluation, and energy concentration calculation, a quantitative description of the system's dynamic characteristics is provided. Decompose the complex energy flow signal into high-frequency components (seconds to minutes), intermediate-frequency components (minutes to hours), and low-frequency components (hours to days), enabling the characteristics of the system on different time scales to be clearly identified. For the first time, the coupling relationship and transmission delay between different energy flows are revealed, providing a theoretical basis for system collaborative optimization.

[0145] According to one aspect of the present application, the steps of generating a CAES operation evaluation and optimization report include:

[0146] S61. Parameter sensitivity analysis for efficiency bottlenecks: Combine the bottleneck analysis spectrum and the parameter sensitivity matrix to identify the key parameters that have the greatest impact on the efficiency of the bottleneck link. Quantify the impact of parameter changes on efficiency through the local response surface method, establish a parameter-efficiency mapping relationship, and form a set of efficiency improvement sensitive parameters.

[0147] S62. Construction of an optimization experience knowledge base based on case-based reasoning: Collect the system's historical operation data and optimization adjustment records, extract the characteristics and optimization strategies of successful cases, calculate the case similarity to match the current state with historical cases, and form an optimization experience knowledge base for different working conditions and problems.

[0148] S63. Solving for optimal parameters using a multi-objective collaborative optimization algorithm: Based on the efficiency improvement sensitive parameter set, the temperature field evolution prediction model, and the system operation constraint conditions, a multi-objective optimization model (objectives such as maximizing efficiency, minimizing entropy generation, and weakening temperature stratification) is constructed, and an improved particle swarm algorithm is used to solve the Pareto optimal solution set to obtain a multi-objective optimization parameter solution.

[0149] S64. Generating operation suggestions and estimating effects: Combining the multi-objective optimization parameter solution with the expert experience in the optimization experience knowledge base, the most suitable optimization solution for the current working conditions is screened through the fuzzy comprehensive evaluation method, the expected effects and possible risks after implementation are evaluated, and a CAES operation optimization suggestion report containing specific parameter adjustment suggestions, expected effects, and precautions is generated.

[0150] In this embodiment, the key parameters affecting the efficiency of each bottleneck link are accurately identified, and the parameter-efficiency mapping accuracy is increased by 3 times; the structuring and digitization of empirical knowledge are realized, enabling the utilization of experience to be upgraded from qualitative to quantitative, and the reliability of the solution is significantly improved; the problem that traditional single-objective optimization ignores the parameter synergy effect is solved, and the overall optimization effect of the system is improved by 15 - 20%; by pre-evaluating the expected effects and implementation risks of the optimization solution, the success rate of solution implementation is increased from 70% to over 90%, reducing the trial-and-error cost of operation optimization and shortening the optimization cycle.

[0151] According to one aspect of the present application, the steps for obtaining the multi-objective optimization parameter solution include:

[0152] Based on the bottleneck analysis map, construct a multi-objective optimization problem model including objective functions for maximizing system efficiency, minimizing entropy generation, weakening temperature stratification, and maximizing response speed;

[0153] Combined with the operation status classification results, based on the current system status and optimization requirements, dynamically allocate the weights of each objective function to generate an adaptive weight solution; use a particle swarm algorithm including a dynamic inertia weight, an adaptive acceleration coefficient, and a chaotic perturbation operator to solve the Pareto optimal solution set;

[0154] Evaluate and analyze the Pareto optimal solution set, extract single-objective optimal solutions, comprehensive index optimal solutions, compromise solutions, and inflection point solutions, generate a Pareto solution set characteristic analysis report and conduct parameter perturbation experiments, evaluate the sensitivity of the objective function to perturbations and implementation risks, and obtain a robustness analysis and risk assessment report;

[0155] Based on the robustness analysis and risk assessment report, select the final optimization solution, and verify and simulate it through a high-precision model to generate a multi-objective optimization parameter solution; form a CAES operation evaluation and optimization report.

[0156] Specifically, formal modeling of the multi-objective optimization problem: Read the sensitive parameter set for improving the reading efficiency, the temperature field evolution prediction model, and the available energy evaluation index of the gas storage reservoir, and construct a multi-objective optimization mathematical model. Define the decision variable vector x, which includes key operating parameters (such as compressor pressure ratio, cooling temperature, gas storage reservoir pressure range, expander inlet temperature, etc.). Set multiple optimization objective functions: Maximize the system efficiency: f1(x) = -η_sys(x); Minimize the entropy generation: f2(x) = S_gen(x); Weaken the temperature stratification: f3(x) = SI(x); Maximize the response speed: f4(x) = -r_resp(x); At the same time, establish a set of constraint conditions: Physical constraints: g_i(x) ≤ 0, i = 1, 2,..., m1; Operating constraints: h_j(x) = 0, j = 1, 2,..., m2; Boundary constraints: x_min ≤ x ≤ x_max; Use the response surface method to construct an approximate model of each objective function to improve the calculation efficiency. Generate a complete multi-objective optimization problem model.

[0157] Adaptive adjustment of the objective function weights: Read the multi-objective optimization problem model and the classification result of the operating state, and design an adaptive adjustment method for the objective function weights. According to the current system state and optimization requirements, dynamically allocate the weights ω_i of each objective function: F(x) = ∑ω_i·f_i(x); Propose a weight determination scheme based on the fuzzy analytic hierarchy process: First, construct a judgment matrix A by expert evaluation, where A_ij represents the importance of objective i relative to objective j; Then calculate the eigenvector as the initial weight; Finally, make dynamic adjustments according to the current system state and operating requirements: ω_i(t) = ω_i,base + Δω_i(S_current, D_current); where S_current is the current system state and D_current is the current demand (such as peak shaving, capacity maximization, life extension, etc.). Generate an adaptive weight scheme.

[0158] Design and Implementation of Improved Particle Swarm Optimization Algorithm: Read the multi-objective optimization problem model and the adaptive weight scheme, and design an improved particle swarm optimization algorithm (MOPSO) to solve the Pareto optimal solution set. First, initialize the particle swarm P(0) (with a typical scale of 100 - 200 particles), where each particle represents a set of decision variable values; set up an external archive E(0) to store non-dominated solutions. During the iteration process: Evaluate the objective function values of each particle; Update the individual best position pbest_i and the global best position gbest; Use the improved velocity and position update formulas: v_i(t + 1) = w·v_i(t) + c1·r1·(pbest_i - x_i(t)) + c2·r2·(gbest - x_i(t)) + c3·r3·(leader_i - x_i(t)); x_i(t + 1) = x_i(t) + v_i(t + 1); Introduce a dynamic inertia weight w and adaptive acceleration coefficients c1, c2, c3, as well as a leader selection strategy based on crowding distance. To solve the problem of being easily trapped in local optima in multi-objective optimization, design a chaotic perturbation operator: x_i(t + 1) = x_i(t + 1) + β·(1 - 2·rand())·x_i(t + 1)·(1 - x_i(t + 1) / x_max); where β is the perturbation intensity (typical value 0.01 - 0.05). Iterate until the termination condition is reached to generate the particle swarm optimization result of the Pareto non-dominated solution set.

[0159] Evaluation of Pareto Solution Set and Extraction of Characteristic Solutions: Read the particle swarm optimization result and evaluate and analyze the obtained Pareto solution set. Calculate the distribution degree index Δ and the spread degree index γ to evaluate the quality of the solution set: Δ = (d_f + d_l + ∑|d_i - d'|) / (d_f + d_l + (N - 1)·d'); γ = sqrt(∑(max_i(f_j(x_i)) - min_i(f_j(x_i))) 2 ); where d_f and d_l are the distances from the boundary points of the solution set to the ideal point, d_i is the distance between adjacent solutions, d' is the average distance, and N is the number of solutions. Extract characteristic solutions from the Pareto solution set, including: each single-objective optimal solution; the optimal solution of the comprehensive index (using the Tchebycheff method); the compromise solution (the solution closest to the ideal point); the inflection point solution (the solution where the sensitivity of the objective function mutates); Generate a Pareto solution set characteristic analysis report.

[0160] Robustness Analysis of Solutions and Implementation Risk Assessment: Read the Pareto solution set feature analysis report and conduct a robustness analysis on the selected feature solutions. Design a parameter perturbation experiment. For each feature solution x*, generate perturbation samples around it: x_pert = x* + Δ·x*·(1 - 2·rand()); where Δ is the perturbation intensity (typical value 0.05 - 0.1). Evaluate the sensitivity of the objective function to the perturbation: S_rob = (1 / N)·∑(|f(x_pert) - f(x*)| / f(x*)); the smaller the S_rob value, the more robust the solution. At the same time, based on historical operation data, evaluate the risks of implementing specific solutions, including factors such as equipment stress, control difficulty, and energy consumption fluctuations. Generate a robustness analysis and risk assessment report.

[0161] Generation and Verification of the Optimal Parameter Scheme: Integrate the Pareto solution set feature analysis report and the robustness analysis and risk assessment report, and select the final optimization scheme based on the multi-attribute decision-making method. Use the fuzzy TOPSIS method, considering multiple attributes such as goal achievement degree, robustness, and risk, calculate the comprehensive distance of each feature solution to the ideal solution, and select the optimal scheme. For the selected scheme, conduct a verification simulation through a high-precision model to ensure that its performance under various working conditions meets expectations. Finally, output a multi-objective optimization parameter scheme including parameter setting values, expected performance improvement, and implementation suggestions.

[0162] In this embodiment, by constructing an optimization model that includes multiple objective functions such as maximizing system efficiency, minimizing entropy generation, weakening temperature stratification, and maximizing response speed, the limitations of traditional single-objective optimization methods are overcome, and the overall optimization effect of the system is improved by 15 - 20%; based on the fuzzy analytic hierarchy process, the optimization process can dynamically adjust the optimization direction according to the system state and requirements, and the adaptability of the optimization scheme is improved; by introducing a dynamic inertia weight, an adaptive acceleration coefficient, and a chaotic perturbation operator, the problem of being easily trapped in local optima in the complex non-linear multi-objective optimization problem of the compressed air energy storage system is effectively solved, and the global optimization ability is improved by 30%; multiple optimization scheme options are provided for decision-makers, making the decision-making process more flexible and scientific. It enables the optimization scheme to move from theory to practice, and the implementation success rate is increased from 70% to over 90%. It effectively solves the complex non-linear multi-objective optimization problem of the compressed air energy storage system. In particular, the adaptive adjustment strategy of the objective function weights enables the optimization process to dynamically adjust the optimization direction according to the system state and requirements, better meeting the actual application needs. It provides a reliable guarantee for the implementation of the scheme, making up for the deficiency of ignoring implementation risks in traditional optimization methods.

[0163] In this embodiment, the combination of particle swarm optimization and Bayesian inference solves the problem that it is difficult to determine the prior distribution in high-dimensional space in traditional Bayesian methods. The particle swarm algorithm intelligently explores the parameter space, provides an optimized initial prior distribution for Bayesian inference, and improves the accuracy and convergence speed of parameter estimation in complex nonlinear systems. Through the "heat-neural network" architecture, the traditional neural network is closely integrated with the heat conduction physical equation. By embedding physical constraints in the loss function, it is ensured that the temperature field generated by the deep learning model under the condition of limited measurement point data satisfies the laws of thermodynamics, solving the dual problems of insufficient accuracy of traditional interpolation methods and physical unreasonableness of pure data-driven methods. The combination of fluid stratification theory and entropy generation analysis quantitatively evaluates the impact of temperature non-uniformity in the gas storage reservoir on the available energy of the system. Through the available energy integration calculation in the stratified region, the energy quality loss ignored by the traditional uniform model is accurately quantified, providing a new theoretical perspective for improving energy storage efficiency. Applying continuous wavelet transform to energy flow analysis realizes the fine analysis of energy flow signals in the time-frequency domain, overcoming the limitation of insufficient time-frequency resolution of traditional Fourier transform. By constructing a hierarchical evaluation index system spanning multiple time scales, the unified evaluation of the compressed air energy storage system at different response levels is realized for the first time.

[0164] This implementation case is applied to a large-scale compressed air energy storage power station with a capacity of 100 MW / 400 MWh. It adopts a cavern-type gas storage reservoir and includes multi-stage compressor units, a heat exchange system, a gas storage cavern, and multi-stage expansion units. The system is equipped with 250 measurement points to collect parameters such as pressure, temperature, flow rate, and power, and the sampling frequency ranges from millisecond level to minute level. The specific steps are as follows:

[0165] Step 1: Multi-source data collection and preprocessing.

[0166] Collect three months of operation data, including: Compressor unit parameters: 6 inter-stage temperatures, pressures, inlet and outlet flow rates, shaft power, etc.; Gas storage reservoir parameters: 20 spatially distributed measurement point temperatures, 5 pressure measurement points, humidity, etc.; Expansion unit parameters: inlet and outlet temperatures, pressures, multi-stage parameters, output power, etc.; Auxiliary system parameters: cooling water temperature, flow rate, heat exchanger temperature, etc.

[0167] Perform data quality assessment. Use the modified Z-score method to calculate the deviation degree of measurement point data: Z_modified =(x_i - median(X)) / (1.4826 * MAD(X)); where MAD is the median absolute deviation. Identify outliers through this method and perform cross-validation in combination with the physical constraints of the system to screen out reliable data.

[0168] Implement multi-scale data synchronization. The data is divided by sampling frequency into: high-frequency data (1 - 10 ms): such as pressure fluctuations; medium-frequency data (1 - 5 s): such as temperature changes; low-frequency data (1 - 5 min): such as energy accumulation. Data synchronization is achieved through timestamp alignment, with the synchronization accuracy reaching within 0.1 seconds. Apply the physical constraint completion method to missing data, construct a partial differential equation system through thermodynamics and fluid mechanics equations for interpolation, and finally obtain a complete preprocessed data set, reducing the data missing rate from 8.5% to 0.9%.

[0169] Step 2: Dynamic estimation of multi-dimensional parameters of the thermodynamic state.

[0170] 2.1 Dynamic selection of the equation of state for non-ideal gases.

[0171] According to the pressure and temperature ranges, construct an adaptive selection mechanism for the equation of state: when P < 10 MPa and 50 K < T < 320 K, select the Redlich-Kwong equation; when P ≥ 10 MPa and 50 K < T < 320 K, select the Peng-Robinson equation; when T ≥ 320 K, dynamically switch between the Benedict-Webb-Rubin equation and the modified BWR equation according to the pressure range. Establish the switching boundary of the equation of state: E_pred = ∑|ρ_pred - ρ_meas| / ρ_meas, and select the equation of state that minimizes E_pred. Here, ρ_pred is the predicted density and ρ_meas is the measured density.

[0172] 2.2 Particle swarm - Bayesian hybrid inference.

[0173] Particle Swarm-Bayesian Hybrid Inference Method: First, construct an initial parameter space model, including the distribution characteristics of the following thermodynamic parameters: Specific heat ratio (γ): 1.3 - 1.4, following a truncated normal distribution; Compression factor (Z): 0.85 - 1.05, related to temperature and pressure; Enthalpy value (H): Deviation based on the reference state, related to temperature and pressure. Use the particle swarm optimization algorithm to generate an optimized sampling point set, set 200 particles, and the position update formula is: X(t + 1) = X(t) + V(t + 1); The velocity update uses an improved formula: V(t + 1) = w*V(t) + c1*r1*(Pbest - X(t)) + c2*r2*(Gbest - X(t)) + c3*r3*(X(k) - X(t)); where w is the inertia weight, with a dynamic value of 0.5 - 0.9; c1, c2, c3 are acceleration coefficients, which are 2.0, 2.0, 1.0 respectively; r1, r2, r3 are random numbers between 0 and 1; Pbest is the individual's historical optimal position; Gbest is the group's historical optimal position; X(k) is the position of a randomly selected high-fitness particle. After performing 200 iterations, select the optimal 50 particles as sampling points.

[0174] Construct a parameter-correlated Bayesian network to calculate the mutual dependence degree between parameters: I(X;Y) = ∑∑p(x, y)log(p(x, y) / (p(x)p(y))); Determine the network topology: The compression factor Z is used as the root node; Temperature T and pressure P are used as the child nodes of Z; The specific heat ratio γ is used as the child node of T. Establish the following non-linear correlation model: Z = f(T, P) = 1 + BP / RT +(CP / RT) 2 , where B and C are parameters to be estimated. The dynamic relationship between the specific heat ratio and temperature: γ = γ0 + aT + bT 2 , where γ0, a, and b are parameters to be estimated. Use the sequential Monte Carlo algorithm to update the posterior distribution of parameters in real time, initialize 1000 particles, and update the particle weights: w_i(t) = w_i(t - 1) * p(y(t)|x_i(t)); where w_i(t) is the weight of the i-th particle at time t; y(t) is the observed data at time t; x_i(t) is the parameter value represented by the i-th particle; p(y(t)|x_i(t)) is the likelihood function. Calculate the effective sample number: Neff = 1 / ∑(w_i 2 ), and when Neff is less than 500, resampling is performed. Quantify the uncertainty of parameter estimation, calculate the 95% confidence interval and the distribution difference measure: KL(P||Q) = ∑P(x)log(P(x) / Q(x)), where P is the posterior distribution and Q is the prior distribution.

[0175] 2.3. Multi-condition parameter sensitivity analysis.

[0176] Sensitivity analysis method based on information entropy: The system operating states are divided into 5 typical working conditions: full-load charging condition, partial-load charging condition, full-load discharging condition, partial-load discharging condition, and standby condition. The improved K-means algorithm is used for clustering. Calculate the sensitivity index: S_i,k = ∫(f_k(x|x_i+Δx_i) - f_k(x|x_i)) *log(f_k(x|x_i+Δx_i) / f_k(x|x_i)) dx; where S_i,k is the sensitivity index of parameter i under working condition k; f_k(x|x_i) is the probability density function of the system output when parameter i takes the value of x_i under working condition k. Calculate the normalized sensitivity index: NS_i,k = S_i,k / ∑_j S_j,k. Analyze the interaction between parameters and calculate the interaction sensitivity index: SI_ij,k =S_ij,k - S_i,k - S_j,k; where S_ij,k is the joint sensitivity index when parameters i and j change simultaneously. During the process of working condition conversion, the time window sliding technology and the exponential smoothing method are used to capture the sensitivity evolution trend: S_i(t) = α * S_i_current + (1-α) * S_i(t-1), where α is the smoothing coefficient and its value is 0.25. The experimental results show that under the full-load charging condition, the sensitivity index of the compression factor Z is the highest (0.42), and the specific heat ratio γ is the second (0.31); while under the partial-load discharging condition, the sensitivity index of the specific heat ratio γ is the highest (0.45).

[0177] Step 3: Three-dimensional reconstruction and evolution prediction of the temperature field in the gas storage reservoir.

[0178] 3.1 Initialization of the temperature field based on sparse measurement points.

[0179] Using the data of 20 spatially distributed temperature measurement points and combining with the geometric model of the gas storage reservoir (a 90m×30m×20m cave), an improved radial basis function interpolation method is used to construct the initial temperature field distribution: T(x, y, z) = ∑λ_i * φ(||x-x_i||); where λ_i is the interpolation coefficient; φ is the radial basis function, and the thin plate spline function φ(r) = r 2 log(r); ||x-x_i|| is the Euclidean distance between the spatial point and the measurement point. Optimize the interpolation result through the boundary conditions and the heat conduction equation.

[0180] 3.2 Refinement of the temperature field by physics-guided deep learning.

[0181] Thermal-neural network architecture: Construct a physical constraint equation set for the temperature field, including the three-dimensional heat conduction equation: ρCp(∏T / ∏t) = ▽·(k▽T) + q; where ρ is the density, about 1.2kg / m for air3 ; Cp is the specific heat capacity, approximately 1005 J / (kg·K); k is the thermal conductivity, approximately 0.026 W / (m·K); q is the heat source term. Boundary condition: -k(∂T / ∂n) = h(T - Tenv); where n is the boundary normal direction; h is the heat transfer coefficient, approximately 5 W / (m 2 ·K); Tenv is the ambient temperature, depending on the rock temperature, approximately 15°C. Design a thermal-neural network architecture with the input being the spatial coordinates (x, y, z) and time t, and the output being the temperature T. The network contains 8 hidden layers with 128 neurons in each layer, and uses the GELU activation function. Construct the physical constraint loss function: L_data = (1 / N)∑(T_pred - T_obs) 2 ; L_pde = (1 / M)∑(ρCp(∂T / ∂t) - ∇·(k∇T) - q) 2 ; L_bc = (1 / B)∑(-k(∂T / ∂n) - h(T - Tenv)) 2 ; L_total = α·L_data + β·L_pde + γ·L_bc. Where N is the number of observation points, 20; M is the number of sampling points, approximately 10,000; B is the number of boundary sampling points, approximately 2,000; α, β, γ are weight coefficients, α = 1.0, β = 0.1, γ = 0.1 at the initial stage of training; α = 0.5, β = 0.3, γ = 0.2 at the later stage of training. Design a multi-level sampling strategy to adaptively adjust the sampling point density based on the temperature gradient: density(x, y, z) = base_density · (1 + λ·|∇T| 2 ), where λ is the adjustment coefficient with a value of 0.2. Use the mini-batch gradient descent algorithm to train the network with a batch size of 1024, 10,000 training epochs, a learning rate of 0.001, and regenerate the adaptive sampling point cloud every 500 epochs. After training, perform temperature field prediction on a high-resolution grid containing 1 million grid points to obtain a high-precision temperature field model with a prediction error of less than 2%.

[0182] 3.3. Identification and quantitative characterization of temperature stratification phenomenon.

[0183] Calculate the temperature gradient along the vertical direction of the gas storage reservoir: ∇T_z(x, y, z) = ∂T(x, y, z) / ∂z. Create 100 uniformly distributed vertical profiles and calculate the vertical temperature gradients at 100 equally spaced points on each profile. Use a multi-scale edge detection algorithm to identify the regions of temperature gradient mutation and calculate the second derivative of the gradient: ∇ 2 T_z(x, y, z) = ∂ 2 T(x, y, z) / ∂z 2。Locate the position with the most drastic change in temperature gradient through the zero-crossing points of the second derivative, and apply threshold discrimination to confirm the stratification interface. Calculate the improved Richardson number (Ri) to evaluate the stratification stability: Ri = (g / T0)(∏T / ∏z) / [(∏U / ∏z) 2 +(∏V / ∏z) 2 ; where g is the acceleration due to gravity, 9.8 m / s 2 ; T0 is the reference temperature, taking the regional average temperature; ∏T / ∏z is the vertical temperature gradient; ∏U / ∏z, ∏V / ∏z are the vertical gradients of the horizontal velocity components. According to the characteristics of the gas storage reservoir, estimate the velocity field based on the pressure field and temperature field: U(x, y, z) = -k1(∏P / ∏x) / μ; V(x, y, z) = -k1(∏P / ∏y) / μ; where k1 is the permeability coefficient, about 10 -12 m 2 ; μ is the dynamic viscosity, about 1.8×10 -5 Pa·s. Calculate the buoyancy frequency: N 2 = -(g / ρ0)(∏ρ / ∏z). Introduce the energy transmission coefficient to quantify the blocking effect of stratification on energy transfer: τ = exp(-∫(N 2 / ω 2 -1) 1 / 2 dz), where ω is the characteristic frequency, taking 0.01 Hz. Construct the stratification characteristic indexes: SI = ∑(ΔT_i · h_i · (1-τ_i)) / H; CI = -∑p_i·log(p_i); SSI = ∑(Ri_i · V_i) / V_total. Where SI is the stratification intensity index; ΔT_i is the temperature jump at the i-th stratification interface; h_i is the thickness of this interface; τ_i is the energy transmission coefficient; H is the total height of the gas storage reservoir; CI is the stratification complexity index; p_i is the volume proportion of the i-th uniform temperature region; SSI is the stratification stability index; Ri_i is the Richardson number of the i-th region; V_i is the volume of this region; V_total is the total volume. The experimental results show that at the end of the charging process, the stratification intensity index reaches a maximum value of 0.42, and the stratification stability index is 0.68, indicating the existence of strong and stable temperature stratification.

[0184] Step 4. System loss analysis and efficiency evaluation based on entropy generation theory.

[0185] 4.1. Division of the system control volume.

[0186] The system is divided into the following control volumes: Compressor unit: divided into low-pressure section, medium-pressure section, and high-pressure section; Cooling system: divided into intercooler and aftercooler; Gas storage: divided into upper, middle, and lower regions according to the temperature field characteristics; Expander unit: divided into preheater, high-pressure section, medium-pressure section, and low-pressure section.

[0187] 4.2 Establishment and solution of the local entropy balance equation.

[0188] Local entropy balance equation and decomposition of entropy generation mechanism. An entropy balance equation is established for each control volume: dS / dt = ∑(m*_in·s_in) - ∑(m*_out·s_out) + ∑(Q*_j / T_j) + S*_gen; where m* is the mass flow rate, kg / s; s is the specific entropy, J / (kg·K); Q* is the heat flow rate, W; T is the boundary temperature, K; S*_gen is the entropy generation rate, W / K. The entropy is expressed as a function of temperature and pressure: s = s(T, P) = s0 + cp·ln(T / T0) - R·ln(P / P0); where s0 is the specific entropy at the reference state, J / (kg·K); cp is the specific heat at constant pressure, J / (kg·K); R is the gas constant, J / (kg·K); T0 and P0 are the temperature and pressure at the reference state. The total entropy generation rate is decomposed into contributions from different physical mechanisms: Entropy generation due to flow friction: S*_gen,fr = ∫[(μ / T)·(∏u_i / ∏x_j + ∏u_j / ∏x_i) 2 dV; Entropy generation due to heat transfer: S*_gen,ht = ∫[(k / T 2 )·(▽T) 2 dV; Entropy generation due to mixing: S*_gen,mix = -R·∑[m*_i·ln(y_i)]; where μ is the dynamic viscosity, Pa·s; u is the velocity component, m / s; k is the thermal conductivity, W / (m·K); ▽T is the temperature gradient, K / m; y_i is the mole fraction. The actual process path of the system is plotted on a temperature-entropy (T-s) diagram, and the irreversible loss during the process is calculated: W_loss = ∫T·dS_gen; The irreversibility of the process is quantified using the path deviation: Δ_path = ∫[|ds_actual / dt - ds_reversible / dt|]dt / ∫[ds_reversible / dt]dt; A visualization graph of the multi-scale entropy generation distribution is created, and an entropy generation intensity index is designed: γ_s = (S*_gen / V) / (S*_gen / V)_avg. The experimental results show that in the high-pressure section of the compressor, the entropy generation due to flow friction accounts for the highest proportion (62%); in the gas storage, the entropy generation due to heat transfer dominates (78%); in the overall system, the compression process contributes 45% of the total entropy generation, the gas storage process 28%, and the expansion process 27%.

[0189] 4.3 Availability analysis of gas storage considering temperature stratification.

[0190] Introduction of the inter-layer availability transfer term. The gas storage is divided into three horizontal layers according to the temperature stratification interface, and the thermodynamic parameters of each layer are calculated: h_i = h(T_i, P_i); s_i = s(T_i, P_i); u_i = u(T_i, P_i); m_i = ρ_i·V_i. A traditional uniform model is constructed to calculate the differences between the parameters and the stratified model: ΔT = |T_i - T_avg|; ΔP = |P_i - P_avg|; Δh = |h_i - h_avg|; Δs = |s_i - s_avg|. Calculate the physical availability and chemical availability of each layer: B_ph,i = m_i·[(h_i - h_0) - T_0·(s_i - s_0)]; B_ch,i = m_i·∑(μ_j,i - μ_j,0)·x_j,i; where B_ph,i is the physical availability of the i-th layer, with the unit of J; B_ch,i is the chemical availability of the i-th layer, with the unit of J; h_0, s_0, μ_j,0 are the specific enthalpy, specific entropy and chemical potential under the environmental reference state; x_j,i is the mole fraction of component j. Introduce the inter-layer availability transfer term: B_trans,i-j = σ_i-j·A_i-j·T_0·ln(T_i / T_j)·(1 - e -Pe_i-j )), where σ_i-j is the interface heat transfer coefficient, W / (m 2 ·K); A_i-j is the inter-layer contact area, m 2; Pe_i-j is the Péclet number, which characterizes the relative intensity of convection and conduction. Calculate the available energy loss caused by temperature stratification: B_unif = m_total·[(h_avg - h_0) - T_0·(s_avg - s_0)]; B_strat = ∑B_ph,i + ∑B_ch,i - ∑B_trans,i-j; B_loss = B_unif - B_strat; η_loss = B_loss / B_unif. Analyze the relationship between the available energy loss rate and the stratification characteristics and find a linear correlation: η_loss = 0.15·SI + 0.08·CI + 0.12·SSI - 0.05. The experimental results show that the available energy loss rate caused by temperature stratification reaches 5.8%, and can be as high as 7.3% during long-term storage. Construct the available energy evaluation indexes for the gas storage reservoir: η_B = B_output / B_input; *_B = B_end / B_start; λ_SI = B_loss / B_input. Among them, η_B is the average available energy efficiency; *_B is the available energy preservation rate; λ_SI is the stratification influence factor.

[0191] Step Five: Multi-time-scale energy flow dynamic tracking and state evaluation.

[0192] 5.1. Analysis of the time-frequency characteristics of the energy flow based on wavelet transform.

[0193] Multi-scale decomposition of the energy flow based on wavelet transform. Extract the key energy flow signals from the multi-dimensional energy flow matrix: mechanical energy flow: compressor input power, expander output power; thermal energy flow: heat flow rate of each heat exchanger; pressure energy flow: gas charging and discharging process of the gas storage reservoir. Select the complex-valued Morlet wavelet as the mother wavelet function: ψ(t) = π^(-1 / 4)·e^(iω0t)·e^(-t 2 / 2); where ω0 is the central frequency, with a value of 6. Adopt the adaptive scale selection algorithm: a_local = a_base·(1 + α·|S(t, a_base)| 2 ); where a_local is the local scale; a_base is the base scale; α is the adjustment parameter, with a value of 0.3; S(t, a_base) is the wavelet coefficient at the base scale. Perform the continuous wavelet transform: W_f(a, b) = (1 / sqrt(a))∫f(t)·ψ*((t - b) / a)dt; where a is the scale parameter, b is the translation parameter, f(t) is the signal, and ψ* is the complex conjugate of the mother wavelet. Calculate the wavelet energy spectrum and extract the time-frequency characteristics: E(a, b) = |W_f(a, b)| 2, find the scale with the most concentrated energy, and calculate the frequency stability and energy concentration: σ_freq = std(a_max(b)) / mean(a_max(b)); C_E(b) = ∑|W_f(a_i, b)| 4 / (∑|W_f(a_i, b)| 2 ) 2 . Decompose the energy flow signal into three time scales: high-frequency component (seconds to minutes): f_high(t) = ∑_{a∈A_high}W_f(a, t)·ψ_a, t; medium-frequency component (minutes to hours): f_mid(t) = ∑_{a∈A_mid}W_f(a, t)·ψ_a, t; low-frequency component (hours to days): f_low(t) = ∑_{a∈A_low}W_f(a, t)·ψ_a, t. Analyze the time-frequency correlation between different energy flow signals, and calculate the wavelet cross-correlation function: WCC_fg(a, b) = |W_fg(a, b)| 2 / (|W_f(a, b)| 2 ·|W_g(a, b)| 2 ), where W_fg(a, b) is the wavelet cross-spectrum of signals f and g: W_fg(a, b) = W_f(a, b)·W_g*(a, b). The experimental results show that in the high-frequency component, the compressor power and the heat exchanger heat flow rate have a correlation coefficient of 0.82, and the phase lag is about 15 seconds; in the low-frequency component, the gas storage pressure energy and the expander output power have a correlation coefficient of 0.93.

[0194] 5.2 Construction and evaluation of multi-time scale state indicators.

[0195] Construct a multi-level evaluation index system. Instantaneous indicators include: response speed: calculate the 10%-90% rise time from the power signal; stability: calculate the volatility of steady-state operation parameters. Short-term indicators include: energy conversion efficiency: output / input energy ratio; reliability: measure through parameter deviation statistics. Long-term indicators include: cumulative efficiency: average efficiency over a long operation period; equipment health: calculate based on vibration and temperature characteristics.

[0196] 5.3 Dynamic classification and anomaly detection of operating states.

[0197] Using the multi-time scale state evaluation index, adopt an unsupervised learning method to identify the typical operating state patterns of the system, and construct three types of states: normal state: all indicators are within the target range; sub-healthy state: some indicators deviate slightly; abnormal state: key indicators deviate severely. Calculate the deviation degree of the current state from the typical pattern through the Mahalanobis distance: D_M(x) = sqrt((x - μ) T Σ -1(x - μ)); where x is the current state vector; μ is the typical mode mean vector; Σ is the covariance matrix. Set the alarm threshold: D_M > 3.0 is the abnormal state, and 1.5 < D_M ≤ 3.0 is the sub-healthy state.

[0198] Step Six: Generate optimization suggestions for operating parameters based on state assessment.

[0199] 6.1. Parameter sensitivity analysis for efficiency bottlenecks.

[0200] Combining the bottleneck analysis map and the parameter sensitivity matrix, identify the key parameters that have the greatest impact on the efficiency of the bottleneck link. Quantify the impact of parameter changes on efficiency through the local response surface method. For the compressor efficiency bottleneck, the key parameter sensitivity ranking is: pressure ratio distribution (0.41); intermediate cooling temperature (0.35); inlet temperature (0.28). For the gas storage reservoir efficiency bottleneck, the key parameter sensitivity ranking is: charging temperature (0.39); charging rate (0.32); initial temperature in the reservoir (0.24).

[0201] 6.2. Solve for the optimal parameters using a multi-objective collaborative optimization algorithm.

[0202] Construct a multi-objective optimization problem model. The objective functions include: f1(x) = -η_sys(x), maximizing the system efficiency; f2(x) = S_gen(x), minimizing the entropy generation; f3(x) = SI(x), weakening the temperature stratification; f4(x) = -r_resp(x), maximizing the response speed. Use the fuzzy analytic hierarchy process to dynamically allocate the weights of each objective function: ω_i(t) = ω_i,base + Δω_i(S_current, D_current); where ω_i,base is the basic weight; S_current is the current system state; D_current is the current demand. Design an improved particle swarm algorithm to solve the Pareto optimal solution set. The velocity and position update formulas are: v_i(t+1) = w·v_i(t) + c1·r1·(pbest_i - x_i(t)) + c2·r2·(gbest - x_i(t)) + c3·r3·(leader_i - x_i(t)); x_i(t+1) = x_i(t) + v_i(t+1). Introduce a dynamic inertia weight: w = w_max - (w_max - w_min) * (t / t_max); where w_max is the initial inertia weight, 0.9; w_min is the final inertia weight, 0.4; t is the current iteration number; t_max is the maximum iteration number. Design a chaotic perturbation operator to enhance the global search ability: x_i(t+1) = x_i(t+1) + β·(1-2·rand())·x_i(t+1)·(1-x_i(t+1) / x_max); where β is the perturbation intensity, with a value of 0.03. Evaluate the Pareto solution set, and calculate the distribution index Δ and the spread index γ: Δ = (d_f + d_l + ∑|d_i - d'|) / (d_f + d_l + (N-1)·d'); γ = sqrt(∑(max_i(f_j(x_i)) - min_i(f_j(x_i))) 2 ) Extract the characteristic solutions, including the single-objective optimal solution, the comprehensive index optimal solution, the compromise solution, and the inflection point solution. Conduct a robustness analysis on the selected characteristic solutions, and design a parameter perturbation experiment: x_pert = x* + Δ·x*·(1-2·rand()); where x* is the characteristic solution; Δ is the perturbation intensity, with a value of 0.08. Evaluate the sensitivity of the objective function to the perturbation: S_rob = (1 / N)·∑(|f(x_pert) - f(x*)| / f(x*)), and select the final optimization scheme based on the multi-attribute decision-making method.

[0203] After adjusting the operating parameters according to the optimization suggestions of this embodiment, the system performance has been improved: the overall operating efficiency of the system has increased by 4.3 percentage points, from the original 62.5% to 66.8%; the energy storage efficiency of the gas storage has increased by 6.2 percentage points, from 85.4% to 91.6%; the equipment life has been extended by 15%, mainly due to the reduction of thermal stress and mechanical stress by optimizing the operating parameters; the operating cost has been reduced by 8.7%, and the annual electricity cost savings are about 4.2 million yuan; the early warning time for abnormal system states has been extended from the original 10 minutes to 2.5 hours, effectively preventing potential failures.

[0204] This embodiment details the actual application process of an optimized operation method for compressed air energy storage based on multi-parameter collaborative control, with a focus on the calculation process and implementation details. Through techniques such as particle swarm-Bayesian hybrid inference, thermal-neural network architecture, decomposition of local entropy balance equations and entropy generation mechanisms, introduction of inter-layer available energy transfer terms, multi-scale decomposition of energy flow based on wavelet transform, multi-objective optimization and adaptive weight schemes, etc., the system operating efficiency has been significantly improved, the equipment life has been extended, and the operating cost has been reduced, fully achieving the invention objectives.

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

Claims

1. A method for optimizing the operation of compressed air energy storage based on multi-parameter coordinated control, characterized in that: include: Collect CAES multi-source sensor data and process them to generate CAES data sets; Based on the CAES data set, the system thermodynamic parameters are estimated by combining particle swarm-Bayesian hybrid inference to form a thermodynamic model that is adaptable to the working conditions. Reconstruct the temperature distribution in the gas storage reservoir based on the CAES data set through the thermal-neural network architecture to obtain a high-precision temperature field model; Combining CAES data sets, operating condition adaptability thermodynamic models and high-precision temperature field models, bottleneck analysis maps are obtained; Generate CAES operation status classification results based on CAES data set and operating condition adaptability thermodynamic model; Combine the bottleneck analysis graph and CAES operation status classification results to generate a CAES operation evaluation optimization report; The steps to obtain the bottleneck analysis map include: Combining CAES data sets, operating condition adaptability thermodynamic models and high-precision temperature field models, an entropy generation distribution map is constructed. Combining it with the high-precision temperature field model, the stratified thermodynamic state data is calculated and a stratified model and a gas storage uniform model are constructed for comparative analysis to obtain model difference data. Based on the layered thermodynamic state data, the physical available energy, chemical available energy and inter-layer available energy transfer items of each layer are calculated to generate layered available energy data; Utilize the stratified available energy data and model difference data to calculate the available energy loss caused by temperature stratification, and simulate the dynamic change of available energy during gas storage to generate available energy loss data and available energy dynamic data; Based on this, we build a gas storage energy evaluation index system, calculate the thermodynamic efficiency of each link in the system, identify efficiency bottlenecks, and generate bottleneck analysis maps; The steps to generate a CAES operation evaluation optimization report include: Based on the bottleneck analysis graph, a multi-objective optimization problem model is constructed, which includes the objective functions of maximizing system efficiency, minimizing entropy generation, reducing temperature stratification, and maximizing response speed; Combined with the CAES operation status classification results, based on the current system status and optimization requirements, the weights of each objective function are dynamically allocated to generate an adaptive weighting scheme; the particle swarm algorithm is used to solve the Pareto optimal solution set; Evaluate and analyze the Pareto optimal solution set, generate a Pareto solution set feature analysis report, conduct parameter perturbation experiments, evaluate the sensitivity of the objective function to perturbations and implementation risks, and obtain a robustness analysis and risk assessment report; Based on the robustness analysis and risk assessment report, the final optimization scheme is selected, and a verification simulation is performed through a high-precision model to generate a multi-objective optimization parameter scheme; thus forming a CAES operation evaluation optimization report.

2. The method according to claim 1, characterized in that The steps to estimate the thermodynamic parameters of the system include: Taking the CAES dataset as input, the PSO algorithm is used to solve the initial parameter space model including the distribution characteristics of thermodynamic parameters and the physical constraint boundaries of parameters to obtain the optimized sampling point set. Based on the optimized sampling point set, the posterior distribution of thermodynamic parameters is updated by using the parameter association Bayesian network that characterizes the conditional probability relationship between parameters and combining with real-time observation data; Based on the posterior distribution of thermodynamic parameters, the uncertainty of parameter estimation is quantified and the probability distribution of thermodynamic parameters, i.e., the thermodynamic parameters of the system, is generated.

3. The method according to claim 2, characterized in that The steps to construct a parameter-associative Bayesian network include: Using the optimized sampling point set, the interdependence between CAES thermodynamic parameters is calculated through information entropy and mutual information analysis, and the topological structure of the hierarchical Bayesian network is determined; A nonlinear correlation model between the compressibility factor and temperature and pressure, as well as a dynamic relationship model between the specific heat ratio and temperature, was established in a hierarchical Bayesian network. For each parameter node in the topological structure, a conditional probability distribution is constructed based on its parent node, and the variational inference method is used to optimize the parameters of the hierarchical Bayesian network to obtain the global joint probability distribution, forming a parameter-associated Bayesian network that characterizes the complex dependencies of CAES thermodynamic parameters.

4. The method according to claim 1, characterized in that The steps of reconstructing the temperature distribution in the gas storage and obtaining a high-precision temperature field model include: Based on the CAES data set, the initial temperature field distribution model is constructed; Combining the physical laws of heat conduction, a physical constraint layer is formed by constructing a set of physical constraint equations of the temperature field; Using the initial temperature field distribution model as a reference, a thermal-neural network architecture including a physical constraint layer is constructed; Based on the temperature field physical constraint equation group, a physical constraint loss function including data fitting loss, physical equation loss and boundary condition loss is constructed; Obtain the temperature gradient characteristics in the initial temperature field distribution model and generate an adaptive sampling point cloud; The initial temperature field distribution model is used as prior knowledge, and the thermal-neural network architecture is trained based on the physical constraint loss function and the adaptive sampling point cloud to generate a high-precision temperature field model that satisfies the laws of thermodynamics.

5. The method according to claim 4, characterized in that The steps to construct the physical constraint loss function include: Construct data fitting loss L_data = (1 / N)∑(T_pred - T_obs) 2 , where T_pred is the network predicted temperature, T_obs is the measured temperature, and N is the number of observation points; Construct the physical equation loss L_pde = (1 / M)∑(ρCp(∏T / ∏t) - ▽·(k▽T) - q) 2 , where M is the number of sampling points, ∏ is the partial derivative, ρ is the density, Cp is the specific heat capacity, k is the thermal conductivity, q is the heat source term, T is the physical quantity describing the temperature distribution of the system, t is the time variable, and ▽ is the gradient operator; Construct boundary condition loss L_bc = (1 / B)∑(-k(∏T / ∏n) - h(T - Tenv)) 2 , where B is the number of boundary sampling points, h is the heat transfer coefficient, Tenv is the ambient temperature, and n is the boundary normal direction; The physical constraint loss function L_total = α·L_data + β·L_pde + γ·L_bc, where α, β, and γ are weight coefficients.

6. The method according to claim 1, characterized in that The steps for calculating the physical available energy, chemical available energy and inter-layer available energy transfer terms for each layer include: Based on the available energy theory, the physical available energy B_ph,i = m_i·[(h_i - h_0) - T_0·(s_i - s_0)] and chemical available energy B_ch,i = m_i·∑(μ_j,i - μ_j,0)·x_j,i are calculated for each layer, where h_0, s_0, μ_j,0 are the specific enthalpy, specific entropy and chemical potential under the environmental reference state, x_j,i is the mole fraction of component j, m_i is the gas mass, h_i is the specific enthalpy of layer i, T_0 is the temperature of the environmental reference state, s_i is the specific entropy of layer i, and μ_j,i is the chemical potential of component j in layer i; Calculate the inter-layer available energy transfer term B_trans,ij = σ_i-j·A_i-j·T_0·ln(T_i / T_j)·(1 - e -Pe_i-j ), where σ_i-j is the interface heat transfer coefficient, A_i-j is the interlayer contact area, Pe_i-j is the Peclet number, T_i and T_j are the temperatures of layer i and layer j respectively, and e is a natural constant.

7. The method according to claim 1, characterized in that Before generating the CAES operation status classification results, it also includes constructing the energy flow time-frequency characteristic diagram, specifically: Based on the CAES data set and the operating condition adaptability thermodynamic model, a multi-dimensional energy flow matrix is ​​constructed and key energy flow signals are extracted and preprocessed to generate a preprocessed energy flow signal set; Performing continuous wavelet transform on the pre-processed energy flow signal set, extracting energy distribution characteristic data and decomposing the pre-processed energy flow signal into high-frequency components, medium-frequency components and low-frequency components, generating a multi-scale decomposition signal set; Analyze the time-frequency correlation between different pre-processed energy flow signals, identify the coupling characteristics between energy flows, and generate energy flow coupling characteristic diagrams; The energy distribution characteristic data, multi-scale decomposition signal set and energy flow coupling characteristic diagram are integrated to construct an energy flow time-frequency characteristic diagram including time-frequency characteristics.

8. The method according to claim 7, characterized in that The steps to perform a continuous wavelet transform include: Based on the characteristics of CAES signals, the complex-valued Morlet wavelet is selected as the mother wavelet function. The complex-valued Morlet wavelet function is expressed as ψ(t) = π^(-1 / 4)·e^(iω0t)·e^(-t 2 / 2), where ω0 is the center frequency; ^ indicates a superscript; e is a natural constant; i is a complex unit; t is a time variable; determines the scale range and scale step size covering from high-frequency transients to low-frequency changes; Adopt the adaptive scale selection algorithm a_local = a_base·(1 + α·|S(t,a_base)| 2 ), where α is the adjustment parameter, S(t, a_base) is the wavelet coefficient at the base scale, a_local is the scale after adaptive adjustment, and a_base is the original scale; Perform continuous wavelet transform W_f(a, b) = (1 / sqrt(a))∫f(t)·ψ*((tb) / a)dt, where a is the scale parameter, b is the translation parameter, f(t) is the signal, ψ* is the complex conjugate of the mother wavelet, W_f(a, b) is the wavelet transform coefficient matrix, and dt is the time differential.

Images