Humidity control method for gas storage reservoir in large-scale compressed air energy storage power station

Abstract

The present invention discloses a method for controlling humidity in a gas storage reservoir of a large-scale compressed air energy storage power station, which belongs to the technical field of intelligent control of energy storage power stations. The method comprises: obtaining grid dispatch instructions and real-time operating data of the power station, using a pre-configured PCE proxy model to perform rolling time domain optimization, and calculating the instantaneous optimal humidity set point; obtaining sparse real measurement data, calculating the prediction error, and identifying the equipment degradation stage based on the statistical characteristics of the error, and generating adaptive update weights; constructing and solving an optimization problem with physical constraints and adaptive weights based on an offline-constructed PCE coefficient coupling constraint matrix containing the law of conservation of thermodynamics, so as to collaboratively correct the coefficients of the PCE proxy model; and issuing the optimal humidity set point for execution. The present invention realizes the adaptability and high fidelity of the control model throughout its entire life cycle, significantly improving the safety and economy of the operation of the energy storage power station.

Description

Technical Field

[0001] The present invention relates to a compressed air energy storage technology, in particular to a humidity control method for a gas storage reservoir in a large-scale compressed air energy storage power station. Background Art

[0002] Compressed air energy storage (CAES), a new energy storage technology with the potential for large-scale, long-term, and low-cost development, plays a vital role in the global acceleration of energy transition and the construction of a new power system dominated by renewable energy. It can effectively mitigate the intermittent and volatile nature of renewable energy sources such as wind and solar power, providing power grids with various services such as peak shaving, frequency regulation, and backup. It is a key supporting technology for ensuring the safe and stable operation of the power grid and enhancing its capacity to absorb renewable energy. In large-scale CAES power plants with a capacity of one million kilowatts, energy is released primarily through the expansion of high-pressure air in turbines (turboexpanders). During this process, the air temperature rapidly drops due to the Joule-Thomson effect, even dropping below freezing under certain operating conditions. If the compressed air in the storage facility is too humid, liquid water will precipitate on turbine blades, valves, and pipe walls, or even condense into solid ice. This not only severely corrodes equipment and shortens system life, but can also cause ice erosion or dynamic imbalance on high-speed turbine blades, leading to catastrophic safety accidents. Therefore, precise, reliable and economical control of the air humidity at the gas storage outlet is the lifeline and core technical difficulty to ensure the safe, efficient and long-term operation of a million-kilowatt compressed air energy storage power station, and is directly related to the return on investment and technical feasibility of the entire power station.

[0003] Currently, various technical approaches exist in engineering practice and academic research for humidity control in compressed air energy storage systems. In most operational projects, a conservative control strategy based on worst-case design is commonly employed. This strategy uses offline thermodynamic calculations to determine the lowest temperature that can occur throughout the system's operating range. Based on this, a fixed, extremely low dew point temperature (e.g., -60°C) is set as the control target for the dehumidification system. Simple on-off or PID (proportional-integral-derivative) control is used to ensure that humidity remains within an absolutely safe range at all times. In academic research, to enhance the sophistication of control, some researchers have begun to explore mechanism-based control models. These models typically use simplified thermodynamic equations (such as the ideal gas equation of state) to estimate the required dew point temperature online, enabling dynamic adjustment of the control setpoint. Furthermore, with the advancement of data science, some cutting-edge research has begun applying data-driven surrogate models, such as offline-trained artificial neural networks (ANNs) or support vector machines (SVMs), to directly predict critical temperatures based on real-time operating inputs, providing a basis for humidity control decisions. These methods together form the basis of current humidity control technology and provide a preliminary solution to ensure the basic safe operation of energy storage power stations.

[0004] However, these existing technologies still face numerous challenges in pursuing deeper applications with higher security, economy, and adaptability, such as the issue of physical consistency and reliability of data-driven models. Summary of the Invention

[0005] The purpose of the invention is to solve the problem of poor physical consistency of traditional data-driven models and provide a humidity control method for gas storage reservoirs in large-scale compressed air energy storage power stations.

[0006] Technical solution: A method for controlling humidity in a large-scale compressed air energy storage power station gas storage reservoir, comprising:

[0007] Obtain grid dispatch instructions and real-time power plant operating data; use a pre-configured PCE agent model to calculate the instantaneous optimal humidity set point and generate a model prediction output trajectory;

[0008] Acquire sparse real measurement data, combine it with the model prediction output trajectory, and use the pre-configured PCE coefficient coupling constraint matrix to correct the coefficients of the PCE proxy model to generate an updated PCE proxy model for the next cycle;

[0009] The instantaneous optimal humidity set point is converted into a dehumidification equipment control instruction and sent for execution, and the execution results are collected as part of the power plant real-time operating condition data for the next control cycle.

[0010] Beneficial effect: By constructing a coefficient collaborative correction mechanism under physical constraints and an adaptive update mechanism of degradation perception, the control model is adaptive and high-fidelity throughout its life cycle, significantly improving the safety and economy of energy storage power station operation. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 It is a flow chart of the present invention.

[0012] Figure 2 It is a flow chart of the present invention for correcting the coefficients of the PCE proxy model.

[0013] Figure 3 It is a flow chart of constructing the PCE coefficient coupling constraint matrix of the present invention.

[0014] Figure 4 It is a flow chart of the present invention for calculating the instantaneous optimal humidity set point. DETAILED DESCRIPTION

[0015] To make the purpose, technical solutions and advantages of the present invention clearer, the following Figures 1 to 4 The specific embodiments of the present invention are described in detail. It should be noted that, unless there is a conflict, the embodiments and features in the embodiments of this application can be combined with each other. The specific embodiments described herein are merely for explaining the present invention, not for limiting the present invention.

[0016] First, during the research process, the applicant discovered that traditional purely data-driven agent models (such as neural networks) are essentially black boxes. Their online correction process relies entirely on measurement data and lacks the inherent constraints of physical laws. When system operating conditions change or equipment degrades prematurely, model corrections may cause its internal parameters to deviate from physical reality, resulting in mathematically correct but physically absurd predictions (for example, predicted system efficiency exceeding 100% or violating the law of conservation of energy). Using such physically untrustworthy models for closed-loop safety control poses significant potential risks. Furthermore, existing adaptive methods (such as standard Kalman filtering) often fail to effectively distinguish the source of prediction errors when correcting the model—whether they arise from random sensor measurement noise or actual performance degradation caused by equipment wear and efficiency. This indiscriminate correction approach can cause the model to overreact to random noise when the system is healthy, resulting in unnecessary jitter in the control output. However, when equipment degrades slowly, the model may be sluggish in responding and unable to track changes in the actual state, ultimately leading to a continuous decline in control accuracy and economy over time.

[0017] To this end, the following technical solutions are provided: grid dispatch instructions and real-time power plant operating data are obtained, and a pre-configured PCE proxy model (Polynomial Chaos Expansion) is used for rolling time domain optimization to calculate the instantaneous optimal humidity set point. Sparse real-world measurement data is obtained, prediction errors are calculated, and equipment degradation stages are identified based on the statistical characteristics of the errors to generate adaptive update weights. Based on an offline-constructed PCE coefficient coupling constraint matrix that incorporates thermodynamic conservation laws, an optimization problem with physical constraints and adaptive weights is constructed and solved to collaboratively modify the coefficients of the PCE proxy model. The optimal humidity set point is then issued for execution.

[0018] By building an intelligent humidity control framework that deeply integrates physical mechanisms and data-driven, and has online self-evolution capabilities, we have abandoned the drawback of traditional control methods where the model gradually becomes disconnected from the real world. By correcting the closed loop through dual constraints, the control core (PCE agent model) can continuously and accurately reflect the real health status of the equipment. Specifically,

[0019] Through sensitivity analysis and mathematical derivation, macroscopic thermodynamic conservation laws are transformed into a direct, rigid set of linear constraint equations for the abstract mathematical coefficients within the PCE proxy model. This ensures that any online model corrections strictly adhere to the fundamental laws of the physical world, achieving coefficient-level coupling of physical constraints and resolving the physical consistency issue of data-driven models. Furthermore, an intelligent mechanism for diagnosing error sources is constructed. By analyzing the time rate of change of the prediction error variance, the device's healthy, slowly degrading, or accelerated degradation phase is identified. Based on this information, different correction weights are dynamically assigned to different model coefficients, achieving a transition from blind correction to precise and efficient self-adaptation.

[0020] Therefore, it is possible to minimize the parasitic energy consumption of the dehumidification system and improve the round-trip efficiency of the power station while ensuring the absolute safety of the gas storage reservoir. At the same time, the model can automatically adapt to the performance degradation caused by equipment aging, ensuring the optimality of the control strategy throughout the life cycle of the power station, greatly enhancing the long-term operational reliability and economic benefits of the system.

