Compressed air energy storage optimized operation method based on multi-parameter cooperative control

Abstract

The invention discloses a compressed air energy storage optimization operation method based on multi-parameter cooperative control, and the method comprises the steps: multi-source data collection and preprocessing, thermodynamic state multi-dimensional parameter dynamic estimation, three-dimensional reconstruction and evolution prediction of a temperature field in a gas storage, and system loss analysis and efficiency evaluation based on an entropy generation theory. And multi-time scale energy flow dynamic tracking and state evaluation, and operation parameter optimization suggestion generation based on state evaluation. According to the method, accurate estimation of non-ideal gas parameters is realized through particle swarm-Bayesian hybrid inference, a temperature field in a gas storage is reconstructed, interlayer available energy transfer items are introduced to quantify temperature stratification influence, the operation efficiency and reliability of the system are improved, the service life of equipment is prolonged, and the operation cost is reduced.

Description

Technical Field

[0001] The present invention belongs to the field of CAES, and in particular to a compressed air energy storage optimization operation method based on multi-parameter coordinated 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 electricity and release energy when needed, effectively solving the problem of absorbing intermittent energy such as wind energy and solar energy. As the proportion of renewable energy continues to increase, how to improve the operating efficiency of the CAES system, reduce energy loss, and extend equipment life has become a key technical challenge in energy transformation. As the basis for optimization control, the CAES system operation status assessment can monitor system performance in real time, predict potential risks, and guide optimization strategies, which is of great significance to improving the economy and reliability of the system. Especially in scenarios such as large-scale applications and long-term peak regulation, accurate status assessment can significantly improve system conversion efficiency, reduce energy storage costs, and provide more reliable support for grid peak and frequency regulation.

[0003] At present, there are five main technical routes for evaluating the operating status of CAES systems. The first category is an evaluation method based on physical models. By establishing a thermodynamic and fluid mechanics model of the system, the ideal gas state equation or a simple modified state equation is used to calculate the system state parameters. The second category is an evaluation method based on data-driven, which uses historical operating data to train statistical models or machine learning models to predict system performance and state changes. The third category is an empirical evaluation method, which establishes a rule base and judgment criteria based on expert experience and operating laws to evaluate the operating status and failure risks of the system. The fourth category is an evaluation method based on energy analysis, which analyzes the system energy conversion process and efficiency loss through energy balance and the first and second laws of thermodynamics. The fifth category is a hybrid evaluation method, which combines physical models with data-driven technology, uses known physical laws to constrain the model structure, and uses measured data to optimize model parameters to improve evaluation accuracy and generalization capabilities.

[0004] The main problems of the existing technology are as follows: First, in terms of thermodynamic state uncertainty, the existing CAES system analysis usually assumes that it follows the ideal state equation, but in fact compressed air exhibits significant non-ideal behavior under high pressure and variable temperature environments. Especially in the multi-stage compression process, the heating, cooling and phase change of the gas lead to dynamic changes in parameters such as specific heat capacity and compression factor. The existing static model or simplified assumption cannot accurately capture these changes, resulting in deviations in the prediction of system operation efficiency. Secondly, regarding the problem of non-uniform temperature distribution in gas storage, there are obvious temperature gradients and stratification phenomena inside the gas storage of large CAES systems (especially underground cave types), resulting in non-uniformity of gas properties in the spatial dimension. The existing control method regards the gas storage as a uniform body, only considering the inlet and outlet parameters, and ignoring the impact of the internal temperature distribution on the system efficiency. Studies have shown that non-uniform temperature in gas storage can lead to increased thermodynamic efficiency losses, and there is a lack of charging and discharging optimization strategies based on the real-time distribution of the internal temperature field. Summary of the invention

[0005] The purpose of the invention is to provide a compressed air energy storage optimization operation method based on multi-parameter coordinated control, in order to solve at least one technical problem existing in the prior art.

[0006] Technical solution: A method for optimizing the operation of compressed air energy storage based on multi-parameter coordinated control, comprising the following steps: 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.

[0007] Beneficial effect: The present invention constructs a complete technical system for evaluating the operating status of a compressed air energy storage system, which overcomes the problems of inaccurate parameter estimation, insufficient temperature field analysis, unclear efficiency loss analysis, and insufficient coordination of multi-time scale characteristics in traditional evaluation methods, thereby improving the overall operating efficiency of the system. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

[0011] Figure 4 It is a flowchart of the steps of reconstructing the temperature distribution in the gas storage provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0012] In order to enable those skilled in the art to better understand the scheme of the present invention, the technical scheme in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of the present invention.

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

[0014] like Figure 1 As shown, a compressed air energy storage optimization operation method based on multi-parameter coordinated control includes the following steps: S1, collect CAES multi-source sensor data and process it to generate a complete pre-processed data set, namely, the CAES data set; Specifically, CAES multi-source sensor data includes compressor unit parameter data, gas storage reservoir parameter data, expansion unit parameter data and auxiliary system parameter data. The compressor unit parameter data can be the inlet and outlet temperature, pressure, flow rate, shaft power, etc. of the multi-stage compressor; the gas storage reservoir parameter data can be the multi-point spatial distribution temperature, pressure, humidity, etc.; the expansion 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.

[0015] S2. Based on the complete preprocessed data set, the system thermodynamic parameters are estimated by combining particle swarm-Bayesian hybrid inference to form a working condition adaptability thermodynamic model; Specifically, particle swarm optimization is similar to a group search strategy and is good at finding the optimal solution, while Bayesian inference can handle uncertainty and the possible distribution of estimated parameters, and combine them to derive and estimate the thermodynamic parameters of the system (such as temperature, pressure, efficiency, etc.) to establish a thermodynamic model that can adapt to different operating conditions (operating conditions). In other words, this model can be adjusted according to actual operating conditions to maintain accuracy and reliability.

[0016] S3, based on the complete preprocessed data set, reconstruct the temperature distribution in the gas storage through the thermal-neural network architecture to obtain a high-precision temperature field model; Specifically, the thermal-neural network architecture is a model that combines thermodynamic knowledge with neural network capabilities. This model can more accurately understand and predict the distribution of heat; using the model to generate accurate temperature distribution maps of gas storage reservoirs helps to more clearly understand the internal thermal characteristics.

[0017] S4. Combining the complete preprocessed data set, the operating condition adaptability thermodynamic model and the high-precision temperature field model, a bottleneck analysis map is obtained; Specifically, by combining a complete preprocessed data set, a condition-adaptive thermodynamic model, and a high-precision temperature field model, a map is generated that can intuitively display key bottlenecks in the system, such as areas with limited performance or parts that need to be optimized.

[0018] S5. Generate CAES operation status classification results based on the complete preprocessed data set and the operating condition adaptability thermodynamic model; Specifically, by inputting the pre-processed data into the condition-adaptive thermodynamic model, the performance of the CAES system under various operating conditions can be analyzed. Then, these operating conditions are classified according to the set classification criteria (such as "high efficiency", "stable", and "low efficiency"). This can help identify the system characteristics under different conditions and facilitate optimization and improvement.

[0019] S6. Combine the bottleneck analysis graph and CAES operation status classification results to generate a CAES operation optimization recommendation report.

[0020] Specifically, by combining bottleneck analysis and operating status classification, specific optimization suggestions 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 status.

[0021] This embodiment can improve the operating efficiency and reliability of the CAES system from the entire process of data collection, model building to optimization suggestion generation; through high-precision models and intelligent analysis methods, it can accurately locate problems and formulate targeted solutions, thereby reducing energy loss, improving energy storage efficiency, and increasing system stability.

[0022] According to one aspect of the present application, the step of processing and generating a CAES data set includes: S11. Layered and graded multi-source sensor data acquisition: Obtain compressor unit parameter data (including inlet and outlet temperatures, pressures, flow rates, shaft power, etc. of multi-stage compressors), gas storage reservoir parameter data (including multi-point spatially distributed temperature, pressure, humidity, etc.), expansion unit parameter data (including inlet and outlet temperatures, pressures, flow rates, output power, etc.) and auxiliary system parameter data (including cooling water temperature, flow rate, heat exchanger temperature, etc.) to form an original multi-source data set.

[0023] S12. Data quality assessment and outlier processing: Perform integrity check, consistency verification and outlier identification on the original multi-source data set, calculate the deviation of the data at each measuring point using the modified Z-score method, and perform cross-validation in combination with the system's physical constraints (such as conservation of mass and conservation of energy) to generate quality assessment indicators and labeled data sets.

[0024] S13. Multi-scale data synchronization and fusion preprocessing: The labeled data set is divided into high-frequency data (millisecond level, such as pressure fluctuation), medium-frequency data (second level, such as temperature change) and low-frequency data (minute level, such as energy accumulation) according to the sampling frequency. Wavelet transform is used to extract multi-scale features, and data synchronization is achieved through timestamp alignment to generate a multi-scale synchronized data set.

[0025] S14. Data completion and reconstruction based on physical constraints: For missing data in multi-scale synchronous data sets, a set of partial differential equations is established in combination with the physical constraints of the system, and data completion is performed through high-order physical interpolation methods. For data that is difficult to complete through physical models, an improved space-time tensor completion algorithm is used to reconstruct the data and obtain a complete preprocessed data set.

[0026] This embodiment achieves data coverage of the entire compressed air energy storage system under all operating conditions and in the entire system through hierarchical and graded multi-source sensor data collection and preprocessing, and adopts a strategy combining the modified Z-score method with system physical constraint cross-validation to improve data quality and increase the outlier detection rate by 33%. Multi-scale data synchronization and fusion preprocessing technology overcomes the data synchronization difficulties caused by different sampling frequencies, enabling accurate alignment of high-frequency data (millisecond level) and low-frequency data (minute level), and reducing the time synchronization error from the original 3-5 seconds to less than 0.1 seconds. Based on physical constraint-based data completion and reconstruction technology, the data missing rate is reduced from the original 7-9% to less than 1%, providing a high-quality data foundation for subsequent status assessment.

[0027] According to one aspect of the present application, the step of forming a working condition adaptability thermodynamic model includes: S21. Dynamic selection of non-ideal gas state equation: Based on the pressure and temperature range in the complete preprocessed data set, the applicable state equation (including Peng-Robinson equation, Redlich-Kwong-Soave equation, etc.) is dynamically selected, and the state equation switching boundary is constructed by minimizing the prediction error to form an adaptive state equation model.

[0028] S22. Real-time parameter estimation based on particle swarm-Bayesian hybrid inference: Using the complete preprocessed data set and the adaptive state equation model, a Bayesian network is constructed to characterize the probabilistic dependency between parameters. The particle swarm optimization is used to initialize the prior distribution, and the sequential Monte Carlo method is used to update the posterior distribution in real time to obtain the probability distribution of thermodynamic parameters.

[0029] S23. Reliability parameter estimation based on uncertainty propagation: Quantify the uncertainty and analyze the propagation of the probability distribution of thermodynamic parameters. Use the multi-level second-order moment method to track the propagation path of uncertainty in the system, establish the credibility evaluation index of parameter estimation, and form a reliability model of thermodynamic parameters.

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

[0031] This embodiment reduces the gas state prediction error from 5-7% of the traditional method to less than 1.5%. The parameter estimation time is shortened from minutes to seconds of the traditional Bayesian method, the parameter accuracy is improved by 2-3 times, and the energy balance closure is improved to more than 98%. By quantifying the uncertainty of parameter estimation, the prediction interval coverage is increased from 85% to more than 95%. It provides a theoretical basis for identifying key parameters in different operating scenarios, making the allocation of computing resources more reasonable, the estimation accuracy of important parameters higher, and the real-time and accuracy of state assessment improved.

[0032] like Figure 2 As shown, according to one aspect of the present application, the step of estimating the thermodynamic parameters of the system includes: Taking the CAES dataset as input, the particle swarm optimization (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, and the optimized sampling point set is obtained; 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.

[0033] like Figure 3 As shown, according to one aspect of the present application, the steps of constructing a parameter-associated 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 were 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.

[0034] Specifically, construct prior knowledge integration and parameter space: read the complete preprocessed data set and adaptive state equation model, analyze historical data through Markov chain Monte Carlo method, extract the distribution characteristics of key thermodynamic parameters (specific heat ratio γ, compression factor Z, enthalpy value H, etc.), set the physical constraint boundaries of parameters based on expert experience, establish a multidimensional parameter space, and generate an initial parameter space model.

[0035] Prior distribution sampling for 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 position and velocity vectors of N particles (typical value N=200) in the parameter space, and each particle represents a set of parameter combinations. The particle position is updated through iteration, where 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 tth iteration; the velocity update adopts 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)); the particle mutual learning term c3r3(X(k)-X(t)) is added here to enhance the parameter space exploration capability, where X(k) is a randomly selected high fitness particle, w is the inertia weight (set to a dynamic value of 0.5-0.9), c1, c2, c3 are acceleration coefficients, r1, r2, r3 are random numbers from 0 to 1, and V(t) is the velocity vector of the particle at the tth iteration; Pbest is the optimal position experienced by the particle itself; Gbest is the optimal position found by the entire particle swarm so far. The iterative calculation is performed until convergence or the maximum number of iterations is reached, and finally an optimized sampling point set is formed.

[0036] Hierarchical Bayesian network construction and parameter association modeling: Using the optimized sampling point set, a hierarchical Bayesian network is constructed to characterize the conditional probability relationship between parameters. First, the degree of mutual dependence between parameters is calculated through information entropy and mutual information analysis to determine the network topology; then, for each parameter node, a conditional probability distribution is constructed based on its parent node; finally, the variational inference method is used to optimize the network parameters to obtain the global joint probability distribution. The network specifically establishes the nonlinear association between the compression factor Z and the temperature T and pressure P, as well as the dynamic relationship between the specific heat ratio γ and the temperature T, forming a parameter association Bayesian network.

[0037] Dynamic update of the posterior distribution of Sequential Monte Carlo: Read the parameter-associated 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 joint distribution of parameters; then, when new observation data is obtained, calculate the likelihood value of each particle and update the particle weight: w_i(t) = w_i(t-1) * p(y(t)|x_i(t)); where y(t) is the observation 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. Next, calculate the number of effective samples: Neff = 1 / ∑(w_i 2); When Neff is less than the threshold (typical value is M / 2), the resampling step is performed to select a new set of particles through polynomial resampling or residual resampling, and the weights of all particles are reset to 1 / M. Finally, the state of the particles is updated to generate the posterior distribution of the thermodynamic parameters.

[0038] Parameter estimation uncertainty quantification and reliability evaluation: Read the posterior distribution of thermodynamic parameters for post-processing analysis, calculate the mean, variance, skewness, kurtosis and other statistics of each parameter, and construct a 95% credible interval. Introduce a distribution difference measure based on KL divergence to evaluate the information gain of the prior distribution and the posterior distribution, and quantify the convergence of parameter estimation. At the same time, the uncertainty level of the estimate is characterized by calculating the entropy value of the posterior distribution, and the inference quality is evaluated in combination with the effective sample number of particles. Finally, a thermodynamic parameter probability distribution containing parameter estimates, uncertainty indicators and convergence evaluation is formed as the output of this step, which will be used in step S23.

[0039] The traditional Bayesian method has low sampling efficiency in high-dimensional parameter space, and 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 the traditional Bayesian method in high-dimensional space, increases the parameter estimation speed by 10-15 times, and has better convergence. In particular, the introduction of the parameter interaction network model enables the nonlinear correlation between the compression factor Z and the temperature T and pressure P, as well as the dynamic relationship between the specific heat ratio γ and the temperature T to be accurately characterized, and the accuracy of thermodynamic parameter estimation is improved by 3 times. At the same time, it can adapt to changes in system states in real time, and the delay time is shortened from 30 seconds to 1-2 seconds, which provides the possibility for rapid response control; it also quantifies the uncertainty of parameter estimation, making the decision-making process more reliable and avoiding control errors caused by parameter estimation errors.

[0040] According to one aspect of the present application, the sensitivity analysis and dynamic adjustment of multi-operating condition parameters are specifically as follows: Extract system operation characteristics from the complete preprocessed data set, use unsupervised clustering method to divide the system operation status into multiple typical working conditions, and generate a working condition feature library; Based on the operating condition feature library and the probability distribution of thermodynamic parameters, the sensitivity analysis method of information entropy is used to calculate the sensitivity index of each parameter under different operating conditions and generate a parameter importance ranking table; Based on the parameter importance ranking table, the interaction sensitivity index between parameters is calculated, and the parameter association network is constructed using graph theory methods to identify strongly correlated parameter groups and generate a parameter interaction network model; Using the operating condition feature library and parameter interaction network model, the dynamic changes of parameter sensitivity during the system operating condition conversion process are analyzed, and the time window sliding technology and exponential smoothing method are used to capture the sensitivity evolution trend and generate a sensitivity evolution model. Integrate the parameter importance ranking table, parameter interaction network model and sensitivity evolution model, dynamically adjust the parameter estimation strategy, use finer grids for highly sensitive parameters, use joint estimation strategy for strongly correlated parameter groups, and generate parameter sensitivity matrix; Based on the parameter sensitivity matrix and the probability distribution of thermodynamic parameters, the model structure and solution method are adjusted according to the characteristics of different working conditions to generate an adaptable thermodynamic model that can adapt to the characteristics of different CAES working conditions.

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

[0042] Quantitative evaluation of parameter importance based on information entropy: read the operating condition feature library and the probability distribution of thermodynamic parameters, and design a sensitivity analysis method based on information entropy. First, for each operating condition k and parameter i, the entropy-based sensitivity index S_i, k is calculated 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 operating condition k; Δx_i is a small disturbance of parameter i. This embodiment takes into account the impact of parameter disturbance on the entire probability distribution, rather than focusing only on mean changes. Then, the normalized sensitivity index NS_i,k is calculated: NS_i,k = S_i,k / ∑_j S_j,k; and according to the size of the sensitivity index, the importance of parameters under each working condition is ranked to generate a parameter importance ranking table.