[0021] First, the main characters (symbols) used in the examples and their unique meanings in this invention are explained: PCE stands for Polynomial Chaos Expansion. KKT stands for Karush-Kuhn-Tucker, a necessary condition for the optimal solution of nonlinear programming. x is the PCE coefficient vector, a set of mathematical coefficients that represent the PCE agent model. new is the updated PCE coefficient vector. old is the PCE coefficient vector before updating. z is the prediction error vector, which refers to the difference between the model prediction output and the actual measurement data. Kmodified is the modified Kalman gain matrix. A eq A is a linear equality constraint matrix that reflects the conservation law of thermodynamics. ineq is the linear inequality constraint matrix, defining the feasible region of the PCE coefficient. ineq is the boundary vector of the linear inequality constraint. L(·) is the Lagrangian function. λ is the Lagrangian multiplier vector corresponding to the equality constraint. μ is the Lagrangian multiplier vector corresponding to the inequality constraint. w j The adaptive update weight assigned to the j-th PCE coefficient. σ 2 is the variance of the prediction error. 2 is the change in the forecast error variance. Δt is the time step when calculating the variance change rate. j is the degradation sensitivity index of the jth PCE coefficient. α, β, γ are the preset adjustment parameters in the error-weight nonlinear mapping function. matrix is the filter state covariance matrix. P old is the state covariance matrix before updating. new is the updated state covariance matrix. constraint is the projection matrix. I is the identity matrix. C is the constraint Jacobian matrix. H is the observation Jacobian matrix. R matrix is the measurement noise covariance matrix. P is the pressure. T is the temperature. T dewpoint is the freezing point of the real gas.

[0022] w(P) is a smooth transition weight function that depends on the pressure P. T1 and T2 are independent freezing point temperature values ​​calculated by two different state equations.

[0023] Example 1: This embodiment provides a method for controlling humidity in a large-scale compressed air energy storage power station. This method can be applied to a computer device deployed in the main control system of the energy storage power station. The device includes a processor and a memory. The memory stores a computer program that, when executed by the processor, implements this method. The method specifically includes the following steps:

[0024] Step S100: Obtain grid dispatch instructions and real-time power station operating data; and use a pre-configured PCE agent model to calculate the instantaneous optimal humidity set point and generate a model prediction output trajectory.

[0025] The PCE proxy model can simulate the behavior of a high-fidelity physical model with extremely high computational efficiency, enabling rapid online prediction of the dynamic response of energy storage systems. The instantaneous optimal humidity setpoint is defined as the humidity control target value that ensures system operational safety (e.g., preventing icing in pipes) and minimizes operating costs (e.g., dehumidification energy consumption) within the current and a very short time step in the future. The model's predicted output trajectory represents the time series forecast of a series of key physical quantities (e.g., temperature and pressure) output by the PCE proxy model over a future prediction horizon (e.g., two hours).

[0026] Specifically, the power plant acquires real-time operating data, such as the current gas storage pressure P, temperature T, and gas flow rate, at a high frequency (e.g., once per second). Simultaneously, it obtains forecasted load instructions for the next N time steps (e.g., one value every 15 minutes for the next four hours) from the upper-level power grid dispatch center. Using this data as input, the currently optimal PCE proxy model is invoked within a rolling optimization framework to predict future system state changes. By solving a constrained optimization problem with operational economy as the objective, a series of future optimal control sequences are derived, and the first value is extracted as the current instantaneous optimal humidity setpoint. Simultaneously, the future state predictions generated during this optimization process are saved as model prediction output trajectories for subsequent model correction steps.

[0027] Traditional humidity control based on fixed thresholds cannot adapt to changing operating conditions and equipment aging, resulting in poor economic efficiency. By introducing a PCE proxy model for rapid prediction and rolling optimization, a shift from passive response to active prediction is achieved. The motivation is to use future predictions to plan the optimal control strategy in advance, thereby minimizing the dehumidification system's energy consumption while ensuring absolute safety, thereby improving the economic efficiency of the entire energy storage power station.

[0028] In some embodiments, the PCE proxy model can be replaced with other efficient proxy models, such as those based on radial basis function (RBF) networks, Kriging models, or deep neural networks (DNNs). Besides the sequential quadratic programming (SQP) algorithm, the rolling optimization solver can also be other nonlinear programming algorithms such as the interior point method or the active set method.

[0029] Step S200 involves acquiring sparse real-world measurement data, combining it with the model's predicted output trajectory, and using a preconfigured PCE coefficient coupling constraint matrix to modify the coefficients of the PCE proxy model to generate an updated PCE proxy model for the next cycle. Sparse real-world measurement data refers to precise measurements collected at a low frequency (e.g., once a minute) from high-precision sensors installed at key locations on the equipment. This frequency is lower than the control frequency, but its accuracy is higher than that of conventional operating condition monitoring data. The PCE coefficient coupling constraint matrix is ​​a mathematical matrix constructed offline and is used to translate macroscopic physical conservation laws (such as conservation of energy and conservation of mass) into mathematical constraints on the relationships between abstract coefficients within the PCE model.

[0030] For example, sparse real-world measurement data such as temperature and pressure at key locations (such as turbine outlets) is obtained from the power plant's Supervisory Control and Data Acquisition (SCADA) system. The predicted value at the same time as the real-world measurement data is extracted from the model prediction output trajectory generated in step S100, and the difference between the two is calculated to form a prediction error vector z. Based on this prediction error vector z and using the physical rules defined by the PCE coefficient coupling constraint matrix as strong constraints, a constrained optimization problem is constructed. Solving this problem yields a set of modified PCE coefficient vectors x. new Finally, the PCE agent model is reconstructed using this new set of coefficients for calculation in the next control cycle.

[0031] Because proxy models tend to lose accuracy over time and with device aging, this step, unlike traditional correction methods, establishes an online self-correction mechanism and introduces a PCE coefficient coupling constraint matrix. This ensures that the model correction process does not violate fundamental physical laws, avoiding mathematically correct but physically absurd correction results. This makes model correction more robust and reliable, accurately reflecting the true health status of the device and achieving full lifecycle adaptability.

[0032] Step S300: convert the instantaneous optimal humidity set point into a dehumidification equipment control instruction and issue it for execution, and collect the execution result as part of the power station real-time operating condition data of the next control cycle.

[0033] Specifically, the instantaneous optimal humidity set point calculated in step S100 (e.g., a specific dew point temperature) is converted into specific operational instructions that can be recognized and executed by the underlying dehumidification equipment. For example, these instructions may include turning a dehumidifier unit on or off, setting the target temperature for a regeneration heater, or adjusting a desiccant switching valve. These instructions are sent via an industrial control bus (e.g., Profibus or Modbus) to the dehumidification system's programmable logic controller (PLC) or distributed control system (DCS) for execution. Simultaneously, the system's execution status and actual changes in the gas storage reservoir's humidity are monitored. This collected data serves as feedback and becomes part of the power plant's real-time operating data for the next control cycle (i.e., the next execution of step S100), thus forming a closed-loop control system.

[0034] This step bridges the gap between upper-level optimization decisions and underlying physical execution, forming the execution and feedback loops in control theory. By accurately translating optimization results into device actions and monitoring execution results in real time, the stability and effectiveness of the entire closed-loop control system are ensured.

[0035] Through the above steps, this embodiment can adaptively modify the model online and perform predictive optimization control based on it, thereby significantly improving the operating economy of the compressed air energy storage power station while ensuring the absolute safety of the gas storage reservoir.

[0036] Embodiment 2: Based on the above embodiment 1, step S200, ie, the process of correcting the PCE proxy model coefficients, is described in more detail.

[0037] First, the steps of correcting the coefficients of the PCE proxy model include: determining the prediction error vector based on the difference between the model's predicted output trajectory and the sparse real measurement data; constructing a constrained optimization problem with the PCE coefficient coupling constraint matrix as the physical constraint condition and the prediction error vector as the optimization basis; solving the constrained optimization problem to produce an updated PCE coefficient vector that satisfies the physical coupling constraints; and reconstructing the updated PCE proxy model using the updated PCE coefficient vector.

[0038] In this embodiment, the process is further optimized, that is, before constructing the constrained optimization problem, it also includes: identifying the current degradation stage of the device based on the statistical characteristics of the prediction error vector within the sliding time window, especially the rate of change of its variance; and, in response to the identified degradation stage, generating a corresponding adaptive update weight vector for each PCE coefficient, and applying the weight vector to the construction of the constrained optimization problem.

[0039] The specific implementation of this process is as follows:

[0040] Step S210: Identify the current degradation stage of the device. This step includes calculating the variance of the prediction error vector within a sliding time window and determining the rate of change of the variance between consecutive time windows to obtain an error variance change rate; comparing the error variance change rate with at least one preset threshold, and classifying the degradation stage into one of a plurality of predefined degradation modes based on the comparison result.

[0041] For example, the process is as follows:

[0042] The system maintains a sliding time window of length L (e.g., storing the forecast error vector z for the past 60 minutes).