[0043] Identification of key parameters and analysis of interaction effects: Read the parameter importance ranking table and select the top N (typical value N=5) parameters for each working condition as key parameters. The interaction between parameters is evaluated by introducing the interaction sensitivity index SI_ij,k: 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 parameter pair, and construct a parameter association network through graph theory methods (such as the minimum spanning tree algorithm), identify strongly correlated parameter groups, and form a parameter interaction network model.

[0044] Analysis of parameter sensitivity evolution during operating condition conversion: Read the operating condition feature library and parameter interaction network model to analyze the dynamic changes of parameter sensitivity during the conversion of the system from one operating condition to another. Design a time window sliding technology to continuously calculate the sensitivity index during the operating condition conversion, and capture the sensitivity evolution trend through exponential smoothing: 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 in 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, intersections and stable areas, and construct a sensitivity evolution model that describes the law of parameter importance changes.

[0045] Sensitivity-based adaptive parameter estimation strategy generation: Integrate the parameter importance ranking table, parameter interaction network model and sensitivity evolution model to design a parameter estimation strategy that is adaptive to the working condition. Dynamically adjust the allocation of computing resources for parameter estimation for key parameters under different working conditions, use finer grids or more Monte Carlo samples for highly sensitive parameters; use a joint estimation strategy for strongly correlated parameter groups; and focus on 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 reliable estimation of important parameters and generate a parameter sensitivity matrix containing parameter estimation optimization strategies under different working conditions.

[0046] Adaptive adjustment of thermodynamic models for multiple working conditions: Read the parameter sensitivity matrix and the probability distribution of thermodynamic parameters to build a thermodynamic model that is adaptable to working conditions. First, based on the sensitivity characteristics of the parameters under different working conditions, customize the thermodynamic parameter combination for each working condition; second, design the working condition identification and smooth switching mechanism to ensure the stability of the model during the working condition conversion 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.), and generate a thermodynamic model that is adaptable to the characteristics of different working conditions.

[0047] Traditional sensitivity analysis is usually based on a single operating condition or assumes that parameters are independent of each other. This embodiment not only quantifies the relative importance of parameters under different operating conditions, but also analyzes the complex interactions between parameters and the dynamic evolution of sensitivity during operating condition conversion. Through the sensitivity index based on information entropy, the impact of parameter disturbances on the entire probability distribution is captured, rather than just focusing on mean changes, and the parameter uncertainty propagation characteristics are more comprehensively characterized. The final operating condition adaptability thermodynamic model can dynamically adapt to changes in the system's operating state, while ensuring accurate estimation of key parameters, optimizing the allocation of computing resources, and improving the accuracy and real-time performance of compressed air energy storage system state assessment. This embodiment accurately identifies key parameters under different operating conditions, and improves the accuracy of parameter importance ranking from 80% to 95%; reveals the complex dependencies between parameters, and improves the system's understanding of the joint impact of parameters from qualitative to quantitative; for the first time, it realizes the quantitative description of the dynamic change law of parameter importance, and provides a theoretical basis for parameter optimization during the operating condition switching period; optimizes the allocation of computing resources, improves the estimation accuracy of important parameters by 40%, and reduces the amount of calculation by 25%; improves the adaptability of the model to different operating conditions, and reduces the prediction error from 5-8% to 1.5-3%, providing a reliable model support for the system's optimized operation under all operating conditions.

[0048] According to one aspect of the present application, the step of obtaining a high-precision temperature field model includes: S31. Temperature field initialization and reconstruction based on sparse measurement points: Using the complete preprocessed data set of limited 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. The interpolation results are optimized through physical constraints (such as boundary conditions and heat conduction equations) to obtain the initial temperature field distribution model.

[0049] S32. Physics-guided deep learning temperature field refinement: Taking the initial temperature field distribution model as prior knowledge, a deep learning network architecture (thermal-neural network) integrating the heat conduction physics equation is designed. The loss function is constrained by the spatiotemporal heat conduction partial differential equation to achieve refined reconstruction of the temperature field and form a high-precision temperature field model.

[0050] S33. Identification and quantitative characterization of temperature stratification phenomenon: Extract stratification features from high-precision temperature field models, identify stratification interfaces by calculating vertical temperature gradients, calculate stratification stability using improved Richardson numbers, quantify stratification strength using density stratification theory, and generate temperature stratification characteristic indicators.

[0051] S34. Prediction of dynamic evolution of temperature field and stability analysis: Based on high-precision temperature field model and temperature stratification characteristic indicators, a multi-physics field coupling model including convection, conduction and radiation is constructed. The evolution of temperature field is predicted through adaptive time step integration, and the stability and critical conditions of temperature field are analyzed to form a temperature field evolution prediction model and stability evaluation indicators.

[0052] This embodiment overcomes the limitation of sparse measurement points in the gas storage reservoir, and the reconstruction accuracy is 3 times higher than that of the traditional interpolation method, with a relative error of less than 2%. The physics-guided deep learning temperature field refinement method integrates the heat conduction physics equation into the neural network, so that the reconstruction results strictly meet the laws of thermodynamics, the physical rationality of the temperature field is significantly improved, and the singular points are reduced by more than 95%. The temperature stratification phenomenon identification and quantitative characterization technology has achieved the precise positioning of the stratification interface in the gas storage reservoir for the first time, and the quantitative evaluation of the stratification intensity has provided a new perspective for understanding the energy flow in the gas storage reservoir. The temperature field dynamic evolution prediction technology extends the temperature field prediction time domain from minutes to hours, and the prediction error is controlled within 3.5%, providing a sufficient decision-making time window for optimizing the operation strategy.

[0053] like Figure 4 As shown, according to one aspect of the present application, 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.

[0054] Specifically, the physical constraint equations of the temperature field are constructed: the initial temperature field distribution model is read, and based on the basic physical laws of heat conduction, a group of partial differential equations describing the temperature field distribution in the gas storage reservoir is established. This includes the three-dimensional heat conduction equation: ρCp(∏T / ∏t) = ▽·(k▽T) + q; where ρ is density, Cp is specific heat capacity, k is 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. The physical constraint equations of the temperature field are generated by combining the geometric shape and material properties of the gas storage reservoir.

[0055] Thermal-neural network architecture design and construction: Read the initial temperature field distribution model and the temperature field physical constraint equations, and design an innovative "thermal-neural network" architecture. The network uses a multi-layer perceptron as the basic structure, with spatial coordinates (x, y, z) and time t as input, and the temperature T of the corresponding point as output. The network contains 8 hidden layers, 128 neurons in each layer, uses the GELU activation function, and adopts an adaptive learning rate optimizer. The innovation of the network is that it embeds a physical constraint layer, calculates the spatial and temporal derivatives of the temperature field through automatic differentiation, and makes the network output conform to the heat conduction equation. After the construction is completed, the initial model of the thermal-neural network is generated.

[0056] Read the thermal-neural network initial model and the temperature field physical constraint equation group, and the steps of constructing the physical constraint loss function include: constructing the 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; 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, ▽ 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. A dynamic adjustment strategy is adopted. At the beginning of training, α is large, and β and γ are increased as the training progresses to ensure that the network satisfies both data fitting and physical constraints.

[0057] Point cloud sampling strategy and adaptive mesh generation: Read the initial model of the thermal-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 perform adaptive encryption based on the temperature gradient and increase the sampling density in areas with drastic temperature changes; finally, perform local fine sampling around the observation point and in the boundary area. 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 adaptive sampling point cloud.

[0058] Thermal-neural network training and temperature field generation: Read the thermal-neural network initial model, physical constraint loss function and adaptive sampling point cloud, and execute the network training process. A small batch gradient descent algorithm is used with a batch size of 1024 and a training round number of 10,000. An early stopping strategy is used to avoid overfitting. During the training process, the adaptive sampling point cloud is regenerated every 500 rounds to ensure that the difficult areas are fully learned. After the training is completed, the trained network is used to predict the temperature field on a high-resolution grid to generate a high-precision temperature field model containing millions of grid points.

[0059] The partial differential equation of heat conduction is directly encoded into the loss function of the neural network, so that the network can generate a physically reasonable temperature field even with limited measurement point data. Compared with the traditional interpolation-based method, this embodiment can not only reconstruct the temperature distribution between the measurement points more accurately, but also ensure that the results meet the laws of thermodynamics, solving the shortcomings of the pure data-driven method in satisfying physical constraints. This embodiment directly encodes the partial differential equation of heat conduction into the loss function of the neural network by designing a "heat-neural network" architecture, so that the network output strictly meets the laws of thermodynamics. Compared with the traditional interpolation-based method, the temperature field reconstruction error in the sparse measurement point area is reduced from 8-10% to 2-3%, and the physically unreasonable singular points and oscillation phenomena are completely eliminated. The introduction of the physical constraint loss function increases the convergence speed of network training by 5 times and reduces the amount of training data required by 70%. The adaptive sampling strategy concentrates computing resources in areas with large temperature gradients, increases computing efficiency by 8 times, and maintains high accuracy. 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 available energy evaluation, and promoting the transformation of the energy storage system from a "black box" to a "transparent box".