[0043] Calculate the variance σ of the prediction error of each physical quantity (such as temperature and pressure) in the window 2 . Further, calculate the variance σ of the current time window 2 The variance σ between (t) and the previous time window 2 The rate of change between (t-Δt), that is, the variance change rate Δσ 2 (t)=[σ 2 (t)-σ 2 (t-Δt)] / Δt.

[0044] The calculated variance change rate Δσ 2 (t) is compared with two preset thresholds (threshold 1, threshold 2).

[0045] If |Δσ 2 (t)∣[cite start ]<threshold 1, the device is determined to be in the initial stable stage.

[0046] If the threshold 1≤|Δσ 2 (t)∣[cite start ]<threshold 2, it is determined to be a mid-term slow degradation stage.

[0047] If |Δσ 2 (t)∣[cite start ]≥threshold 2, it is determined to be a late accelerated degradation stage.

[0048] Simply modifying the model cannot distinguish whether errors are caused by random noise or true device degradation. This step analyzes the changing trends in the statistical characteristics of the errors (i.e., the rate of change of the variance) to determine the evolutionary pattern of the device's health. This makes model modifications more targeted, allowing the system to respond promptly and make more significant adjustments to the model when accelerated device degradation occurs.

[0049] Step S220: Generate corresponding adaptive update weight vectors for each PCE coefficient. Specifically, the steps of generating corresponding adaptive update weight vectors for each PCE coefficient include: quantifying the degradation sensitivity index of each PCE coefficient based on the identified degradation stage; using a preset non - linear mapping function to calculate the weight value for each PCE coefficient according to the change rate of the degradation sensitivity index and variance; and normalizing the weight values to form an adaptive update weight vector.

[0050] In this embodiment, the specific process is as follows:

[0051] Based on the degradation stage identified in step S210 and combined with the coefficient - physical quantity mapping table constructed offline, calculate the degradation sensitivity index S j of each PCE coefficient x sensitivity_j . This index is used to quantify the degree of influence of the change of a certain coefficient on the physical quantity that is most significant for the current degradation.

[0052] Adopt a preset non - linear arctangent hyperbolic function (tanh) to calculate the original weight value: w j =α·tanh(β·S sensitivity_j ·∣Δσ j 2 ∣)+γ. Where α, β, and γ are preset adjustment parameters.

[0053] Normalization: Normalize all the calculated original weight values to ensure that the sum of all weights is a constant, and form the final adaptive update weight vector.

[0054] Principle / motivation / effect elaboration: The traditional filtering algorithm has the same correction strength for all states (coefficients). The motivation of this step is to achieve differential and intelligent correction. For those coefficients that are highly correlated with the current device degradation phenomenon (high sensitivity) and have a large prediction error (large variance change rate), higher update weights are assigned, enabling the model correction to be targeted, thus improving the efficiency and accuracy of the correction.

[0055] Step S230: Solve a constrained optimization problem to update the PCE coefficient. Specifically, it includes: for the constrained optimization problem, introducing Lagrange multipliers to construct a Lagrangian function; based on the Lagrangian function, establishing a set of KKT condition equations representing the necessary conditions for the optimal solution; and calculating the updated PCE coefficient vector by iteratively solving the KKT condition equations. In this embodiment, the specific process is as follows:

[0056] Construct the following constrained optimization problem:

[0057] Minimize the objective function: min丨丨x new -(x old +K modified ·z)丨丨2 ;

[0058] The constraint conditions are: A eq ·x new = 0 and A ineq ·x new ≤ b ineq where K modified is the modified Kalman gain considering the adaptively updated weight vector calculated in step S220.

[0059] Introduce Lagrange multipliers λ and μ for the above problem, and construct the Lagrangian function L(x new , λ, μ) = ||x new -(x old + K modified ·z) || 2 + λT(A eq x new ) + μT(A ineq x new - b ineq ).

[0060] By taking the partial derivatives of the Lagrangian function and setting them to zero, a set of KKT condition equations is established. Use iterative algorithms such as sequential quadratic programming (SQP) to solve this system of equations to obtain the updated PCE coefficient vector xnew that finally satisfies all physical constraints.

[0061] This step makes the update of the PCE coefficients carried out within the framework of strictly adhering to physical laws, and its solution is the one closest to the Kalman filter prediction value among all feasible solutions that satisfy physical constraints, taking into account both data-driven corrections and physical rule constraints.

[0062] Step S24: Update the filter state covariance matrix. In this embodiment, modifying the coefficients of the PCE surrogate model also includes updating the filter state covariance matrix, specifically: based on the PCE coefficient coupling constraint matrix, construct a projection matrix that can project uncertainties onto the physical constraint subspace; use the projection matrix to modify the standard covariance update process to generate an updated state covariance matrix that reflects the influence of physical coupling constraints.

[0063] In this embodiment, the specific process is as follows:

[0064] Based on the constraint matrices A eq and A ineq combine them into the constraint Jacobian matrix C. Construct the projection matrix P constraint = I - C T (CC T ) -1 C

[0065] Use this projection matrix to modify the covariance update formula of the standard Kalman filter: P new =P constraint ·(IK modified H)·P old ·P T constraint .

[0066] Traditional covariance updates fail to consider the physical constraints between state variables (PCE coefficients), leading to an overestimation of their uncertainty. This step, through the projection matrix, constrains the updated uncertainty (covariance) to a physically valid space, resulting in a more accurate estimate of the model state and providing more reliable prior information for the next round of corrections.

[0067] Example 3: This example describes in detail the offline construction method of the PCE coefficient coupling constraint matrix, which is the prerequisite for implementing constraint correction in Example 2. This method can be completed during the energy storage power station design or model initialization stage.

[0068] Specifically, the preconfigured PCE coefficient coupling constraint matrix is ​​generated through an offline construction method, including: performing sensitivity analysis on the coefficients and basis functions of the initial PCE proxy model, identifying the physical quantities dominated by each PCE coefficient, and establishing a coefficient-physical quantity mapping relationship; based on the coefficient-physical quantity mapping relationship and the thermodynamic conservation law of the system, deriving the mathematical constraint relationship that must be satisfied between the PCE coefficients; and formally constructing the mathematical constraint relationship into a PCE coefficient coupling constraint matrix.

[0069] Among them, the mathematical constraints include: a set of equality constraints derived from the conservation laws of thermodynamics and reflected as linear combinations of PCE coefficients; and a set of inequality constraints derived from the physical boundary conditions of the system and reflected as the range of values ​​that the PCE coefficients must satisfy.

[0070] In this embodiment, the specific calculation process is as follows:

[0071] Step S310: Identifying the Mapping Relationship Between PCE Coefficients and Physical Quantities. After the initial PCE proxy model training is complete, a sensitivity analysis is performed on its coefficients x and basis functions. Specifically, by perturbing individual PCE coefficients one by one, the model output physical quantity (such as turbine efficiency, heat transfer coefficient, pressure loss coefficient, etc.) is observed to see which changes most significantly. A mapping table is then established between each (or each group of) PCE coefficients and a dominant physical quantity. For example, the analysis revealed that coefficients x5 and x6 primarily affect turbine efficiency, while coefficient x8 primarily affects pipeline pressure loss.

[0072] Step S320: Thermodynamic conservation constraint relationships are constructed based on the mapping relationship established in step S310, converting macroscopic physical laws into constraints between coefficients.

[0073] Equality constraints: For example, consider the energy conservation law for a system: input energy = output work + heat loss. Since the mapping table already knows which coefficients represent efficiency and which represent heat loss, this energy conservation equation can be converted into one or more linear combinations of PCE coefficients that must satisfy the following equation: for example, c1x5 + c2x6 - c3x8 = 0. Combining these relationships creates a system of linear equality constraints.

[0074] Inequality constraints: The physical boundaries or inherent properties of a system constitute inequality constraints. For example, the PCE coefficient (or combination thereof) corresponding to any efficiency term (such as turbine efficiency) must be greater than 0 and less than 1 (or less than the Carnot cycle efficiency limit) under all operating conditions. Similarly, the pressure loss coefficient must be greater than or equal to 0. These constitute inequality constraints.

[0075] Step S330: Mathematical construction of coupling constraint matrix to formalize the mathematical constraint relationship derived in step S320.

[0076] Write all linear equality constraint equations into matrix form A eq x=0, thus obtaining the equality constraint matrix A eq .

[0077] Write all inequality constraints in standard form A ineq x≤b ineq , thus obtaining the inequality constraint matrix A ineq and the boundary vector b ineq .

[0078] Finally, the complete PCE coefficient coupling constraint matrix is ​​composed of Aeq, A ineq and b ineq The matrix is ​​constructed once and for all in the offline phase and stored in the system for online correction.

[0079] Example 4: This example describes in detail how to calculate the instantaneous optimal humidity set point.

[0080] Specifically, the step of calculating the instantaneous optimal humidity set point is implemented through rolling horizon optimization. This process includes: at each prediction step of the rolling horizon optimization, based on the current pressure value, adaptively selecting the optimal real gas state equation from a preset pressure-state equation mapping table; using the selected real gas state equation to calculate the corresponding real gas freezing point temperature; using the real gas freezing point temperature as a safety constraint, solving the rolling horizon optimization to obtain the instantaneous optimal humidity set point.