[0060] According to one aspect of the present application, the identification and quantitative characterization of temperature stratification phenomenon are specifically as follows: S331. Calculation and analysis of vertical temperature gradient: Read the high-precision temperature field model and calculate the temperature gradient along the vertical direction of the gas storage (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 at the spatial coordinate position (x, y, z). Create uniformly distributed vertical profiles (typically 100) in the entire gas storage space, calculate the vertical temperature gradient of 100 equally spaced points on each profile, and determine the gradient distribution characteristics through statistical analysis, including mean, variance, maximum value and distribution form. At the same time, calculate the rate of change of temperature with height to generate a vertical temperature gradient distribution map.

[0061] S332, layered interface identification and positioning: read the vertical temperature gradient distribution map, and use a multi-scale edge detection algorithm to identify the temperature gradient mutation area. First, apply Gaussian smoothing to eliminate noise, and then calculate the second-order derivative of the gradient: ▽ 2 T_z(x,y,z) = ∏ 2 T(x,y,z) / ∏z 2 ; The location where the temperature gradient changes most dramatically is located by the zero crossing point of the second-order derivative. For the identified potential interface, a threshold judgment is applied: when the average temperature difference between adjacent areas exceeds the preset threshold (typically 1.5°C) and the continuous distance exceeds the minimum scale (typically 5% of the reservoir height), it is confirmed as a stratified interface. Combined with the spatial clustering method, the adjacent interfaces are integrated to finally generate a stratified interface location map.

[0062] S333, Improved Richardson number calculation and stratification stability evaluation: Read the high-precision temperature field model and stratification interface position map, calculate the improved Richardson number (Ri) to evaluate stratification stability: Ri = (g / T 0 )(∏T / ∏z) / [(∏U / ∏z) 2 + (∏V / ∏z) 2 ]; where g is the acceleration due to gravity, T 0 is the reference temperature, ∏T / ∏z is the vertical temperature gradient, ∏U / ∏z and ∏V / ∏z are the vertical gradients of the horizontal velocity components. Considering the difficulty in measuring the airflow velocity in the gas storage, a velocity field estimation method based on the pressure field and temperature field is designed: U(x, y, z) = -k 1 (∏P / ∏x) / μ; V(x, y, z) =-k 1 (∏P / ∏y) / μ; where k 1where \(k\) is the permeability coefficient and \(\mu\) is the dynamic viscosity. The stratification stability is judged by the Richardson number: \(Ri > 0.25\) indicates strong stable stratification, \(0 < Ri < 0.25\) indicates weak stable stratification, and \(Ri < 0\) indicates unstable stratification. The Richardson number is calculated for each identified stratification interface to generate a stratification stability evaluation map.

[0063] 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: \(\rho(x,y,z)=P(x,y,z)\cdot M / (R\cdot 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 / \rho 0 )(\partial\rho / \partial z)\); where \(g\) is the acceleration due to gravity and \(\rho\) 0 is the reference density; the larger the value of \(N\), the stronger the stratification. Introduce the energy transmission coefficient (\(\tau\)) to quantify the energy barrier effect of stratification on energy transfer: \(\tau=\exp(-\int(N 2 / \omega 2 - 1) 1 / 2 dz)\); where \(\omega\) is the characteristic frequency and the integration interval is the stratification region. The value of \(\tau\) 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 an energy barrier effect map.

[0064] S335. Comprehensive Characterization of Stratification Characteristics 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 characteristic index system. It includes: Stratification Intensity Index (SI): \(SI=\sum(\Delta T_i\cdot h_i\cdot(1 - \tau_i)) / H\); where \(\Delta T_i\) is the temperature jump at the \(i\)-th stratification interface, \(h_i\) is the thickness of this interface, \(\tau_i\) is the energy transmission coefficient, and \(H\) is the total height of the gas storage reservoir; Stratification Complexity Index (CI): \(CI =-\sum p_i\cdot\log(p_i)\); where \(p_i\) is the volume proportion of the \(i\)-th uniform temperature region; Stratification Stability Index (SSI): \(SSI=\sum(Ri_i\cdot 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.

[0065] Combining these indicators, the temperature stratification characteristic index is generated as a quantitative basis for evaluating the impact of non-uniform temperature distribution in gas storage. Combining the stratification theory in fluid mechanics with the specific environment of gas storage, not only the location of the stratification interface is identified, but also the stratification intensity and stability are quantitatively characterized through the improved Richardson number and energy permeability coefficient. For the first time, a comprehensive index system for evaluating the impact of gas storage temperature stratification on system performance is established.

[0066] This embodiment achieves a positioning accuracy of the temperature stratification interface within 2% of the gas storage reservoir height, which is 5 times higher than the traditional method. It overcomes the difficulty of measuring the airflow 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 permeability coefficient, the barrier 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 indicator system including stratification intensity index, stratification complexity index and stratification stability index is established, which upgrades the temperature stratification evaluation from qualitative description to quantitative characterization, providing a scientific basis for optimizing operation strategies.

[0067] According to one aspect of the present application, the step of obtaining a bottleneck analysis map comprises: S41. Identification and quantification of multi-dimensional energy flows in the system: Utilize a complete pre-processed data set and an adaptable thermodynamic model to the working conditions to identify various energy flows (mechanical energy, thermal energy, pressure energy, etc.) in the system, quantify the size of each energy flow through the energy balance equation, establish an energy flow conversion relationship network, and generate a multi-dimensional energy flow matrix.

[0068] S42. Regional entropy generation calculation and irreversible loss location: Based on the multi-dimensional energy flow matrix and high-precision temperature field model, the system is divided into multiple control volumes. The local entropy balance equation is used to calculate the entropy generation rate in each control volume, identify the main irreversible loss locations and types (such as flow friction, heat transfer, chemical reaction, etc.), and form an entropy generation distribution map.

[0069] S43. Analysis of available energy of gas storage considering temperature stratification: Combine the temperature stratification characteristic index and the entropy generation distribution map to construct a calculation model for available energy of gas storage considering temperature non-uniformity. The energy quality of gas storage is evaluated by integrating the available energy in the stratified area, and the available energy loss caused by temperature stratification is calculated to obtain the evaluation index of available energy of gas storage.

[0070] S44. Comprehensive system efficiency evaluation and bottleneck identification: Based on the entropy generation distribution diagram and the available energy evaluation index of the gas storage reservoir, the thermodynamic efficiency of each link of the system is calculated, and the system performance is comprehensively evaluated through generalized efficiency indicators (such as available energy efficiency, energy efficiency, entropy efficiency, etc.), the efficiency bottleneck points are identified, and the system efficiency evaluation report and bottleneck analysis map are generated.

[0071] This embodiment achieves accurate measurement of different forms of energy flow (mechanical energy, thermal energy, pressure energy, etc.), and the energy balance closure is increased from 95% to 99%. It accurately identifies the location and type of the main irreversible losses of the system, and improves the accuracy of locating the loss source from the system level to the component level, providing a clear direction for targeted optimization. The available energy analysis technology of gas storage considering temperature stratification accurately reflects the energy quality in the gas storage than the traditional uniform model, and the accuracy of efficiency evaluation is improved by 7-12%. Through comprehensive evaluation of multiple efficiency indicators, the system optimization direction is clarified, the optimization potential is quantified, and accurate data support is provided for system improvement.

[0072] According to one aspect of the present application, the calculation of entropy generation and location of irreversible losses in different regions are specifically as follows: S421. System control volume division and boundary definition: Read the multi-dimensional energy flow matrix and high-precision temperature field model, and divide the compressed air energy storage system into multiple control volumes based on the system physical structure and energy flow conversion characteristics. Mainly including: compressor unit (further subdivided into multi-stage compression unit), cooling system, gas storage (subdivided into multiple sub-areas according to temperature field characteristics), expansion unit (subdivided into preheater and multi-stage expansion unit). For each control volume, clearly define the boundary surface and boundary conditions, determine the energy and mass inflow / outflow points, and generate the system control volume model.

[0073] S422. Establishment and solution of 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 thermodynamic relations and state equations, entropy is expressed as a function of measurable parameters such as temperature and pressure: s = s(T, P) = s 0 +cp·ln(T / T 0 ) - R·ln(P / P 0 );where s 0 is the reference specific entropy; cp is the specific heat capacity at constant pressure; T 0 is the reference temperature; R is the gas constant; P 0 The entropy balance equation of each control volume is solved by numerical integration method, the entropy generation rate of each control volume is calculated, and the local entropy generation rate data are generated.

[0074] S423, Entropy generation mechanism decomposition and contribution quantification: Read the local entropy generation rate data and decompose the total entropy generation rate into the contributions of different physical mechanisms. Mainly including: Flow friction entropy generation: S*_gen, fr = ∫[(μ / T)·(∏u_i / ∏x_j+ ∏u_j / ∏x_i) 2 ]dV; Heat transfer entropy generation: S*_gen, ht = ∫[(k / T 2 )·(▽T) 2 ]dV; Chemical reaction entropy generation (if any): S*_gen,ch = -∫(1 / T)·∑(μ_i·r*_i)dV; Mixing entropy generation: 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. The entropy generation rate of each mechanism is calculated for each control volume, and its relative contribution is analyzed to generate the entropy generation mechanism distribution data.

[0075] S424, Process path analysis in temperature-entropy (Ts) space: Read local entropy generation rate data and operating condition adaptability thermodynamic model, draw the actual process path of the system on the temperature-entropy (Ts) 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 path deviation to quantify 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 Δ_path value, the more irreversible the process. Analyze the irreversibility of key processes of the system (such as compression, expansion, storage, etc.) and generate process irreversibility evaluation data.

[0076] S425. Visualization of multi-scale entropy generation distribution and hotspot identification: Integrate local entropy generation rate data, entropy generation mechanism distribution data, and process irreversibility assessment data to create a multi-scale entropy generation distribution visualization. Adaptive grid refinement technology is used to increase grid density in areas with high entropy generation rates to ensure accurate capture of entropy generation hotspots. Design 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 hotspot. For each hotspot, analyze its physical causes and impacts, and generate an entropy generation hotspot analysis report.

[0077] S426. Construction of spatiotemporal distribution map of entropy generation: read local entropy generation rate data, entropy generation mechanism distribution data and entropy generation hot spot analysis report, and construct a complete spatiotemporal distribution map of entropy generation. Use three-dimensional visualization technology to display the spatial distribution of entropy generation caused by different physical mechanisms, and combine time series analysis to display the evolution characteristics of entropy generation in 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.

[0078] 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 compressed air energy storage systems. By decomposing the total entropy generation into contributions from different physical mechanisms, the root cause of efficiency loss is deeply revealed, providing precise guidance for system optimization. In particular, the introduction of path deviation quantification of process irreversibility provides a new perspective for evaluating the degree of deviation between the actual process and the ideal process. This embodiment increases the accuracy of loss mechanism identification from 70% to 95%; quantifies process irreversibility through path deviation, and improves the assessment of irreversible losses from qualitative to quantitative; accurately locates efficiency loss hotspots, providing a clear direction for targeted optimization; achieves high-resolution spatiotemporal mapping of system losses, providing a solid theoretical foundation for system optimization design and operation strategy formulation.

[0079] According to one aspect of the present application, considering the available energy analysis of the temperature-stratified gas storage, the steps of obtaining 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 analyze the relationship between the available energy loss rate and stratification intensity; Simulate the dynamic changes of available energy during gas storage, identify the process stage with the most serious available energy loss, and generate available energy loss data and available energy dynamic data; Based on this, a gas storage available energy evaluation index system is constructed, the thermodynamic efficiency of each link in the system is calculated, the efficiency bottlenecks are identified, and a bottleneck analysis map is generated.

[0080] Specifically, calculation of thermodynamic state parameters of stratified gas storage: read high-precision temperature field model and temperature stratification characteristic index, and divide the gas storage space into multiple horizontal layers according to the temperature stratification interface. For each layer, calculate the average temperature T_i, pressure P_i, density ρ_i and volume V_i. Based on the operating condition adaptability thermodynamic model, 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); combined with mass conservation, calculate the gas mass of each layer m_i = ρ_i·V_i, and generate stratified thermodynamic state data.

[0081] Comparative analysis between traditional uniform model and stratified model: Read the stratified thermodynamic state data, construct the traditional uniform model, and assume that the temperature and pressure in the gas storage are uniformly distributed. Calculate the average temperature T_avg, pressure P_avg and corresponding thermodynamic parameters under the uniform model. By comparing and analyzing the differences between the two models: ΔT = |T_i - T_avg|; ΔP = |P_i - P_avg|; Δh = |h_i - h_avg|; Δs = |s_i - s_avg|; quantitatively evaluate the impact of temperature stratification on thermodynamic state calculations, and generate a model difference evaluation report.

[0082] Accurate calculation of available energy under stratified conditions: Read the stratified thermodynamic state data, and calculate the physical available energy B_ph,i and chemical available energy B_ch,i of each layer based on the available energy theory: 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 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. The interlayer available energy transfer term B_trans,ij is introduced to represent the available energy capture caused by the temperature stratification hindering the transfer of heat and mass in the vertical direction: 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 (a measure of the relative strength of convection and conduction), T_i and T_j are the temperatures of layer i and layer j respectively, and e is a natural constant. The available energy of each layer is summarized to generate the calculation results of the available energy of each layer.

[0083] Quantification of available energy loss caused by temperature stratification: Read the calculation results of stratified available energy and the model difference evaluation report, and calculate the available energy loss caused by temperature stratification. First, calculate the total available energy 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 available energy B_strat considering stratification: B_strat = ∑B_ph, i + ∑B_ch, i - ∑B_trans, ij; the available energy loss is defined as the difference between the two: B_loss = B_unif - B_strat; further analyze the relationship between the available energy loss rate and the stratification intensity: η_loss = B_loss / B_unif = f(SI, CI, SSI); where SI, CI, SSI are the stratification intensity, complexity and stability index in the temperature stratification characteristic index. Through regression analysis, a quantitative relationship model between loss rate and stratification characteristics is established to generate an available energy loss assessment report.

[0084] Analysis of the dynamic evolution of available energy during gas storage: Read the calculation results of stratified available energy and the temperature field evolution prediction model to simulate the dynamic changes of available energy during the gas storage process (inflation, storage, and deflation). 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 available energy at the current moment, compile the available energy time series data, analyze the available energy change rate during key periods (such as the initial stage of inflation, long-term storage, and the beginning of deflation), identify the process stage with the most serious available energy loss, and generate a dynamic evolution diagram of available energy.

[0085] Generation of gas storage optimization indicators based on available energy analysis: Integrate the available energy loss assessment report and the available energy dynamic evolution diagram to build an optimization-oriented gas storage available energy assessment indicator system. Including: average available energy efficiency η_B = B_output / B_input; available energy preservation rate *_B = B_end / B_start (storage period); stratified impact factor λ_SI = B_loss / B_input; formulate differentiated evaluation weights for different application scenarios (short-term peak load regulation, medium-term load balancing, long-term seasonal energy storage) to generate a comprehensive score. The final output includes a gas storage available energy assessment indicator with detailed available energy analysis results and optimization suggestions.

[0086] This embodiment deeply combines the fluid stratification theory with the available energy analysis for the first time. By introducing the interlayer available energy transfer term, the impact of temperature non-uniformity on the energy quality of the system is quantitatively evaluated. Compared with the traditional uniform model, it can more accurately capture the energy loss mechanism in the gas storage process, providing a new theoretical perspective and quantitative tool for improving 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 that ignores temperature non-uniformity is overcome, and the accuracy of available energy calculation is improved by 7-12%; the traditional uniform model and the stratified model are compared and analyzed, and the impact of temperature stratification on the calculation of thermodynamic state is quantitatively evaluated, providing a theoretical basis for model selection. By introducing the interlayer available energy transfer term, the available energy calculation considering vertical heat transfer obstacles is realized for the first time, which reflects the energy quality and availability in the gas storage more accurately than the traditional method. A quantitative relationship model between the available energy loss rate and the stratification characteristics is established, which improves the impact of stratification on system efficiency from qualitative understanding to quantitative evaluation, and provides a clear direction for slowing down stratification and improving system efficiency.

[0087] According to one aspect of the present application, the step of generating a CAES operating status classification result includes: S51. Analysis of time-frequency characteristics of energy flow based on wavelet transform: Perform continuous wavelet transform on the multidimensional energy flow matrix, extract the spectral characteristics of energy flow on 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.

[0088] S52. Construction and evaluation of multi-time scale status indicators: Based on the time-frequency characteristic diagram of energy flow, a multi-level evaluation system is constructed, including instantaneous indicators (response speed, stability, etc.), short-term indicators (energy conversion efficiency, reliability, etc.) and long-term indicators (cumulative efficiency, equipment health, etc.). The hierarchical analysis method is used to determine the indicator weights and generate multi-time scale status evaluation indicators.

[0089] S53. Dynamic classification of operating status and anomaly detection: Utilize multi-time scale status evaluation indicators and adopt unsupervised learning methods to identify the typical operating status patterns of the system, build a dynamic classification model that includes three types of status: normal, sub-healthy, and abnormal, calculate the deviation between the current status and the typical pattern through Mahalanobis distance, and form the CAES operating status classification results and anomaly detection alarms.

[0090] S54. System status trend prediction considering time correlation: Based on the CAES operation status classification results and historical data, a long short-term memory network is constructed to capture the time dependency of status indicators. The perception of key time points is enhanced through the attention mechanism, and the system status evolution trend in the future is predicted to generate a status trend forecast report.

[0091] This embodiment expands the energy flow fluctuation characteristics from a single time domain analysis to a joint analysis of the time and frequency domains, and the feature recognition accuracy is improved by 45%; a multi-level evaluation system from instantaneous to long-term is established to make the system status assessment more comprehensive and systematic; the early warning time for abnormal conditions is extended from minutes to hours in traditional methods, and 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 trend prediction time domain is extended from 1-2 hours to 24 hours, providing sufficient response time for preventive maintenance and optimized operation.

[0092] According to one aspect of the present application, the step of forming an energy flow time-frequency characteristic diagram includes: Based on the complete preprocessing 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 preprocessing 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.

[0093] Specifically, multi-scale energy flow signal extraction and preprocessing: read the multi-dimensional energy flow matrix and extract key energy flow signals from it, 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 (gas storage filling and deflation process). Preprocess each signal: remove outliers and noise (using median filter); interpolate missing data (using cubic spline interpolation); normalize the signal (make the average value of each signal 1); generate a complete set of preprocessed energy flow signals with a time series.

[0094] Continuous wavelet transform parameter optimization and execution: 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 compressed air energy storage system, complex-valued Morlet wavelet is usually the most suitable because of its advantages in time-frequency positioning): ψ(t) = π^(-1 / 4)·e^(iω 0 t)·e^(-t 2 / 2); where ω 0is 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 to cover frequencies 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 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.

[0095] 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: Main frequency band identification: For each time point b, find the scale a_max(b) where the energy is most concentrated; Frequency stability evaluation: Calculate the rate of change of the main 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.

[0096] Time-frequency domain decomposition of energy flow signal: read the wavelet transform coefficient matrix and energy distribution feature data, and perform time-frequency domain decomposition of the signal. Based on the dynamic characteristics of the system, the energy flow signal is decomposed into three time scales: high-frequency component (seconds to minutes): f_high(t) = ∑_{a∈A_high}W_f(a, t)·ψ_a, t; corresponding to system rapid response, transient process and disturbance. Medium-frequency component (minutes to hours): f_mid(t) = ∑_{a∈A_mid}W_f(a, t)·ψ_a, t; corresponding to load changes, working condition conversion and other processes. Low-frequency component (hours to days): 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 frequency through the wavelet variance peak and valley values, 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.

[0097] Time-frequency correlation analysis and energy flow coupling characteristics 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 coefficients 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 measurement to quantify signal complexity and uncertainty, and provide a more comprehensive description of the time-frequency characteristics of energy flows. Generate energy flow coupling characteristic diagrams.

[0098] Energy flow time-frequency characteristic spectrum generation: Integrate energy distribution characteristic data, multi-scale decomposition signal set and energy flow coupling characteristic diagram to construct a complete energy flow time-frequency characteristic spectrum. Use three-dimensional visualization technology, with time and frequency as coordinate axes and energy intensity as color or height, to intuitively display the time-frequency distribution characteristics of the system energy flow. Design feature extraction algorithms 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, an energy flow time-frequency characteristic diagram containing time-frequency characteristics is generated, providing a basis for multi-time scale state assessment.

[0099] This embodiment breaks through the inherent limitations of traditional Fourier transform in time-frequency resolution and realizes the fine analysis of energy flow signals in the time-frequency domain. In particular, through adaptive scale selection and band boundary determination, the analysis method can be automatically adjusted according to the signal characteristics to better capture the dynamic behavior of the compressed air energy storage system at different time scales. The introduction of wavelet cross-correlation and wavelet entropy provides a new means of quantification for the coupling relationship between energy flows and system complexity. This embodiment overcomes the limitations of insufficient time-frequency resolution of traditional Fourier transform, realizes the fine analysis of energy flow signals in the time-frequency domain, and improves the accuracy of spectral feature recognition by 45%. The analysis method can be automatically adjusted according to the signal characteristics, improving the ability to capture short-term transients and long-term trends at the same time. Through main frequency band identification, frequency stability evaluation and energy concentration calculation, a quantitative description of the dynamic characteristics of the system is provided. The complex energy flow signal is decomposed into high-frequency components (seconds to minutes), medium-frequency components (minutes to hours) and low-frequency components (hours to days), so that the characteristics of the system at different time scales can 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.

[0100] According to one aspect of the present application, the step of generating a CAES operation evaluation optimization report includes: S61. Parameter sensitivity analysis for efficiency bottlenecks: Combine the bottleneck analysis map and 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 sensitive parameters for efficiency improvement.

[0101] S62. Construction of optimization experience knowledge base based on case reasoning: Collect historical system operation data and optimization adjustment records, extract the characteristics and optimization strategies of successful cases, match the current status with historical cases through case similarity calculation, and form an optimization experience knowledge base for different working conditions and problems.

[0102] S63. Multi-objective collaborative optimization algorithm to solve the optimal parameters: Based on the efficiency improvement sensitive parameter set, temperature field evolution prediction model and system operation constraints, a multi-objective optimization model (maximizing efficiency, minimizing entropy generation, reducing temperature stratification, etc.) 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.

[0103] S64. Operation suggestion generation and effect estimation: Combine the multi-objective optimization parameter scheme with the expert experience in the optimization experience knowledge base, select the optimization scheme that best suits the current working conditions through the fuzzy comprehensive evaluation method, evaluate the expected effects and possible risks after implementation, and generate a CAES operation optimization suggestion report containing specific parameter adjustment suggestions, expected effects and precautions.

[0104] This embodiment accurately identifies the key parameters that affect the efficiency of each bottleneck link, and the accuracy of parameter-efficiency mapping is improved by 3 times; it realizes the structuring and digitization of empirical knowledge, so that the utilization of experience is upgraded from qualitative to quantitative, and the reliability of the solution is significantly improved; it solves the problem of traditional single-objective optimization ignoring the synergistic effect of parameters, 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 plan, the success rate of the plan implementation is increased from 70% to more than 90%, reducing the trial and error cost of operation optimization and shortening the optimization cycle.

[0105] According to one aspect of the present application, the step of obtaining a multi-objective optimization parameter solution includes: 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 results of the operation status classification, based on the current system status and optimization requirements, the weights of each objective function are dynamically allocated to generate an adaptive weighting scheme; a particle swarm algorithm including dynamic inertia weights, adaptive acceleration coefficients and chaotic perturbation operators is used to solve the Pareto optimal solution set; Evaluate and analyze the Pareto optimal solution set, extract the single-objective optimal solution, comprehensive indicator optimal solution, compromise solution and inflection point solution, generate a Pareto solution set feature analysis report and conduct parameter perturbation experiments to 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.

[0106] Specifically, the multi-objective optimization problem is formally modeled: the efficiency improvement sensitive parameter set, temperature field evolution prediction model and gas storage available energy evaluation index are read, and a multi-objective optimization mathematical model is constructed. The decision variable vector x is defined, which contains key operating parameters (such as compressor pressure ratio, cooling temperature, gas storage pressure range, expander inlet temperature, etc.). Multiple optimization objective functions are set: System efficiency maximization: f 1 (x) = -η_sys(x); Entropy generation is minimized: f 2 (x) = S_gen(x); Temperature stratification weakens: f 3 (x) = SI(x); maximize response speed: f 4 (x) = -r_resp(x); At the same time, establish a set of constraints: physical constraints: g_i(x) ≤ 0, i=1, 2, ..., m 1 ; Operation constraints: h_j(x) = 0, j=1,2,...,m 2 ; Boundary constraints: x_min ≤ x ≤ x_max; Use response surface method to build approximate models of each objective function to improve computational efficiency. Generate a complete multi-objective optimization problem model.

[0107] Adaptive adjustment of objective function weights: Read the multi-objective optimization problem model and the classification results of the operating status, and design an adaptive adjustment method for the objective function weights. According to the current system status 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 target i relative to target j; then calculate the eigenvector as the initial weight; finally, dynamically adjust according to the current system status and operating requirements: ω_i(t) = ω_i, base + Δω_i(S_current, D_current); where S_current is the current system status, and D_current is the current demand (such as peak load regulation, capacity maximization, life extension, etc.). Generate an adaptive weight scheme.

[0108] Improved particle swarm algorithm design and implementation: Read the multi-objective optimization problem model and adaptive weighting scheme, and design an improved particle swarm algorithm (MOPSO) to solve the Pareto optimal solution set. First, initialize the particle swarm P(0) (typical size is 100-200 particles), 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 value of each particle; update the individual optimal position pbest_i and the global optimal position gbest; use the improved speed and position update formula: v_i(t+1) = w·v_i(t) + c1 ·r 1 (pbest_i - x_i(t)) + c 2 ·r 2 (gbest - x_i(t)) + c 3 ·r 3 (leader_i - x_i(t)); x_i(t+1) = x_i(t) + v_i(t+1); dynamic inertia weight w and adaptive acceleration coefficient c are introduced 1 、c 2 、c 3 , and a leader selection strategy based on crowding distance. To solve the problem of being easily trapped in local optimality in multi-objective optimization, a chaotic perturbation operator is designed: 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, and generate the Pareto non-dominated solution set particle swarm optimization result.

[0109] Pareto solution set evaluation and characteristic solution extraction: Read the particle swarm optimization results, evaluate and analyze the obtained Pareto solution set. Calculate the distribution index Δ and the expansion index γ to evaluate the solution set quality: Δ = (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: optimal solutions for each single objective; optimal solutions for comprehensive indicators (using the Tchebycheff method); compromise solutions (solutions closest to the ideal point); inflection point solutions (solutions with sudden changes in the sensitivity of the objective function); generate Pareto solution set characteristic analysis reports.

[0110] Robustness analysis of solutions and implementation risk assessment: Read the characteristic analysis report of the Pareto solution set and perform robustness analysis on the selected characteristic solutions. Design a parameter perturbation experiment, and for each characteristic 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 perturbations: 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 operating 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.

[0111] Generation and verification of optimal parameter solutions: Integrate the Pareto solution set feature analysis report and robustness analysis and risk assessment report, and select the final optimization solution based on the multi-attribute decision-making method. The fuzzy TOPSIS method is used to consider multiple attributes such as target achievement, robustness and risk, calculate the comprehensive distance from each feature solution to the ideal solution, and select the optimal solution. The selected solution is verified and simulated through a high-precision model to ensure that its performance under various working conditions meets expectations. The final output includes a multi-objective optimization parameter solution with parameter setting values, expected performance improvements and implementation recommendations.

[0112] This embodiment overcomes the limitations of the traditional single-objective optimization method by constructing an optimization model that includes multiple objective functions such as maximizing system efficiency, minimizing entropy generation, reducing temperature stratification, and maximizing response speed, and the overall optimization effect of the system is improved by 15-20%; based on the fuzzy hierarchical analysis method, the optimization process can dynamically adjust the optimization direction according to the system state and demand, and the adaptability of the optimization scheme is improved; by introducing dynamic inertia weights, adaptive acceleration coefficients, and chaotic perturbation operators, the problem of being easily trapped in the local optimum in the complex nonlinear multi-objective optimization problem of the compressed air energy storage system is effectively solved, and the global optimization ability is improved by 30%; it provides decision makers with a variety of optimization schemes to make the decision-making process more flexible and scientific. The optimization scheme is moved from theory to practice, and the implementation success rate is increased from 70% to more than 90%. It effectively solves the complex nonlinear multi-objective optimization problem of the compressed air energy storage system. In particular, the adaptive adjustment strategy of the objective function weight enables the optimization process to dynamically adjust the optimization direction according to the system state and demand, and better meet the actual application needs. It provides a reliable guarantee for the implementation of the scheme and makes up for the deficiency of ignoring the implementation risk in the traditional optimization method.

[0113] In this embodiment, particle swarm optimization is combined with Bayesian inference to solve the problem that the prior distribution of the traditional Bayesian method is difficult to determine in high-dimensional space. 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 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 fluid stratification theory is combined with entropy generation analysis to quantitatively evaluate the impact of temperature non-uniformity of gas storage on the available energy of the system. Through the calculation of the available energy integral in the stratified area, the energy mass loss ignored by the traditional uniform model is accurately quantified, providing a new theoretical perspective for improving energy storage efficiency. The continuous wavelet transform is applied to energy flow analysis to achieve a 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 across multiple time scales, a unified evaluation of compressed air energy storage systems at different response levels was achieved for the first time.

[0114] This implementation case is applied to a large compressed air energy storage power station with a capacity of 100MW / 400MWh. It uses a cavernous gas storage facility, including a multi-stage compressor unit, a heat exchange system, a gas storage cave, and a multi-stage expansion unit. The system has 250 measurement points to collect parameters such as pressure, temperature, flow, and power, with sampling frequencies ranging from milliseconds to minutes. The specific steps include: Step 1: Multi-source data collection and preprocessing.

[0115] Three months of operating data were collected, including: compressor unit parameters: 6 inter-stage temperatures, pressures, inlet and outlet flows, shaft power, etc.; gas storage reservoir parameters: 20 spatially distributed measuring point temperatures, 5 pressure measuring points, humidity, etc.; expansion unit parameters: inlet and outlet temperatures, pressures, multiple inter-stage parameters, output power, etc.; auxiliary system parameters: cooling water temperature, flow, heat exchanger temperature, etc.

[0116] Perform data quality assessment and use the modified Z-score method to calculate the deviation of the measurement point data: Z_modified =(x_i - median(X)) / (1.4826 * MAD(X)); where MAD is the median absolute deviation. This method is used to identify outliers and perform cross-validation in combination with the system's physical constraints to screen out reliable data.

[0117] Implement multi-scale data synchronization. Divide the data according to the 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. Achieve data synchronization through timestamp alignment, with the synchronization accuracy reaching within 0.1 seconds. Apply a 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%.

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

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

[0120] 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. Where ρ_pred is the predicted density and ρ_meas is the measured density.

[0121] 2.2 Particle swarm-Bayesian hybrid inference.

[0122] Particle swarm-Bayesian hybrid inference method: First, an initial parameter space model is constructed, which includes the distribution characteristics of the following thermodynamic parameters: specific heat ratio (γ): 1.3-1.4, obeying 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. The particle swarm optimization algorithm is used to generate the optimized sampling point set, 200 particles are set, and the position update formula is: X(t+1) = X(t) + V(t+1); the speed update adopts 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)); where w is the inertia weight, the dynamic value is 0.5-0.9; c1, c2, c3 are acceleration coefficients, which are 2.0, 2.0, 1.0 respectively; r1, r2, r3 are 0-1 random numbers; Pbest is the individual historical optimal position; Gbest is the group historical optimal position; X(k) is the randomly selected high fitness particle position. After executing 200 iterations, the best 50 particles are selected as sampling points.