[0081] To solve the calculation jump problem that may exist in the above method, if the current pressure value is in the boundary area of ​​the adjacent pressure interval, the corresponding real gas freezing point temperature is calculated using the selected real gas state equation, which is achieved through a smooth transition algorithm, including:

[0082] Two independent freezing point temperature values ​​are calculated using two real gas state equations corresponding to adjacent pressure intervals. A weight coefficient is determined based on the relative position of the current pressure value within the boundary area. The weight coefficient is used to perform a weighted average of the two independent freezing point temperature values ​​to obtain the final real gas freezing point temperature after a smooth transition.

[0083] In this embodiment, the calculation process is as follows:

[0084] Step S410: Pressure interval division and state equation mapping In the offline phase, the operating pressure range is divided into multiple intervals according to the design pressure of the gas storage and the physical properties of the gas. For example:

[0085] Low pressure section (<1MPa): The gas behavior is close to that of an ideal gas, and the ideal gas state equation is selected for its simple calculation.

[0086] Medium pressure section (1-8MPa): The gas non-ideality is significant, so the Redlich-Kwong (RK) equation with higher accuracy is selected.

[0087] High-pressure section (>8MPa): The gas is in a high-density state, so the Peng-Robinson (PR) equation, which is suitable for high-pressure environments, is selected. The mapping relationship between this pressure and state equation is stored as a lookup table.

[0088] Step S420: Adaptive state prediction within the rolling time domain is performed at each prediction time step (for example, a step length of 1 minute and a prediction time domain of 2 hours) of the online rolling optimization calculation. For example, the following is an example:

[0089] Get the predicted pressure value P for the current time step.

[0090] The mapping table of step S410 is queried to determine the optimal real gas state equation corresponding to the pressure value P.

[0091] Call the state equation to calculate the actual gas freezing point temperature T at the pressure P and predicted humidity dewpoint .

[0092] At the same time, the PCE proxy model predicts the turbine exhaust temperature T at this time step exhaust .

[0093] Step S430: Smooth transition processing In step S420, if the predicted pressure value P falls exactly within the boundary area between two adjacent pressure intervals (for example, within the range of ±0.1 MPa near the boundary of the medium pressure section and the high pressure section of 8 MPa), the smooth transition algorithm is started, as follows:

[0094] The RK equation and the PR equation are called at the same time to calculate two independent freezing point temperature values, T1=TRK and T2=TPR.

[0095] Design a weight function w(P) that varies linearly with pressure P. For example, in the range of 7.9 MPa to 8.1 MPa, when P = 7.9 MPa, w(P) = 1, and when P = 8.1 MPa, w(P) = 0.

[0096] The final freezing point temperature is obtained by weighted average: T dewpoint =w(P)·T1+(1-w(P))·T2.

[0097] This embodiment avoids a step jump in the freezing point temperature calculation result caused by a hard switch of the state equation when crossing the pressure range boundary, thereby ensuring the continuity of the constraint boundary and the stability of the optimization solution.

[0098] Step S440: constrained optimization solution is used to construct a function with the minimum total dehumidification energy consumption (or the weighted sum related to adsorbent consumption and system pressure loss) in the predicted time domain as the optimization target. exhaust Must be higher than the corresponding real gas freezing point temperature T dewpoint A safety margin is added as a hard constraint. The sequential quadratic programming (SQP) algorithm is used to solve the constrained optimization problem, the optimal control sequence is obtained, and the first value is extracted as the instantaneous optimal humidity set point.

[0099] Through the above steps, this embodiment can accurately calculate the optimal humidity set point while ensuring absolute safety, and ensure the stability and reliability of the control process through a smooth transition algorithm.

[0100] Embodiment 5: Based on embodiment 3, this embodiment describes the process of offline construction and publishing of an adaptive physical agent model to demonstrate how to complete the entire offline preparation work.

[0101] Step S500: Offline construction and publishing of the adaptive physical agent model. This step is completed on a server with high-performance computing capabilities, providing the basic model and key parameters for subsequent online control. Specifically, this step includes the following sub-steps:

[0102] Step S510: Acquire basic data and define the operation space.

[0103] Equipment performance graph data: refers to a table or set of curves provided by equipment manufacturers (such as turbine and compressor manufacturers) that describes the performance parameters (such as efficiency and power consumption) of their products under different operating conditions (such as different speeds, pressure ratios, and flow rates).

[0104] Operational state space: A multi-dimensional vector space whose dimensions correspond to key input parameters that affect the system's operating state (such as ambient temperature, initial pressure, load instructions, etc.). The boundaries of the space are determined by the system's design limits and operating procedures.

[0105] First, complete design parameters for the target million-kilowatt compressed air energy storage power plant are collected, such as the effective volume of the gas storage reservoir (e.g., 300,000 cubic meters), geometric parameters, pipeline materials and dimensions, turbine / compressor ratings, and heat exchanger design parameters. Performance data for key equipment is also imported, such as an efficiency-operating condition mapping table provided by a turbine manufacturer for a specific model. Based on these parameters and the power plant's operating procedures (e.g., maximum and minimum operating pressures, maximum load rate), an operating state space covering all possible operating conditions is mathematically defined.

[0106] Step S520: Construct a high-fidelity training dataset.

[0107] Specific implementation method: Within the operating state space defined in step S510, to avoid data point clustering, a Latin hypercube sampling method is used to generate tens of thousands (e.g., 50,000) evenly distributed operating condition input vectors. Each input vector represents a possible operating condition. Subsequently, a high-fidelity physical model (e.g., a mechanistic model based on AspenPlus or Modelica) that couples the device performance map with a high-precision real-world gas equation of state (e.g., the PR equation) is invoked to solve each of these 50,000 operating condition input vectors, obtaining high-precision system output results (e.g., turbine exhaust temperature, pressure, system efficiency, etc.) for each operating condition. Finally, the input operating condition vectors and the corresponding output results are combined to form a large offline training dataset containing 50,000 samples.

[0108] Step S530: training an initial PCE proxy model.

[0109] Specific implementation method: Based on the high-fidelity training data set generated in step S520, a polynomial chaos expansion (PCE) method is used to train the proxy model.

[0110] Basis function selection: Analyze the statistical distribution characteristics of each input parameter in the training dataset. If the parameters are approximately Gaussian, choose Hermite polynomials as the basis function; if they are approximately uniformly distributed, choose Legendre polynomials as the basis function.

[0111] Coefficient solution: A non-invasive method is used to solve a set of optimal PCE coefficients that minimize the error between the PCE model output and the high-fidelity output in the training dataset through the least squares method enhanced by QR decomposition.

[0112] Model solidification: The optimal orthogonal polynomial basis function is combined with the solved optimal coefficient vector x and solidified to form the initial PCE proxy model. This model can complete a prediction calculation in milliseconds.

[0113] Step S540: Constructing a physical coupling constraint matrix for PCE coefficients. This step is a key innovative preparation that provides a physical rule basis for subsequent online corrections.

[0114] Sub-step S541: Identify the mapping relationship between PCE coefficients and physical quantities. Read the optimal PCE model coefficients and basis functions output from step S530. Combined with the thermodynamic parameters of the energy storage system, sensitivity analysis is performed to identify which macroscopic physical quantity each PCE coefficient has the greatest impact on. For example, by perturbing the jth coefficient xj and observing the change in the model output, if it is found that it primarily affects turbine efficiency, a mapping relationship from xj to turbine efficiency is established. Ultimately, a complete coefficient-physical quantity mapping table is generated.

[0115] Sub-step S542: Construction of thermodynamic conservation constraint relationship. Based on the mapping table of sub-step S541 and the energy conservation and mass conservation laws of the system, the mathematical relationship that must be satisfied between the PCE coefficients is derived. For example, based on the law of conservation of energy ▽·(ρvh)=▽·(k▽T), it is reflected as a linear relationship in the lumped parameter model. If the coefficient xa dominates the efficiency and the coefficient xb dominates the heat loss, a set of linear constraint equations of the form kaxa+kbxb=constant can be derived. At the same time, nonlinear physical constraints are identified, such as the efficiency must be less than the theoretical upper limit determined by the Carnot cycle efficiency, which constitutes a set of nonlinear constraint functions.

[0116] Sub-step S543: Mathematical construction of the coupling constraint matrix. Formalize the derived mathematical relations. Convert all linear constraint equations into the standard form of equality constraint matrix Aeq and vector beq, ensuring that Aeq·x=beq (in some cases, the homogeneous linear equation system can be expressed as Aeq·x=0). Convert the set of nonlinear constraint functions (or their linearized forms) into the inequality constraint matrix A ineq and constraint boundary vector b ineq , ensure A ineq x≤b ineq Finally, Aeq, beq, A ineq 、b ineq Combine to form a complete PCE coefficient coupling constraint matrix and store it for future use.

[0117] Sub-step S544: Verify the validity of the constraint matrix and initialize the filter state. Use the initial optimal PCE coefficient xinitial obtained in step S530 to verify and ensure that it meets all constraints, i.e. Aeq·x initial =beqandA ineq ·x initial ≤b ineq , to verify the mathematical validity of the constraint matrix. After verification, the initial coefficient vector x initial Set as the initial state estimation vector of the Kalman filter during subsequent online correction. At the same time, according to the uncertainty of the physical quantity corresponding to each PCE coefficient, set the diagonal filter initial state covariance matrix P matrix , where a larger initial variance is assigned to the high-order coefficients corresponding to those physical quantities that are easily affected by equipment aging (such as efficiency and heat transfer coefficient).

[0118] Through the above steps, a computationally efficient, physically constrained initial PCE proxy model and its correction system were fully constructed and released, fully preparing for subsequent online applications.

[0119] Example 6. Based on Examples 1 and 4, this example provides a detailed description of the complete closed-loop process of online rolling optimization and instruction execution.

[0120] Step S600: Online rolling optimization and instantaneous optimal humidity set point calculation. This step is periodically executed by the main control computer when the energy storage power station is in operation.

[0121] Sub-step S610: Real-time Data Acquisition and Preprocessing. Real-time power plant operating data, including key parameters such as gas storage pressure, temperature, and gas flow, is collected via the fieldbus at a high frequency (e.g., once per second). Simultaneously, forecasted load instructions for the next N time steps (e.g., the next two hours) are received from the power grid dispatching center. The collected high-frequency data is filtered (e.g., using a median filter or low-pass filter) and de-noised to obtain stable and reliable pre-processed operating data, which serves as the current state input for the optimization calculation.

[0122] Sub-step S620: predicting the adaptive state of the voltage division section.

[0123] Constructing pressure ranges and mapping tables: Based on the design parameters of the gas storage facility, the operating pressure is divided into three ranges: low pressure (<1MPa), medium pressure (1-8MPa), and high pressure (>8MPa). For each range, the optimal real gas equation of state is matched and a mapping table is constructed: the low pressure range corresponds to the ideal gas equation, the medium pressure range corresponds to the RK equation, and the high pressure range corresponds to the PR equation.

[0124] Switching smooth transition algorithm: To solve the problem of calculation result jump when switching across pressure intervals, a smooth transition algorithm is designed. For example, within the range of ±0.1MPa of the pressure interval boundary, a weighted average method is used. The weight function w(P) changes linearly with the pressure P in this area, and the final freezing point temperature T dewpoint =w(P)*T 方程1 +(1-w(P))*T 方程2 .

[0125] Calculate the state over the forecast horizon step by step: Set the forecast horizon for the rolling optimization (e.g., 2 hours, 1-minute steps). Starting from the current moment, perform a step-by-step forecast based on the forecasted load command. At each time step, the system queries the mapping table and weighting function to determine the applicable state equation based on the predicted pressure value for that step. The system then calls the latest adaptive PCE proxy model to calculate the predicted output for that step (e.g., turbine exhaust temperature and pressure). The corresponding real gas freezing point temperature is then calculated using the selected state equation.

[0126] Prediction trajectory verification and output: Combine the results of all prediction time steps to form a future state prediction trajectory that includes a series of predicted exhaust temperature and a series of actual gas freezing point temperatures. Check the integrity and physical plausibility of this trajectory data (e.g., whether the temperature and pressure trends conform to operating rules). Once verified, output to the next step.

[0127] Sub-step S630: Constrained optimization solves the optimal control strategy. An optimization objective function is constructed to minimize the weighted sum of adsorbent consumption and system pressure loss over the predicted time domain. Hard constraints are established to ensure that at any time during the predicted future state trajectory, the predicted exhaust temperature is always higher than the corresponding actual gas freezing point (a safety margin may be added). A sequential quadratic programming (SQP) algorithm is used to solve this nonlinear constrained optimization problem, resulting in an optimal humidity control sequence. Finally, only the control value at the first time step in the sequence is extracted as the current instantaneous optimal humidity setpoint.

[0128] Step S640: Humidity control command execution and feedback. This step completes the closed loop of control.

[0129] Sub-step S641: Control Instruction Conversion and Distribution. The instantaneous optimal humidity setpoint (a logical value) calculated in step S630 is converted into specific, executable control instructions for the underlying dehumidification equipment. For example, if the setpoint is -40°C dew point, this will cause dehumidifier unit A to start and the regeneration temperature to be set to 150°C. These instructions are then distributed via the control bus to the dehumidification system's PLC or DCS actuators.

[0130] Sub-step S642: Execute performance monitoring and feedback. Sensors continuously monitor the actual operating status of the dehumidification equipment and the actual changes in the humidity in the gas storage. This monitored and collected data serves as control feedback and is retrieved and processed as part of the power plant's real-time operating data at the start of the next control cycle (usually after one minute), forming a complete closed-loop control system.

[0131] Embodiment 7: This embodiment describes a collaborative PCE coefficient correction process based on coupling constraints.

[0132] Step S700: Cooperative PCE coefficient modification based on coupling constraints. This step is performed once in a specific time period (eg, every minute) to update the PCE agent model to adapt to changes in device status.

[0133] Sub-step S710: Acquisition of sparse measurement data and error calculation. Sparse real-world measurement data at key locations is acquired from the power plant's measurement system (e.g., high-precision thermocouples and pressure transmitters) at a low frequency (1 minute). For example, the accurately measured value of the turbine outlet temperature is 35.2°C. Simultaneously, the model-predicted value at the same moment in time is extracted from the future state prediction trajectory generated in step S620, for example, a predicted temperature of 35.8°C. The difference between the two values ​​is calculated to obtain the prediction error vector z = [0.6, ...] for that moment in time.

[0134] Sub-step S720: Degradation stage identification and weight adaptive calculation.

[0135] Specifically, the current degradation stage of the device is identified based on the statistical characteristics of the prediction error vector within the sliding time window, including the rate of change of the variance;

[0136] In response to the identified degradation stage, a corresponding adaptive update weight vector is generated for each PCE coefficient and applied to the formulation of the constrained optimization problem.

[0137] In one embodiment, the process is as follows:

[0138] Maintain a sliding time window that stores the prediction error vector z for the past 60 minutes. Calculate the mean and variance σ of the prediction error sequence of each physical quantity (temperature, pressure, etc.) in the window. i 2 , skewness and other statistical characteristics.

[0139] Calculate the rate of change Δσ of the error variance of each physical quantity between consecutive time windows i 2 (t)=[σ i 2 (t)-σ i 2(t-Δt)] / Δt. According to the change rate matrix, the degradation stage of the current device is identified by comparing it with the preset threshold, that is, the degradation mode is identified. For example, if the temperature error variance change rate |Δσ temp 2 ∣>[cite start ] threshold 2, it is determined to be a late acceleration degradation stage.

[0140] Based on the identified degradation stage and the coefficient-physical quantity mapping table constructed in Example 5, the sensitivity index S of each PCE coefficient is calculated. sensitivity_j =∑ i (|dPCE coefficient j / d physical quantity i ∣×Degradation degree i ), and obtain the coefficient degradation sensitivity analysis results.

[0141] Construct and apply the nonlinear mapping function w j =α·tanh(β·S sensitivity_j ·|Δσ j 2 ∣)+γ. Substitute the sensitivity index and error change rate corresponding to each coefficient, calculate the original weight value, normalize it, and finally output the adaptive updated weight vector, that is, the error-weight mapping.

[0142] Sub-step S730: Correction of the synergy coefficient under coupling constraints.

[0143] Modified Kalman gain construction: The adaptive update weight vector obtained in sub-step S720 is constructed into a diagonal weight matrix W adapt . Calculate the Jacobian matrix H of the model output about the PCE coefficient. Modify the traditional Kalman gain formula: K modified =W adapt ·P old ·H T ·(H·Pold·H T +R matrix ) -1 .

[0144] Constrained optimization problem construction: Read the PCE coefficient coupling constraint matrix (Aeq, Aeq) constructed in Example 5 ineq , b ineq ), combined with the modified Kalman gain and the prediction error vector, the following constrained optimization problem is constructed:

[0145] Minimize |||xnew-(xold+Kmodified·z)|| 2 ;

[0146] Constrained to A eq ·x new =beqandAineq ·x new ≤b ineq ;

[0147] Lagrange multiplier method solution: The Lagrange multiplier method is used to solve this problem and construct the Lagrange function L(x new ,λ,μ). Establish its KKT conditional equations and use the sequential quadratic programming (SQP) algorithm to iteratively solve them to obtain the updated PCE coefficient vector x that satisfies all physical constraints. new .

[0148] Covariance matrix constraint update: construct the constraint Jacobian matrix C and projection matrix P constraint =IC T (CC T ) -1 C. Update formula P using projection covariance new =P constraint ·(IK modified H)·P old ·P T constraint To update the covariance matrix, we ensure that the updated uncertainty also reflects the influence of physical constraints.