[0123] Construct a parameter-associated Bayesian network and calculate the degree of mutual dependence between parameters: I(X;Y) = ∑∑p(x, y)log(p(x, y) / (p(x)p(y))); determine the network topology: compression factor Z as the root node; temperature T and pressure P as child nodes of Z; specific heat ratio γ as the child node of T. Establish the following nonlinear association model: Z = f(T, P) = 1 + BP / RT +(CP / RT) 2 , where B and C are the parameters to be estimated. Dynamic relationship between specific heat ratio and temperature: γ = γ 0 + aT + bT 2 , where γ 0 , a, b are the parameters to be estimated. Use the Sequential Monte Carlo algorithm to update the parameter posterior distribution 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 number of effective samples: Neff = 1 / ∑(w_i 2 ), resampling is performed when Neff is less than 500. To quantify the uncertainty of parameter estimates, 95% credible intervals and distribution difference measures are calculated: KL(P||Q) = ∑P(x)log(P(x) / Q(x)), where P is the posterior distribution and Q is the prior distribution.

[0124] 2.3. Parameter sensitivity analysis under multiple operating conditions.

[0125] Sensitivity analysis method based on information entropy: The system operating state is divided into five 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. The interaction between parameters is analyzed and the interactive sensitivity index is calculated: 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 operating condition conversion process, the time window sliding technology and 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, which is 0.25. The experimental results show that under full-load charging conditions, the sensitivity index of the compression factor Z is the highest (0.42), followed by the specific heat ratio γ (0.31); while under partial-load discharge conditions, the sensitivity index of the specific heat ratio γ is the highest (0.45).

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

[0127] 3.1. Temperature field initialization based on sparse measurement points.

[0128] Using the data of 20 spatially distributed temperature measurement points and the geometric model of the gas storage reservoir (90m×30m×20m cave), the 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-xi|| is the Euclidean distance between the spatial point and the measurement point. The interpolation result is optimized through boundary conditions and heat conduction equations.

[0129] 3.2. Physics-guided deep learning temperature field refinement.

[0130] Thermal-neural network architecture: Constructs a set of physical constraint equations for the temperature field, including the three-dimensional heat conduction equation: ρCp(∏T / ∏t) = ▽·(k▽T) + q; where ρ is the density, which is about 1.2kg / m for air. 3 ; Cp is the specific heat capacity, about 1005 J / (kg·K); k is the thermal conductivity, about 0.026 W / (m·K); q is the heat source term. Boundary conditions: -k(∏T / ∏n) = h(T - Tenv); where n is the boundary normal direction; h is the heat transfer coefficient, about 5 W / (m 2 ·K); Tenv is the ambient temperature, which depends on the rock temperature and is about 15°C. Design a thermal-neural network architecture with inputs of spatial coordinates (x, y, z) and time t, and output of temperature T. The network contains 8 hidden layers, 128 neurons per layer, and uses the GELU activation function. Construct a 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, about 10,000; B is the number of boundary sampling points, about 2,000; α, β, γ are weight coefficients, α=1.0, β=0.1, γ=0.1 in the early stage of training; α=0.5, β=0.3, γ=0.2 in the late stage of training. Design a multi-level sampling strategy and 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 and takes a value of 0.2. The network is trained using a small batch gradient descent algorithm with a batch size of 1024, a training round number of 10,000, a learning rate of 0.001, and an adaptive sampling point cloud is regenerated every 500 rounds. After training, the temperature field is predicted 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%.

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