[0149] Sub-step S740: Model reconstruction and release. Update the PCE coefficient vector x obtained in sub-step S730 new The updated adaptive PCE proxy model is recombined with the original orthogonal polynomial basis function to construct the updated adaptive PCE proxy model. The new model is quickly verified for physical consistency (such as energy conservation error, efficiency rationality) and numerical stability (such as condition number). After verification, the updated model is combined with the updated filter state covariance matrix P new Released together for the next round of online optimization and prediction.

[0150] Example 8: Systematic determination method of key parameters.

[0151] According to one aspect of the present application, a method for determining a degradation identification threshold is as follows:

[0152] The threshold is set to distinguish normal measurement noise fluctuations from real error trends caused by device performance degradation.

[0153] Using an offline simulation platform, the established high-fidelity model of the energy storage system was subjected to thousands of Monte Carlo simulations. Some simulations were run under healthy operating conditions, while others were simulated with faults of varying severity and rate (e.g., linear decrease in turbine efficiency over time).

[0154] The specific process is as follows: 1) Record the rate of change of the prediction error variance Δσ in all simulations 22) Plot Δσ 2 Probability density distribution diagram of . 3) Δσ under healthy working conditions 2 The distribution will be concentrated in a smaller range, and the Δσ under the fault condition 2 4) Select Δσ that can distinguish the two distributions with 95% confidence. 2 The value is used as the judgment threshold for the mid-term slow degradation stage (threshold 1). Similarly, the inflection point value that can distinguish slow degradation from accelerated degradation is selected as threshold 2.

[0155] According to one aspect of the present application, a method for tuning the nonlinear mapping function parameters α, β, and γ is specifically as follows:

[0156] These three parameters together determine the sensitivity and magnitude of weight updates. α controls the base weight, β controls the sensitivity to error changes, and γ controls the upper limit of weight saturation. A grid search and cross-validation approach based on a validation dataset is used.

[0157] The specific steps are as follows: 1) Separate an independent validation dataset from the offline training dataset. 2) Set the possible value ranges and step sizes for α, β, and γ, for example, α∈[0.5, 1.0, 1.5], β∈[1, 5, 10], and γ∈[1, 2, 3]. 3) Iterate through all parameter combinations. Run a complete model correction and prediction simulation for each combination. 4) Minimize the long-term root mean square error (RMSE) of the predictions on the validation dataset as the optimization objective to find the optimal α, β, and γ parameter combination.

[0158] According to one aspect of the present application, the optimal selection of the sliding time window size L is specifically as follows:

[0159] The window size L represents a trade-off between response speed and statistical stability. If the window is too small, the variance calculation is susceptible to noise, leading to misjudgment; if the window is too large, the response to emerging degradation trends will be slow.

[0160] The specific process is: 1) Obtain a long error time series that includes a typical degradation process. 2) Calculate the autocorrelation function of this series. 3) Find the time delay τ when the autocorrelation function first drops to a low value (such as 0.1). 4) In theory, choosing L to be approximately between 2τ and 3τ ensures statistical independence while obtaining sufficient samples for a stable variance estimate.

[0161] Example 9: System exception handling and safety assurance mechanism.

[0162] Scenario 1: Failure in solving the constrained optimization problem.

[0163] In step S730 , the SQP solver fails to converge within a preset maximum number of iterations (eg, 200), or returns a flag indicating no solution or infeasibility.

[0164] 1) Immediate measures: The system abandons the coefficient update and the PCE proxy model maintains the state of the previous cycle. 2) Control instructions: To ensure safety, the controller does not use the optimization results that may be risky, but instead executes a predefined safety and conservative instruction, for example, temporarily adjusting the humidity set point to an absolute safety value 5°C lower than the historical lowest temperature. 3) Alarm and recording: The system immediately generates a level 2 alarm indicating that the model correction has failed, notifying the operator to pay attention. At the same time, all input data before the failure (such as x old ,z,K modified 4) Recovery Mechanism: The system continuously attempts to find a solution over the next few control cycles. If failures exceed three consecutive times, the alarm level is raised to level 1, the adaptive correction module is automatically suspended, and the system switches to non-adaptive control mode based on the initial PCE model, awaiting manual intervention.

[0165] Scenario 2: Key sensor data is missing or faulty.

[0166] In step S710, the key real measurement data (such as turbine outlet temperature) used to calculate the prediction error vector z is not updated within a predetermined time, or its value exceeds a reasonable physical range (for example, the temperature is lower than -50°C or higher than 200°C), or its rate of change is abnormal.

[0167] 1) Data replacement: The system attempts to use the data from a backup sensor (if available). If there is no backup, the system will temporarily use the valid measurement value from the previous cycle in a short period of time (e.g., 1-2 cycles) using the previous value hold method. 2) Module isolation: If the data cannot be restored for 5 consecutive minutes, the system will determine that the sensor is faulty. At this time, the adaptive correction of the PCE coefficient related to the sensor data will be frozen (i.e., its corresponding weight w j 3) Alarm and Mode Switching: The system generates an XX sensor fault alarm. If the faulty sensor is a core safety-related sensor, the entire adaptive correction module will be suspended and switched to non-adaptive control mode.

[0168] Scenario 3: The system status is predicted to cross the safety boundary.

[0169] In the prediction phase of step S620, it is found that even under the optimal control strategy, a certain point in the future state prediction trajectory (such as the predicted exhaust temperature) will still reach or fall below the set safety boundary (such as the actual gas freezing point temperature + safety margin).

[0170] Through an independent, high-priority Safety Veto mechanism, specifically:

[0171] 1) Overwrite optimization results: The system will immediately discard the economic optimal control strategy calculated in step S630.

[0172] 2) Execute the emergency plan: Instead, execute a pre-programmed emergency control instruction with safety as the sole goal, for example, immediately open the dehumidification system to maximum power to force the humidity to the lowest level.

[0173] 3) Lock and Alarm: In this emergency mode, the system will maintain maximum safety measures until it detects that the system status has been stable within the safety zone for multiple consecutive periods (e.g., 10 minutes). Only then will it unlock and resume normal optimized control. At the same time, it will issue the highest level of safety alarm to the operator.

[0174] Example 10. Performance comparison between the present application and the prior art.

[0175] Define a traditional humidity control method as a baseline, for example, an ON / OFF control strategy with a fixed dew point temperature (e.g., -45°C, calculated offline based on the worst-case operating conditions).

[0176] A dynamic simulation platform for a million-kilowatt compressed air energy storage system built on MATLAB / Simulink and calibrated with actual data.

[0177] Design a 72-hour variable load operating curve and simulate slow equipment degradation during this period (for example, starting at the 24th hour, heat exchanger efficiency decreases linearly at a rate of 0.05% per hour). Select a snapshot at the 48th hour of the simulation. Inputs at that moment include the grid dispatch instructions, gas storage pressure, temperature, and other specific input values.

[0178] The cumulative performance indicators of the method of the present invention and the traditional baseline method over the entire 72-hour simulation period are as follows.

[0179] Performance indicators Traditional baseline method Method of the present invention Performance improvements Dehumidification system cumulative power consumption (kWh) 12,580 10,693 Save 15.0% Minimum safety margin between exhaust gas temperature and freezing point temperature (°C) 2.8°C 5.5°C 96.4% increase Number of times the safety boundary is touched 0 0 - Humidity control mean absolute error (°C dew point) 1.8°C 0.4°C Reduced by 77.8% Model prediction RMSE (after degradation occurs) (°C) 2.1°C 0.6°C Higher precision

[0180] Example 11: This example, based on Example 7, provides an enhanced method for detecting and self-repairing spatiotemporal coupling instability in PCE proxy model coefficients. This method is performed after the original constrained optimization solution step (S730) but before the covariance matrix update to address the issue of high-frequency oscillations in the time dimension of the updated PCE coefficients, which can lead to physical distortion of the model.

[0181] In other words, before reconstructing the updated PCE proxy model using the updated PCE coefficient vector, the following steps are also included:

[0182] Collect the updated PCE coefficient vector and generate PCE coefficient time series data by combining it with the PCE coefficient history record;

[0183] Perform time-frequency domain analysis on the PCE coefficient time series data to generate a two-dimensional time-frequency information containing phase information for each PCE coefficient;

[0184] Identifying physical coupling coefficient pairs based on a preset coupling constraint matrix, and calculating a physical coherence degree used to characterize the phase coordination between the physical coupling coefficient pairs based on the two-dimensional time-frequency information;

[0185] comparing the physical coherence with a physical coherence threshold to generate an instability detection result;

[0186] In response to the instability detection result indicating the presence of instability, generating a phase repair signal and converting the phase repair signal into a coefficient correction amount;

[0187] The coefficient correction is integrated into a constrained optimization problem, and the constrained optimization problem is solved to obtain the repaired PCE coefficient vector.