[0132] Calculate the temperature gradient along the vertical direction of the gas storage: ▽T_z(x, y, z) = ∏T(x, y, z) / ∏z. Create 100 evenly distributed vertical profiles and calculate the vertical temperature gradient of 100 equally spaced points on each profile. Use a multi-scale edge detection algorithm to identify the temperature gradient mutation area and calculate the second-order derivative of the gradient: ▽2 T_z(x,y,z) = ∏ 2 T(x,y,z) / ∏z 2 The most drastic temperature gradient change is located by the zero crossing point of the second-order derivative, and the stratification interface is confirmed by applying the threshold discrimination. The improved Richardson number (Ri) is calculated to evaluate the stratification stability: Ri = (g / T 0 )(∏T / ∏z) / [(∏U / ∏z) 2 +(∏V / ∏z) 2 ]; where g is the acceleration due to gravity, 9.8 m / s 2 ; T 0 is the reference temperature, taking the average temperature of the region; ∏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, the velocity field is estimated based on the pressure field and temperature field: U(x, y, z) = -k 1 (∏P / ∏x) / μ; V(x, y, z) = -k 1 (∏P / ∏y) / μ; where k 1 is the permeability coefficient, about 10 -12 m 2 ; μ is the dynamic viscosity, about 1.8×10 -5 Pa·s. Calculate buoyancy frequency: N 2 = -(g / ρ 0 )(∏ρ / ∏z). The energy transmission coefficient is introduced to quantify the barrier effect of stratification on energy transfer: τ = exp(-∫(N 2 / ω 2 -1) 1 / 2 dz), where ω is the characteristic frequency, which is taken as 0.01Hz. Construct the stratification characteristic index: 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 of the i-th stratification interface; h_i is the thickness of the 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 area; SSI is the stratification stability index; Ri_i is the Richardson number of the i-th area; V_i is the volume of the area; 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 that there is a strong and stable temperature stratification.

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

[0134] 4.1. System control volume division.

[0135] 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 areas according to the temperature field characteristics; expansion unit: divided into preheater, high-pressure section, medium-pressure section, and low-pressure section.

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

[0137] The local entropy balance equation and entropy generation mechanism decomposition. The 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. Express entropy as a function of temperature and pressure: s = s(T, P) = s 0 + cp·ln(T / T 0 ) - R·ln(P / P 0 );where s 0 is the reference state specific entropy, J / (kg·K); cp is the constant pressure specific heat, J / (kg·K); R is the gas constant, J / (kg·K); T 0 , P 0 are the reference state temperature and pressure. Decompose the total entropy generation rate into contributions from different physical mechanisms: Flow friction entropy generation: S*_gen,fr =∫[(μ / T)·(∏u_i / ∏x_j + ∏u_j / ∏x_i) 2 ]dV; Heat transfer entropy generation: S*_gen, ht = ∫[(k / T 2 )·(▽T) 2]dV; Mixed entropy generation: 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. Plot the actual process path of the system on the temperature-entropy (Ts) diagram and calculate the irreversible loss of the process: W_loss = ∫T·dS_gen; Use path deviation to quantify process irreversibility: Δ_path = ∫[|ds_actual / dt - ds_reversible / dt|]dt / ∫[ds_reversible / dt]dt; Create a multi-scale entropy generation distribution visualization diagram and design the entropy generation intensity index: γ_s =(S*_gen / V) / (S*_gen / V)_avg. Experimental results show that in the high-pressure section of the compressor, flow friction entropy generation accounts for the highest proportion (62%); in the gas storage, heat transfer entropy generation is dominant (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%.

[0138] 4.3. Analysis of available energy of gas storage considering temperature stratification.

[0139] Introduction of interlayer available energy transfer terms. 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, and the difference between the calculated parameters and the stratified model is: ΔT = |T_i - T_avg|; ΔP = |P_i - P_avg|; Δh = |h_i - h_avg|; Δs = |s_i - s_avg|. Calculate the physical available energy and chemical available energy 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 available energy of the i-th layer, in J; B_ch,i is the chemical available energy of the i-th layer, in 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 molar fraction of component j. Introduce 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, W / (m 2 ·K); A_i-j is the interlayer contact area, m 2;Pe_i-j is the Peclet number, which characterizes the relative strength 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,ij;B_loss = B_unif - B_strat;η_loss = B_loss / B_unif. Analyzing the relationship between the available energy loss rate and the stratification characteristics, a linear correlation was found: η_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 index of gas storage: η_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 stratified impact factor.

[0140] Step 5: Dynamic tracking and status assessment of energy flow at multiple time scales.

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

[0142] Multi-scale decomposition of energy flow based on wavelet transform. Extract key energy flow signals from the multidimensional 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 storage filling and deflation process. Select the complex-valued Morlet wavelet as the mother wavelet function: ψ(t) = π^(-1 / 4)·e^(iω 0 t)·e^(-t 2 / 2); where ω 0 is the center frequency, which takes the value of 6. Adaptive scale selection algorithm is used: 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, which takes a value of 0.3; S(t, a_base) is the wavelet coefficient at the base 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, and ψ* is the complex conjugate of the mother wavelet. Calculate the wavelet energy spectrum and extract the time-frequency features: 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 . The energy flow signal is decomposed 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 has a correlation coefficient of 0.82 with the heat exchanger heat flow rate, and the phase lag is about 15 seconds; in the low-frequency component, the gas storage pressure can have a correlation coefficient of 0.93 with the expander output power.

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

[0144] 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 fluctuation rate of steady-state operating parameters. Short-term indicators include: energy conversion efficiency: output / input energy ratio; reliability: through parameter deviation statistical measurement. Long-term indicators include: cumulative efficiency: average efficiency of long-term operation cycle; equipment health: calculated based on vibration and temperature characteristics.

[0145] 5.3. Dynamic classification and anomaly detection of operating status.

[0146] Using multi-time scale state evaluation indicators, unsupervised learning methods are used to identify the typical operating state mode of the system, and three types of states are constructed: normal state: all indicators are within the target range; sub-health state: some indicators are slightly deviated; abnormal state: key indicators are seriously deviated. The deviation between the current state and the typical mode is calculated by 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.

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

[0148] 6.1 Parameter sensitivity analysis for efficiency bottlenecks.

[0149] Combining the bottleneck analysis graph 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).

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

[0151] Construct a multi-objective optimization problem model, and the objective functions include: f 1 (x) = -η_sys(x), to maximize the system efficiency; f 2 (x) = S_gen(x), to minimize the entropy generation; f 3 (x) = SI(x), to weaken the temperature stratification; f 4 (x) = -r_resp(x), to maximize 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: v_i(t + 1) = w·v_i(t) + c 1 ·r 1 ·(pbest_i - x_i(t)) + c 2 ·r 2 ·(gbest - x_i(t)) + c 3 ·r 3·(leader_i - x_i(t)); x_i(t+1) = x_i(t) + v_i(t+1). Introduce 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 number of iterations; t_max is the maximum number of iterations. Design a chaotic perturbation operator to enhance the global search capability: 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 expansion 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 characteristic solutions, including single-objective optimal solution, comprehensive indicator optimal solution, compromise solution and inflection point solution. Perform robustness analysis on the selected characteristic solutions and design parameter perturbation experiments: x_pert = x*+ Δ·x*·(1-2·rand()); where x* is the characteristic solution; Δ is the perturbation intensity, which is 0.08. Evaluate the sensitivity of the objective function to perturbations: S_rob = (1 / N)·∑(|f(x_pert) - f(x*)| / f(x*)), and select the final optimization solution based on the multi-attribute decision-making method.

[0152] After adjusting the operating parameters according to the optimization suggestions of this embodiment, the system performance was improved: the overall operating efficiency of the system was improved by 4.3 percentage points, from the original 62.5% to 66.8%; the energy storage efficiency of the gas storage reservoir was improved by 6.2 percentage points, from 85.4% to 91.6%; the equipment life was extended by 15%, mainly due to the optimization of operating parameters to reduce thermal stress and mechanical stress; the operating cost was reduced by 8.7%, saving approximately RMB 4.2 million in electricity bills annually; the advance warning time for system abnormal conditions was extended from the original 10 minutes to 2.5 hours, effectively preventing potential failures.

[0153] This implementation case describes in detail the actual application process of a compressed air energy storage optimization operation method based on multi-parameter collaborative control, focusing on the calculation process and implementation details. Through particle swarm-Bayesian hybrid inference, thermal-neural network architecture, local entropy balance equation and entropy generation mechanism decomposition, introduction of inter-layer available energy transfer terms, multi-scale decomposition of energy flow based on wavelet transform, multi-objective optimization and adaptive weighting scheme and other technologies, the system operation efficiency is significantly improved, the equipment life is extended, the operation cost is reduced, and the invention goal is fully achieved.

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

Claims

1. A method for 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.

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 were 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 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, a gas storage available energy evaluation index system is constructed, the thermodynamic efficiency of each link in the system is calculated, the efficiency bottlenecks are identified, and a bottleneck analysis map is generated.

7. The method according to claim 6, 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.

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

9. The method according to claim 8, 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.

10. The method according to claim 1, characterized in that 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.

Images