[0188] Among them, a two-dimensional time-frequency information containing phase information is generated for each PCE coefficient, including:

[0189] Short-time Fourier transform technology is used to process the PCE coefficient time series data to obtain time-frequency two-dimensional information, which includes the frequency component of each PCE coefficient and its phase trajectory evolving over time.

[0190] The calculation of a physical coherence degree includes:

[0191] Extracting the actual phase difference of the physical coupling coefficient within a preset frequency range from the time-frequency two-dimensional information;

[0192] The actual phase difference is compared with the expected phase difference derived based on the coupling constraint matrix to obtain the phase deviation;

[0193] The phase deviations of all physical coupling coefficient pairs are accumulated to form the physical coherence.

[0194] The generating of the phase repair signal includes: constructing a virtual master oscillator having a preset frequency, and setting a desired phase offset relative to the virtual master oscillator for each PCE coefficient;

[0195] Calculating a phase error between a current phase of the unstable PCE coefficient and a desired phase offset;

[0196] A phase error feedback controller is used to generate a phase repair signal based on the accumulated value of the phase error.

[0197] The coefficient correction is integrated into a constrained optimization problem, and the constrained optimization problem is solved, including:

[0198] The coefficient correction is added as a penalty term to the objective function of the constrained optimization problem;

[0199] While keeping the physical constraints defined by the coupling constraint matrix unchanged, the modified constrained optimization problem is solved to produce the repaired PCE coefficient vector.

[0200] After obtaining the repaired PCE coefficient vector, the method further includes: recalculating a verification physical coherence based on the repaired PCE coefficient vector;

[0201] The verified physical coherence is compared with the physical coherence threshold. When the verified physical coherence is lower than the physical coherence threshold, the repaired PCE coefficient vector is confirmed to be valid.

[0202] According to one aspect of the present application, the method specifically comprises the following steps:

[0203] Step S1110, preprocessing of coefficient time series data.

[0204] PCE coefficient time series data is a two-dimensional data matrix that records the PCE coefficient vector for each minute within a preset time window (for example, 120 minutes). Long-term degradation trends primarily refer to slow, unidirectional coefficient changes over days or months, caused by factors such as equipment wear and material aging.

[0205] After each PCE coefficient update, the system obtains the latest coefficient vector and combines it with the past 119 coefficient vectors stored in the historical database to form a sliding window of data covering 120 time points. This data is first normalized using the Z-score to eliminate the impact of magnitude differences between coefficients. Subsequently, a low-pass filter (such as a sliding average or Butterworth low-pass filter) is used to extract and remove long-term degradation trends in the data, retaining only the mid- and high-frequency components that reflect the system's dynamic characteristics. This forms the "residual coefficient time series data" for subsequent analysis.

[0206] Step S1120: time-frequency domain analysis based on short-time Fourier transform.

[0207] Short-time Fourier transform (STFT) analyzes the frequency components and phase information of a signal in different time periods by adding windows.

[0208] The time-frequency two-dimensional information contains a complex data structure with four dimensions: time, frequency, amplitude, and phase, which is used to describe how the frequency component of each PCE coefficient evolves over time.

[0209] The STFT technique is used to analyze the preprocessed residual coefficient time series data. Specifically, a 30-minute Hanning window is set as the analysis window and slid forward in steps of 7.5 minutes. For each PCE coefficient time series, its spectral characteristics within each analysis window are calculated, generating two-dimensional time-frequency information containing both amplitude and phase. This information reveals the specific frequencies where energy is concentrated for each coefficient and how the phase of these frequency components changes over time.

[0210] Step S1130: Physical constraint-based coherence measurement. A physical coupling coefficient pair refers to two or more PCE coefficients that co-occur in the same physical constraint equation (e.g., the energy conservation equation). Physical coherence is a dimensionless metric used to quantify the degree of phase coordination between physical coupling coefficient pairs. Lower values ​​indicate better phase synchronization.

[0211] First, based on the PCE coefficient coupling constraint matrix constructed in Example 5, all coefficient pairs with physical coupling relationships are identified. Subsequently, within the dominant frequency range of typical operating condition changes of the energy storage system (for example, 0.1 to 0.5 cycles per hour), the phase time series of each coupling coefficient pair within this frequency band is extracted. The "actual phase difference" between the two at each moment is calculated and compared with the "expected phase difference" derived based on the physical constraint relationship (for example, it should be 0 or a fixed constant in stable coupling). The phase difference deviation degree of all coupling coefficient pairs is weighted and summed to obtain the final "physical coherence" indicator.

[0212] Step S1140 : instability detection, diagnosis and virtual master oscillator design.

[0213] Based on a large amount of historical stable operation data, the distribution of physical coherence is statistically analyzed, and warning thresholds (such as 1.5 times the standard deviation of the mean) and critical thresholds (such as 3 times the standard deviation of the mean) are set. When the real-time calculated physical coherence exceeds the critical threshold, the system determines that "space-time coupling instability" has occurred. Once instability is detected, the system will further locate which specific coupling coefficient pairs have phase disorder. At the same time, in order to repair it, the system designs a "virtual master oscillator" whose frequency is set to the dominant frequency of the energy storage system's charge and discharge cycle (for example, 0.25 cycles / hour), and assigns each PCE coefficient a desired phase offset relative to the master oscillator as a benchmark for phase synchronization.

[0214] Step S1150: Repair and re-optimization based on phase error feedback.

[0215] For the coefficients diagnosed as unstable, the error between their current phase and the expected phase is calculated.

[0216] Based on the integral control principle, a phase error feedback controller is designed, and the intensity of the "phase repair signal" it outputs is proportional to the accumulated value of the phase error.

[0217] The abstract phase repair signal is converted into a specific numerical correction value based on the physical meaning and typical change amplitude of each coefficient, forming a "coefficient correction value vector".

[0218] The coefficient correction vector is used as a new penalty term and added to the objective function of the original constrained optimization problem. The new objective function becomes: minimize { || x new - (x old + K modified ·z)|| 2 + w·||x new - (x old + correction vector)|| 2}, where w is the penalty weight. While keeping all the original physical constraints unchanged, the optimization problem is re-solved to obtain the “repaired PCE coefficient vector”.

[0219] Step S1160: Verify the repair effect.

[0220] The repaired PCE coefficient vector temporarily replaces the last data point in the historical sequence, and the physical coherence is quickly recalculated. If the calculated new coherence is significantly lower than the warning threshold, the repair is confirmed to be successful, and the repaired coefficient vector is officially adopted and passed to subsequent steps. Otherwise, the system abandons the repair result and uses the original optimization result that only satisfies the static constraints. A repair failure event is recorded to ensure the safety of the repair operation itself.

[0221] Example 12: Describes how to establish a clear and definite physical correspondence for the PCE coefficient.

[0222] In this embodiment, it is assumed that the simplified PCE proxy model is used to predict the turbine outlet temperature T out , its expression is: T out =c0Ψ0+c1Ψ1(P in )+c2Ψ2(T in )+c3Ψ1(P in )Ψ1(T in )+c4Ψ2(η turb ) where P in is the inlet pressure, T in is the inlet temperature, η turbis the turbine isoentropic efficiency. Using the sensitivity analysis method described in Example 5, the following table of physical meanings, dimensions, and typical value ranges was established for each coefficient: c0 is the base outlet temperature (average value under baseline operating conditions); c1 is the linear influence factor of inlet pressure on outlet temperature; c2 is the linear influence factor of inlet temperature on outlet temperature; c3 is the cross-coupling effect factor of pressure and temperature; and c4 is the sensitivity factor to the impact of turbine efficiency degradation on outlet temperature.

[0223] For example, a negative value for c1 indicates that higher inlet pressure leads to lower outlet temperature after expansion, which is consistent with the laws of thermodynamics. c4 is directly related to the health of the equipment, and changes in its value can be directly interpreted as a decrease in turbine efficiency, making the model's adaptive correction results highly interpretable in engineering.

[0224] Example 13: Supplementing step S540 in Example 5, describing the constraint matrix construction process.

[0225] In this embodiment, the Sobol global sensitivity analysis method is used to map PCE coefficients to physical quantities. The Sobol method can quantify the contribution of a single input (here, the PCE coefficient) and the interaction between inputs to the variance of the output (physical quantity).

[0226] For the trained initial PCE model, all its coefficients c i As the input variable of the Sobol method, the predicted physical quantity (such as turbine efficiency η) is used as the output. Within the reasonable range of each coefficient (see Example 12), tens of thousands of input coefficient sample combinations are generated.

[0227] Substitute each sample combination into the PCE model and calculate the corresponding physical quantity output value.

[0228] Using the Saltelli sampling strategy, calculate each coefficient c i The first-order sensitivity index S for the output η i and the global sensitivity index S Ti . S i Indicates the contribution of this coefficient alone to the output variance, S Ti It represents the contribution of this coefficient alone and in interaction with other coefficients to the variance of the output.

[0229] Set the sensitivity index threshold, for example, 0.1. For a certain physical quantity η, if a certain coefficient c j The global sensitivity index S Tj >0.1, and its first-order sensitivity index S j The coefficients of the physical quantity are ranked in the top (for example, the top three), and c is established. j ->η strong mapping relationship.

[0230] Through this method, the key PCE coefficients that dominate various physical quantities can be quantitatively and objectively identified, providing a reliable foundation for the subsequent construction of accurate constraint matrices based on physical laws.

[0231] In short, the physical coupling constraint correction mechanism ensures the physical authenticity and predictive reliability of the data-driven adaptive agent model throughout its entire life cycle, solving the fatal flaw of traditional data-driven models in closed-loop control, which may cause physical errors due to data fitting. Specifically, in the offline stage, a mapping relationship is established between the abstract PCE coefficient and specific physical quantities (such as efficiency and heat loss) through sensitivity analysis, and then universal physical laws such as energy conservation are converted from macroscopic partial differential equations into a specific linear constraint matrix A composed of PCE coefficients. eq Then, in the online correction phase, the constraint matrix is ​​used as a non-violated hard constraint of the constrained optimization problem. This means that no matter what form the input sparse measurement data z takes, the updated coefficient vector x given by the solver (such as the SQP algorithm) new , its solution space is constrained from the outset to a subspace satisfying the laws of physics. Therefore, in the context of compressed air energy storage power plants, where safety and redundancy requirements are extremely high, this method ensures that even when sensor data deviates or equipment status fluctuates dramatically, the modified model will never output absurd predictions that violate the first law of thermodynamics (such as efficiency greater than 1). This provides an absolutely reliable prediction trajectory for upper-level optimization control, greatly improving the robustness and safety of the entire humidity control system.

[0232] The variable weight adaptive correction method gives the model correction process intelligent diagnostic capabilities, enabling it to distinguish error sources and make differentiated and efficient corrections, thereby achieving both control stability when the system is healthy and rapid adaptability when the equipment is degraded. This effect is achieved through the following logical chain: by calculating the rate of change Δσ of the variance of the prediction error vector z within the sliding time window 2 , transforming the abstract error sequence into a quantifiable indicator that characterizes the error trend. 2 This means that the error is mainly due to random measurement noise, and the continuously growing Δσ 2 This points to a systematic model mismatch caused by equipment performance degradation. This diagnostic result is consistent with the degradation sensitivity index S of each coefficient. sensitivity They are substituted into the nonlinear mapping function to calculate a unique update weight w for each coefficient j In the actual operation of the compressed air energy storage power station, it means that when the system operates smoothly, the weights of all coefficients w jThe model is insensitive to tiny sensor noise, which avoids frequent jitter of control instructions and saves unnecessary dehumidification energy. Once the turbine efficiency begins to decline due to wear, the system will immediately detect the corresponding physical quantity error Δσ 2 It continues to increase and quickly assigns an extremely high weight to the PCE coefficient that controls turbine efficiency, allowing it to be quickly and significantly corrected to match reality, thereby ensuring the power generation efficiency and economy of the power plant throughout its life cycle.

[0233] The smooth transition switching algorithm of the state equation ensures the smoothness and continuity of the safety constraint boundary in the rolling time domain optimization, thereby ensuring the numerical stability of the optimization solver and the smoothness of the final control instruction output, and avoiding oscillation of the control system due to model switching. In the million-kilowatt compressed air energy storage system, the rolling time domain optimization solver (such as SQP) is highly sensitive to the characteristics of the constraint boundary. If the traditional hard switching method is used, when the system pressure P crosses the switching point of the state equation (such as 8MPa), the freezing point temperature T dewpoint The state equation will change from RK equation to PR equation instantly, resulting in T as a safety constraint. dewpoint The value produces a step-like jump. This discontinuous constraint boundary will seriously interfere with the convergence process of the optimization solver, which may lead to solution failure or the output of a sharply changing and non-smooth control sequence. This method performs a weighted average of the calculation results of the two equations in the boundary area near the switching point, and the weight w(P) is a continuous function of the pressure P, thereby gluing the segmented constraint boundaries into a smooth curve. This allows the optimization solver to find the optimal solution stably and reliably, and the humidity control command finally output is also continuous and smoothly changing, which is crucial for driving real industrial systems with huge physical inertia. It effectively avoids the impact on actuators (such as valves and dehumidifiers) and improves the overall operation quality of the system.

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

Claims

1. A method for controlling humidity in a large-scale compressed air energy storage power station gas storage reservoir, characterized in that: include: Obtain grid dispatch instructions and real-time power plant operating data; use a pre-configured PCE agent model to calculate the instantaneous optimal humidity set point and generate a model prediction output trajectory; Acquire sparse real measurement data, combine it with the model prediction output trajectory, and use the pre-configured PCE coefficient coupling constraint matrix to correct the coefficients of the PCE proxy model to generate an updated PCE proxy model for the next cycle; Convert the instantaneous optimal humidity set point into a dehumidification equipment control instruction and issue it for execution, and collect the execution results as part of the power plant real-time operating condition data for the next control cycle; Modify the coefficients of the PCE proxy model, including: Determine the prediction error vector based on the difference between the model's predicted output trajectory and the sparse true measurement data; The constrained optimization problem is constructed with the PCE coefficient coupling constraint matrix as the physical constraint condition and the prediction error vector as the optimization basis; Solve the constrained optimization problem to obtain the updated PCE coefficient vector that satisfies the physical coupling constraints; Using the updated PCE coefficient vector, reconstruct and generate an updated PCE proxy model; The construction process of the preconfigured PCE coefficient coupling constraint matrix includes: A sensitivity analysis is performed on the coefficients and basis functions of the initial PCE proxy model to identify the physical quantities dominated by each PCE coefficient and establish a coefficient-physical quantity mapping relationship. Based on the coefficient-physical quantity mapping relationship and the system's thermodynamic conservation law, the mathematical constraints that must be satisfied between PCE coefficients are derived; The mathematical constraint relationship is formally constructed as a PCE coefficient coupling constraint matrix; The mathematical constraints include: equality constraints derived from the conservation laws of thermodynamics, and inequality constraints derived from the physical boundary conditions of the system; Calculate the instantaneous optimal humidity set point using a rolling horizon optimization process, including: At each prediction step of the rolling horizon optimization, the optimal real gas state equation is adaptively selected from the preset pressure-state equation mapping table based on the current pressure value; Using the selected real gas equation of state, calculate the corresponding real gas freezing point temperature; Taking the real gas freezing point temperature as a safety constraint, the receding horizon optimization is solved to obtain the instantaneous optimal humidity set point.

2. The method according to claim 1, characterized in that Solve constrained optimization problems, including: For constrained optimization problems, Lagrange multipliers are introduced to construct Lagrange functions; Based on the Lagrangian function, a set of KKT conditional equations is established to characterize the necessary conditions for the optimal solution; Iteratively solve the KKT conditional equations and calculate the updated PCE coefficient vector.

3. The method according to claim 1, characterized in that Correcting the coefficients of the PCE proxy model also includes updating the filter state covariance matrix, specifically: Based on the PCE coefficient coupling constraint matrix, a projection matrix is ​​constructed that can project uncertainty into the physical constraint subspace; The projection matrix is ​​used to modify the standard covariance update process to generate an updated state covariance matrix that reflects the effects of physical coupling constraints.

4. The method according to claim 1, wherein Before constructing the constrained optimization problem, also include: Identify the current degradation stage of the device based on the statistical characteristics of the prediction error vector within the sliding time window, including the rate of change of the variance; In response to the identified degradation stage, a corresponding adaptive update weight vector is generated for each PCE coefficient and applied to the formulation of the constrained optimization problem.

5. The method according to claim 4, characterized in that Generate the corresponding adaptive update weight vector for each PCE coefficient, including: Based on the identified degradation stage, the degradation sensitivity index of each PCE coefficient is calculated; Using the preset nonlinear mapping function, the weight value of each PCE coefficient is calculated according to the degradation sensitivity index and the rate of change of the variance; The weight values ​​are normalized to form an adaptively updated weight vector.

6. The method according to claim 4, characterized in that Identify the current stage of degradation of the equipment, including: Calculate the variance of the prediction error vector within the sliding time window and determine the rate of change of the variance between consecutive time windows to obtain the error variance change rate; The error variance change rate is compared with at least one preset threshold, and the degradation stage is classified into one of a plurality of predefined degradation modes according to the comparison result.

7. The method according to claim 6, characterized in that If the current pressure value is in the boundary area of ​​the adjacent pressure interval, the corresponding real gas freezing point temperature is calculated using the selected real gas state equation, which is achieved through a smooth transition algorithm, including: Two independent freezing point temperature values ​​are calculated using two real gas state equations corresponding to adjacent pressure intervals; Determine a weight coefficient based on the relative position of the current pressure value within the boundary area; The weight coefficient is used to perform a weighted average on the two independent freezing point temperature values ​​to obtain the final true gas freezing point temperature after a smooth transition.

Images