Multi-mode constant pressure gas storage experiment equivalent simulation method

Abstract

The application discloses a kind of multi-mode constant pressure gas storage experimental equivalent simulation methods, containing leakage characteristic evaluation process, specifically: the closed loop sensitivity of system is spectrum shaping, and constructs sensitivity notch in target perturbation frequency band;Generate the perturbation signal satisfying mass-enthalpy-time first moment triple constraint, and after its anti-resonance spectrum pre-shaping, inject in notch frequency band, to obtain the system response data of high signal-to-noise ratio;Online estimation of structured leakage is carried out using mode-aware instrumental variable method, and the continuity of cross-mode estimation is guaranteed through bridging consistency mechanism.The application can obtain unbiased, continuous leakage estimate, and provides effective tool for accurate evaluation of constant pressure energy storage system.

Description

TECHNICAL FIELD

[0001] The present application relates to compressed air energy storage, in particular, a multi-mode constant pressure gas storage experiment equivalent simulation method. BACKGROUND

[0002] Under the background of building new power systems, large-scale long-time energy storage technology is a key support to improve the flexibility of the power grid and accommodate a high proportion of renewable energy. Among them, compressed air energy storage (CAES) has become one of the most promising technology routes due to its large capacity, long service life, and strong geographical adaptability. The traditional constant-volume compressed air energy storage system continuously decreases the pressure of the gas storage chamber during power generation, which causes the turbine expander to operate under variable conditions (sliding pressure) for a long time, deviating from the design point and affecting the efficiency. To address this challenge, constant pressure energy storage technology has emerged, which maintains constant pressure during the main energy release period through hydraulic compensation or flexible airbags, which can improve system operation efficiency and stability. However, these advanced constant pressure energy storage schemes are mostly in the theoretical design and simulation stage, and the construction of large-scale physical entities is costly and time-consuming. Therefore, in the laboratory environment, by building a high-fidelity physical scale experiment platform, equivalent simulation and experimental verification of the dynamic characteristics, control strategy and potential risks of the constant pressure energy storage system is the only way to accelerate technology iteration, reduce engineering risks and promote industrialization.

[0003] Currently, experimental research on complex fluid and thermal systems mainly uses conventional system identification and state monitoring methods. In terms of system dynamic characteristics testing, researchers usually use step signals, pulse signals or single-frequency sinusoidal signals as excitations to obtain the transfer function or frequency response characteristics of the system. These methods are well established in open-loop or linear systems and can effectively identify the main modalities and response gains of the system. In terms of equipment health state monitoring, especially leakage detection, traditional offline methods dominate, such as the pressure drop method, which estimates the total leakage rate by monitoring the rate of pressure drop over time in a sealed container after the system is shut down. For online monitoring, some methods infer the existence of leakage based on the measurement deviation of input and output flow, pressure, temperature and other steady-state parameters by establishing macroscopic energy or mass balance equations. In terms of signal processing, time or frequency domain filtering techniques are commonly used to improve signal-to-noise ratio, and basic least squares methods are used to fit the simplified leakage model parameters.

[0004] However, directly applying the above conventional experimental and monitoring techniques to the near-constant pressure energy storage experimental system based on closed-loop feedback control in this paper will face deep technical problems caused by its inherent physical characteristics and interrelated. These difficulties mainly manifest in three aspects: effective excitation of weak signals under the background of closed-loop strong inhibition, biased parameter estimation caused by endogenous system, and the lack of estimation continuity under multi-mode dynamic switching. SUMMARY

[0005] To achieve the object, the present application provides a multi-mode constant pressure gas storage experiment equivalent simulation method, which aims to solve the technical problem that under closed-loop constant pressure control, weak leakage signals are difficult to identify and parameter estimation is biased due to the strong inhibition characteristics of the controller and the endogenous nature of the system.

[0006] Technical scheme, according to one aspect of the present application, a multi-mode constant pressure gas storage experiment equivalent simulation method, comprising a leakage characteristic evaluation process, which is specifically:

[0007] Determine the target perturbation frequency band; read the perturbation signal and perform spectral pre-shaping to make the main energy fall into the target perturbation frequency band, inject the perturbation signal after spectral pre-shaping into the system to obtain system response data; based on the system response data, estimate the structured leakage related to the system pressure online to obtain the structured leakage estimation value.

[0008] As a possible implementation manner of one aspect of the present application, generating the perturbation signal includes: generating an initial perturbation signal, applying a pre-configured online small signal model to map it to a mass flow rate increment sequence; obtaining a real-time temperature sequence to calculate a corresponding enthalpy flow rate increment sequence; within a preset time window, simultaneously applying zero integral constraints to the mass flow rate increment sequence and the enthalpy flow rate increment sequence to determine the perturbation signal that satisfies the mass constraint and the enthalpy constraint and store it.

[0009] As a possible implementation manner of one aspect of the present application, within a preset time window, simultaneously applying zero integral constraints to the mass flow rate increment sequence and the enthalpy flow rate increment sequence includes: constructing a time-weighted mass flow rate increment sequence; within a preset time window, applying zero integral constraints to the time-weighted mass flow rate increment sequence to make the perturbation signal satisfy the mass constraint, the enthalpy constraint and the time first moment constraint.

[0010] As a possible implementation manner of one aspect of the present application, making the perturbation signal satisfy the mass constraint, the enthalpy constraint and the time first moment constraint includes: generating a corresponding initial mass flow rate increment sequence based on a pre-configured signal template; using a constrained optimization solving method, under the condition of satisfying the triple constraints, solving to obtain a constrained mass flow rate increment sequence; applying a valve-flow inverse mapping model to inversely solve the constrained mass flow rate increment sequence to the perturbation signal.

[0011] As a possible implementation manner of one aspect of the present application, the spectral pre-shaping of the perturbation signal includes: obtaining the fluid main mode frequency of the system; using a filter to process the perturbation signal to suppress the energy components of the perturbation signal at the main mode frequency and its harmonic frequencies.

[0012] As a possible implementation manner of one aspect of the present application, the online estimation of the structured leakage related to the system pressure comprises: identifying a current operation mode of the system based on system response data, the operation mode comprising a charging mode, a discharging mode and a holding mode; and dividing the system response data into mode data windows corresponding to each operation mode according to the switching of the operation mode, for subsequent structured leakage estimation.

[0013] As a possible implementation manner of one aspect of the present application, for subsequent structured leakage estimation, the method comprises: starting a disturbance-free window to suspend the injection of the perturbation signal at the mode switching between the mode data windows; defining a bridge window by using the data around the mode switching in adjacent mode data windows; and correcting the leakage model parameters obtained from the previous mode data window in the bridge window to realize the continuity of the cross-mode estimation.

[0014] The beneficial effects are that the above technical solutions can obtain unbiased and continuous leakage estimation values, solve the technical problems that the weak leakage signal is difficult to identify and the parameter estimation is biased under closed-loop constant pressure control, and provide an effective tool for accurate evaluation of the constant pressure energy storage system. BRIEF DESCRIPTION OF DRAWINGS

[0015] Figure 1 is a general flowchart of a multi-mode constant pressure gas storage experiment equivalent simulation method.

[0016] Figure 2 is a flowchart of generating a perturbation signal.

[0017] Figure 3 is a step flowchart of simultaneously applying zero integral constraints to the mass flow rate increment sequence and the enthalpy flow rate increment sequence within a preset time window.

[0018] Figure 4 is a step flowchart of making the perturbation signal satisfy the mass constraint, the enthalpy constraint and the time first moment constraint.

[0019] Figure 5 is a step flowchart of performing spectral pre-shaping on the perturbation signal. DETAILED DESCRIPTION

[0020] In order to solve the problems existing in the prior art, the applicant has conducted in-depth research and found that high-fidelity system dynamic identification requires that the test signal have sufficient excitation persistence and spectral richness, but the traditional step or large disturbance injection method will impact the stable working condition of the system, and may even induce safety risks, which makes online fine identification impractical; and offline testing cannot capture the state-dependent dynamics of the system under real operating conditions.

[0021] The lack of identification ability directly leads to the difficulty in quantifying the small, pressure-dependent leakage. In a closed-loop constant pressure control system, the compensatory effect of the controller actively masks the slow pressure drop caused by the small leakage, making the leakage characteristics submerged in the background noise and control action. Without a persistent probe signal that is statistically decoupled from the control action, the leakage characteristic parameters cannot be accurately separated from the closed-loop data, resulting in a serious deviation in the estimation results.

[0022] Existing diagnostic methods are generally based on static or time-invariant system models, ignoring the fact that the compressed air energy storage system is a complex thermodynamic system whose dynamic characteristics (such as the liquid column resonance frequency in a hydraulic drive system) will drift with changes in operating conditions (such as temperature and pressure). This model mismatch problem is particularly pronounced during the strong transient process of charging and discharging mode switching, further weakening the long-term effectiveness and reliability of the diagnostic results.

[0023] The main task of the system is to maintain constant pressure, and the pressure closed-loop controller is designed to actively suppress the effects of internal and external disturbances on system pressure. This leads to an inherent contradiction: in order to detect the system dynamic changes caused by small leakage, a small disturbance excitation signal needs to be injected into the system; but the injected small disturbance signal will be identified by the controller as a disturbance and will be greatly weakened and suppressed. This makes the system response signal containing leakage information extremely weak, submerged in strong background noise and controller activity, with extremely low signal-to-noise ratio, and traditional excitation and identification methods directly fail in this scenario, which is the excitation signal being controlled annihilation problem.

[0024] Even if a weak response signal can be obtained, the accuracy of parameter estimation is also facing severe challenges. In a closed-loop system, the output of the controller (i.e. part of the excitation signal) is a function of the system state (pressure), while the structured leakage itself is also a function of the system state (pressure). This inherent correlation between input and state is referred to as endogeneity in statistics. If the traditional least squares method is used for parameter identification without considering this characteristic, the basic assumption of the method (i.e. the explanatory variable and the residual term are not related) is seriously violated, and the estimated leakage parameters will inevitably have systematic bias and lack consistency, and cannot reflect the true leakage level. This is the closed-loop system estimation bias problem, which is more fundamental than the simple signal-to-noise ratio problem.

[0025] Constant pressure energy storage system needs to switch between multiple modes such as charging, discharging and pressure maintaining. Under different modes, the fluid flow direction, dynamic characteristics and even leakage mechanism of the system may change. The existing identification methods are usually based on a single steady state condition, lacking the ability to perceive the running mode. If the data across modes are mixed, model mismatch error will be introduced; if segmented processing is used, due to the lack of constraints across mode boundaries, the estimated leakage value will produce a sharp jump at the mode switching point, which is not physically realistic, and destroys the continuity and credibility of state estimation. This is the cross-mode estimation discontinuity problem. These three problems are interrelated and constitute the technical bottleneck for online and accurate leakage diagnosis in advanced energy storage experimental systems.

[0026] In order for those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be described clearly and completely below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only some of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of the present application.

[0027] The terms "first", "second", etc. in the specification and the above-mentioned drawings are used to distinguish different objects, not to describe a specific order. In addition, the terms "include" and "have" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, device, or product that includes a series of steps or units is not limited to the listed steps or units, but can optionally include steps or units not listed, or can optionally include other steps or units inherent to the process, method, product or end.

[0028] In this document, the term "embodiment" means that the specific features, structures or characteristics described in connection with the embodiment can be included in at least one embodiment of the present application. The appearance of this phrase in various places in the specification does not necessarily refer to the same embodiment, nor is it mutually exclusive or alternative to other embodiments. It is explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0029] Embodiment one, an equivalent simulation method is proposed, which can be applied to a multi-mode constant pressure gas storage experimental device. The device aims to simulate the high pressure constant pressure gas storage process in deep water environment through the method of gas pressure supplementing water head in laboratory environment.

[0030] Specifically, the experimental device can include water-gas coexistence tanks and water-gas coexistence tanks, and valves and pipeline systems connected therebetween. Including motor M, compressor, heat exchanger, turbine, generator G, high-temperature heat storage tank, low-temperature heat storage tank. In the system preset stage, the water in the water storage pool is pressed into the two water-gas coexistence tanks by the water pump / water turbine, the water-gas coexistence tank can be used to simulate the upper reservoir, and the water-gas coexistence tank is used to simulate the lower reservoir, to realize the constant pressure gas storage in the water-driven mode. Flexible air bags can also be arranged underwater in the water-gas coexistence tank, and the deep water head is simulated by high-pressure air above the water surface in the tank to realize the constant pressure gas storage simulation of the underwater flexible air bag.

[0031] In this embodiment, the principle of air pressure supplementing water head is to use the principle of communicating vessels to replace the huge actual water depth by adjusting the pressure of the gas phase space at the upper part of the two water-gas coexistence tanks. For example, the pressure relationship between the two tanks can be represented as: P1+ρgh1=P2+ρgh2; wherein: P1 is the air pressure at the upper part of the first water-gas coexistence tank; h1 is the liquid level height of the first water-gas coexistence tank; P2 is the air pressure at the upper part of the second water-gas coexistence tank; h2 is the liquid level height of the second water-gas coexistence tank; ρ is the density of water; g is the acceleration of gravity. By actively controlling P1 and P2, a smaller liquid level difference Δh=h1-h2 can be used to simulate the equivalent water head generated by a large pressure difference (P2-P1).

[0032] The device can flexibly switch between multiple operating modes by controlling the opening and closing combination of different valves. For example, in the simulation of water-driven constant pressure gas storage mode, part of the valves are always kept closed. When charging, open the air inlet valve and the water outlet valve, store compressed air in the gas storage tank (such as the second water-gas coexistence tank), and at the same time press the water from the gas storage tank into the upper reservoir (such as the first water-gas coexistence tank); when generating electricity, open the water inlet valve and the air outlet valve to make the water in the upper reservoir flow into the gas storage tank, and the high-pressure air is squeezed to the expander to do work. In the simulation of constant volume compressed air energy storage working condition, another part of the valves are closed throughout the process, only the air inlet valve is opened when charging, and the air outlet valve is switched when generating electricity. In the simulation of underwater constant pressure compressed air energy storage mode, the upper reservoir needs to be isolated, the underwater air inlet valve is opened when charging, and the underwater air outlet valve is opened when generating electricity.

[0033] In such a multi-mode, high-pressure complex fluid system, there is a structured leakage related to the system pressure and operating mode, such as valve leakage, loose flange sealing, etc. Such leakage can seriously affect the evaluation of energy storage efficiency and the judgment of system safety. Therefore, how to accurately and online identify the tiny structured leakage in the context of near-constant pressure closed-loop control is a key technical problem for the experimental device to provide effective experimental data, and also constitutes the realistic background of the technical problems to be solved by the present invention.

[0034] Embodiment two, the embodiment provides a general flow of a multi-mode constant pressure gas storage experimental equivalent simulation method, as shown in Figure 1As shown, the flow solves the problem that the leakage is difficult to accurately estimate due to the mutual coupling of pressure fluctuation and micro-leakage signal under the condition of near constant pressure closed loop control.

[0035] The method comprises the following steps:

[0036] Step one, determine the target perturbation frequency band,

[0037] In this embodiment, the closed loop sensitivity of the system refers to the closed loop transfer characteristic from the mass flow input m* of the system to the pressure output P. Since the system needs to set a closed loop pressure compensation controller to suppress the influence of external disturbance on the system pressure in order to maintain constant pressure. The above suppression effect is manifested in the frequency domain as the closed loop sensitivity of the system is very low in a certain frequency band. The energy of the directly injected perturbation signal will be greatly weakened by the controller, resulting in extremely low signal-to-noise ratio, which cannot be used for subsequent parameter estimation.

[0038] Therefore, this step selects a frequency band far from the main noise frequency and resonance frequency of the system as the target perturbation frequency band Ω d . Further, by superimposing a notch filter on the pressure compensation loop or introducing a feedforward / lag element in the controller design environment, etc., the closed loop sensitivity function is actively spectrum shaped to selectively and locally raise the sensitivity in the target perturbation frequency band Ω d , form a sensitivity notch, so that the perturbation signal injected into this frequency band can be reflected in the system pressure response with a higher signal-to-noise ratio, creating conditions for subsequent accurate estimation.

[0039] In another embodiment of the present application, this step can also include spectrum shaping the closed loop sensitivity of the system to construct a sensitivity notch in the target perturbation frequency band.

[0040] Step two, read the pre-generated perturbation signal, and spectrum pre-shape the perturbation signal so that the main energy of the perturbation signal falls within the target perturbation frequency band determined by the foregoing step. Inject the spectrum pre-shaped perturbation signal into the system to obtain system response data.

[0041] Specifically, when generating the perturbation signal, in order to minimize its disturbance to the overall operating state of the system, a signal that satisfies a certain integral constraint is preferably generated, such as the mass-enthalpy-time first moment triple constraint signal which will be detailed in subsequent embodiments. After generating the signal, in order to make its energy be efficiently utilized and not excite potential fluid resonance of the system, it needs to be spectrum pre-shaped. In some examples, it can be realized by a band-pass filter, and the passband is strictly aligned with the sensitivity notch frequency band Ω dAfter the spectrum pre-shaping, the main energy of the perturbation signal is concentrated in the frequency band where the system response is most sensitive. The signal is then loaded to the actuator of the system (e.g. the valve position command of the control valve) and injected into the running system. The response data of the system are collected synchronously at high frequency, mainly including time series data such as system pressure P[k], mass flow rate m*[k], temperature T[k] and valve position u[k] (where k is the time sampling index).

[0042] Step three, based on the system response data, estimate the structured leakage related to the system pressure online to obtain the structured leakage estimation value.

[0043] Specifically, after obtaining the injection-response data pair, the structured leakage estimation algorithm based on pattern perception is used for leakage estimation. According to the information such as the sign of the mass flow rate m*[k], it is identified that the system is currently in what kind of operation mode such as charging, discharging or pressure maintaining, and the data is divided into corresponding windows. The leakage is modeled as a low-order parameterized term related to the pressure, for example, a function vector containing a constant bias term, a term linearly related to the pressure difference and a term related to the square root of the pressure difference. Since in the closed-loop system, there is a correlation between the perturbation signal as input and the pressure as state, direct least squares estimation will produce biased results. In order to solve this problem, the instrumental variable (IV) or two-stage least squares (2SLS) method is used in this step. This method constructs one or a group of instrumental variables that are strongly correlated with the perturbation input but statistically independent of the noise and pressure state in the system, and uses the instrumental variables to eliminate endogenous interference and obtain unbiased and consistent estimates of the leakage parameters. In this way, even in the case that the small leakage signal is overwhelmed by the background noise of the closed-loop control, the estimation value of the structured leakage can be extracted relatively stably.

[0044] Optionally, according to an aspect of the present application, the steps of baseline identification under the premise of near constant pressure and residual sensitivity shaping can also be:

[0045] The original sequences of pressure P[k], mass flow rate m*[k], temperature T[k] and valve position u[k] are retrieved from the real-time buffer, and are de-biased and unit-normalized according to the calibration coefficients and time stamps of each sensor to obtain clean data stream Clean0 for modeling. In Clean0, the steady-state segment close to the target pressure P set is intercepted, and the drift points are removed by using the sliding window mean and variance threshold, and the ambient pressure P ext and the ambient temperature T ext are fixed, to obtain the steady-state baseline Baseline0 for subsequent shaping.

[0046] Further, inject small amplitude frequency sweep into steady baseline Baseline0, do cross-spectrum and phase regression on {P, m*} to identify liquid column / pipeline equivalent length L and sound speed a, and convert main mode frequency f1, to generate structure parameter package Struct1. Read structure parameter package Struct1 and current bandwidth of controller, select perturbation angular frequency band Ω d according to resonance avoidance principle, and set upper bound ε max for sensitivity of closed loop from mass flow rate m* to pressure P, form frequency band target set Band Target .

[0047] Further, inject frequency band target set Band Target into controller design environment, superimpose notch and feedforward / lag elements on pressure / volume compensation loop and perform phase margin check, output achievable shaping scheme Plan BRSS . Deliver Plan BRSS to execution controller, after online loading, retest closed loop response with perturbation sweep, get shaped sensitivity curve and record residual sensitivity ε(Ω d ). Compare residual sensitivity ε(Ω d ) with ε max , if constraint is not met, roll back parameters and automatically fine-tune notch center and depth; if condition is met, freeze configuration, get verified shaping state BRSS on . Optionally, run 1-2 reference periods with BRSS on , re-summarize steady state statistics of {P, m*, T, u} from latest Clean0, solidify as baseline snapshot Baseline1 (including Ω d , ε max , f1), as boundary and constraint for next step.

[0048] Optionally, in some embodiments, further comprising step four, based on the structured leakage estimation value, performing equivalent evaluation on energy storage efficiency of the constant pressure gas storage system; the calculation of the energy storage efficiency considers the leakage loss: η storage = (E out - ∫m leak · h· dt) / E in ; wherein E out is the released energy, E in is the stored energy, m leak is the estimated leakage flow rate, h is the specific enthalpy, and the integral interval is the entire energy storage period. Under constant pressure conditions, the energy loss rate caused by leakage can be represented as: λ loss = m leak• Cp • T / (P • V • y); where Cp is the specific heat at constant pressure, taken as 1005 J / (kg•K); T is the system temperature; P is the system pressure; V is the gas storage volume; and y is the specific heat ratio, taken as 1.4. For example, when the leakage flow rate is 0.001 kg / s, the system pressure is 10 MPa, the temperature is 300 K, and the volume is 100 m3, the energy loss rate λ 3 = 0.001 x 1005 x 300 / (10 x 10 loss x 100 x 1.4) = 0.0215% / hour. 6

[0049] The laboratory-scale leakage estimates are converted to equivalent leakage characteristics at the target engineering scale by a scaling relationship.

[0050] The scaling relationship is based on similarity criteria:

[0051] m leak,target = m leak,lab • (P target / P lab ) α • (D target / D lab ) 2 • sqrt(T lab / T target ); where the subscript target denotes the target engineering scale and lab denotes the laboratory scale; P is pressure, D is characteristic dimension, and T is temperature; and a is the pressure exponent, taken as 1.0 for laminar flow leakage and 0.5 for turbulent flow leakage. For the three coefficients of the structured leakage model, the scaling relationship is: target = θ1,lab• (D target / D lab ) 2 • sqrt(T lab / T target ); θ 2,target = θ2,lab• (D target / D lab ) 2 • sqrt(T lab / T target ) / P scale ; and θ3, target = θ3,lab• (D target / D lab ) 2 • sqrt(T lab / T target ) / sqrtP scale ; where P scale = P target / P lab is the pressure scaling ratio. ​

[0052] For example, when the laboratory pressure is 1 MPa, the characteristic size is 0.1 m, and the temperature is 300 K, the target engineering pressure is 10 MPa, the characteristic size is 1 m, and the temperature is 280 K, the scale-down factor of the leakage flow rate is: (10 / 1) 0.5 ×(1 / 0.1) 2 ×sqrt(300 / 280)≈103.8, that is, the leakage flow rate at the engineering scale is about 104 times the laboratory measurement value.

[0053] Preferably, the overall flow of a multi-mode constant pressure gas storage experiment equivalent simulation method can also be: determining a target perturbation frequency band; reading a perturbation signal and performing spectral pre-shaping to make the main energy fall into the target perturbation frequency band, injecting the perturbation signal after spectral pre-shaping into the system to obtain system response data; based on the system response data, estimating the structured leakage related to the system pressure online to obtain a structured leakage estimation value.

[0054] Embodiment three, this embodiment describes generating and injecting a perturbation signal to solve the technical problem of how to design an excitation signal with rich information and used for accurate parameter identification without interfering with the normal baseline of the system and exciting potential instability of the system.

[0055] In an optional implementation, the step of generating the perturbation signal includes: generating an initial perturbation signal, applying a pre-configured online small signal model to map the initial perturbation signal to a mass flow rate increment sequence; obtaining a real-time temperature sequence to calculate a corresponding enthalpy flow rate increment sequence; within a preset time window, applying zero integral constraints to the mass flow rate increment sequence and the enthalpy flow rate increment sequence to determine the perturbation signal satisfying the mass constraint and the enthalpy constraint and storing, as shown in Figure 2 .

[0056] Specifically, the online small signal model is a linearized model describing the relationship between the valve position change δu and the mass flow rate change δm* around the steady-state baseline (such as pressure P0, temperature T0, and valve position u0) of the current system operation, which is obtained by a small-range step or sweep experiment. The model can be a simple proportional gain model δm*[k]=K v •δu[k], or more accurately, a first-order or second-order transfer function model.

[0057] In this embodiment, in order not to introduce additional energy and mass accumulation effects to interfere with the constant pressure state of the system, the designed perturbation signal is within a preset time window W, and the cumulative sum of the mass flow increment δm*[k] caused by it should be zero, that is, it satisfies the mass constraint. In order to avoid introducing heat accumulation effects to interfere with the heat balance of the system, the cumulative sum of the enthalpy flow increment h(T[k])•δm*[k] caused by it should also be zero, that is, it satisfies the enthalpy constraint. The enthalpy function h(T[k]) here is obtained by querying the thermodynamic property table of the fluid (for example, air) or calling the property calculation function library according to the real-time collected temperature sequence T[k].

[0058] As a preferred embodiment, the step of simultaneously imposing a zero integral constraint on the mass flow increment sequence and the enthalpy flow increment sequence within the preset time window further comprises: constructing a time-weighted mass flow increment sequence; imposing a zero integral constraint on the time-weighted mass flow increment sequence within the preset time window, so that the perturbation signal satisfies the mass constraint, the enthalpy constraint and the time first moment constraint, as shown in Figure 3

[0059] The time first moment constraint here is to eliminate the low-frequency drift effect that the perturbation signal may introduce. A signal that only satisfies the mass and enthalpy double constraints, although its total integral quantity is zero, may continue to be positive in the first half of a window and negative in the second half. This unbalanced injection mode may cause a slow, difficult-to-recover transient deviation of the system state. By introducing the time first moment constraint, that is, requiring the window integral of the time-weighted mass flow increment k•δm*[k] (where k is the time sampling index) to be zero, the distribution of the perturbation on the time axis can be more balanced, with the action center located at the geometric center of the window. This effectively suppresses the low-frequency interference to the system operating baseline, so that the system can return to the steady state more quickly after the perturbation ends, and improves the signal-to-noise ratio of the experimental data and the estimation accuracy.

[0060] The determination process of the perturbation signal satisfying the mass constraint, the enthalpy constraint and the time first moment constraint described above specifically comprises the following steps:

[0061] generating a corresponding initial mass flow increment sequence based on a preconfigured signal template; using a constrained optimization solving method, solving a constrained mass flow increment sequence under the condition of satisfying the triple constraints; and applying a valve-flow inverse mapping model to inversely solve the constrained mass flow increment sequence into the perturbation signal, as shown in Figure 4

[0062] For example, a time-domain signal template U template , such as a symmetric double pulse signal with an amplitude of ±A and a pulse width of τ. Using the online small signal model described above, U template ​​Mapped to the initial mass flow rate increment sequence δm* init [k]. Construct a quadratic programming optimization problem with equality constraints: its objective function is to minimize the constrained sequence δm*. con [k] and the initial sequence δm* init The Euclidean norm between [k], i.e., min||δm* con -δm* init || 2 This ensures that the final signal, while satisfying the constraints, retains as much of the waveform characteristics of the initial template as possible. The constraints are the following three equations: Σ k∈W δm*[k]•Δt=0;Σ k∈W h(T[k])•δm*[k]•Δt=0; Σ k∈W k•δm*[k]•Δt=0; where: Σ is the discrete summation operator; k is the time sampling index; W is the index set of the perturbation window; δm*[k] is the mass flow increment generated by the perturbation at time k; Δt is the sampling period; h(T[k]) is the value of the enthalpy function at temperature T[k]; 0 is the zero value of the equality constraint; • is the multiplication operator; [] is the index symbol; ∈ indicates the set inclusion relationship.

[0063] By solving this optimization problem, the mass flow rate increment sequence δm* that satisfies the triple constraints can be obtained. con [k], apply the inverse mapping of the small-signal model (e.g., if the positive direction is δm*=K). v •δu, then the inverse mapping is δu=(1 / K v )•δm*), will δm* con [k] The inverse solution yields the final executable valve position perturbation signal U that satisfies the triple constraints. pair Optionally, the signal template U template It can also be a pseudo-random binary sequence (PRBS), a swept-frequency chirp signal, or other sequences with good autocorrelation characteristics. This invention is not limited to dual-pulse signals.

[0064] The perturbation signal U that satisfies the constraints is obtained. pair Subsequently, to ensure the safety and efficiency of the injection process, in this embodiment, the step of spectral preshaping of the perturbation signal includes: acquiring the system's main fluid mode frequency; and processing the perturbation signal using a filter to suppress the energy components of the perturbation signal at the main mode frequency and its harmonic frequencies, such as... Figure 5 As shown.

[0065] Specifically, the fluid main mode frequency f1(corresponding to the angular frequency Ω1) is the first-order resonance frequency of the fluid pressure wave in the system pipeline, which is related to the physical parameters such as the equivalent length of the pipeline and the fluid sound speed, and can be obtained through the baseline identification step of embodiment two. If the injected perturbation signal contains energy at the frequency f1or its integer harmonic (such as 2f1, 3f1) frequencies, the system resonance is easily excited, and severe pressure oscillation (i.e. water hammer effect) is generated, which endangers the safety of equipment. Therefore, it is necessary to design a digital filter to process U pair Preferably, the filter is a band-pass filter, and the frequency band selection needs to satisfy: Ω d ∈ [Ω low ,Ω high ] and {Ω low ,Ω high}∩{α•Ω1,2•Ω1,3•Ω1}={} d ; Wherein: Ω low is the perturbation main frequency band angular frequency; [Ω high ,Ω d ] is a closed interval; {} is a set;∩ is the intersection operation of the set; {} is an empty set; α is a bandwidth coefficient between 0.8 and 1.2; Ω1is the first-order resonance angular frequency; • is multiplication.

[0066] That is, the passband Ω max of the filter not only needs to be aligned with the sensitivity notch in embodiment two, but also needs to strictly avoid the first-order resonance frequency Ω1and its second and third-order harmonic frequencies. In one embodiment, a finite impulse response (FIR) filter with zero-phase filtering (such as the filt filt function) can be used to achieve this, so as to avoid introducing phase distortion.

[0067] Further, in order to make the signal accurately reproduced by the physical actuator, the method of the embodiment further comprises:

[0068] Before injecting the perturbation signal after spectral pre-shaping into the system, the perturbation signal is slope-limited according to the actuator slope upper limit; and the perturbation signal is beat-aligned according to the actuator minimum dwell time.

[0069] The actuator slope upper limit |Δu / Δt| Shape_draw is the maximum opening rate that the control valve can actually achieve. The filtered signal U maxThe minimum dwell time refers to the minimum time that the valve position command needs to stay in a new position in order to overcome the nonlinear effects such as spool sticking, friction, etc. The beat equalization process is to combine or broaden the pulses in the signal sequence whose duration is less than the dwell time, so that each valve position step is effective. Through the two steps of processing, the final perturbation signal can be physically executed by the actuator with high fidelity, ensuring the compliance of the experimental input and the theoretical design.

[0070] According to an aspect of the present application, the initial stage of the baseline identification and residual sensitivity shaping step can also include a cleaning and standardization process of the original data: after retrieving the original sequence of pressure P[k], mass flow m*[k] and the like from the real-time buffer, the accurate time domain alignment of the multi-channel signals is performed according to the factory calibration coefficients and real-time time stamps of the sensors, and unit normalization and sensor zero point deviation processing are performed, forming a clean data stream Clean0 that is time-synchronized and dimensionally consistent, which is relied on for all subsequent calculation steps.

[0071] According to an aspect of the present application, the process of determining the steady-state baseline can also specifically include screening using a sliding window statistics: in the Clean0 data stream, the mean and variance of the data within the sliding window calculation window are calculated. By setting reasonable mean drift threshold and variance fluctuation threshold, non-steady-state data points with process drift or severe fluctuation are automatically removed, and finally the steady-state data segment close to the target pressure P set is intercepted as the steady-state baseline Baseline0.

[0072] According to an aspect of the present application, the initial identification process of obtaining the main modal frequency can also specifically include active frequency sweeping and cross-spectrum analysis: in the stable operating point represented by Baseline0, a small-amplitude sweep excitation signal is jointly injected into the system. By performing cross-spectrum and phase regression analysis on the injected signal and the pressure P and mass flow m* responses of the system, the equivalent length L and fluid sound speed a of the system pipeline are identified, and the main modal frequency f1 is converted from them.

[0073] According to an aspect of the present application, the process of constructing the sensitivity notch and verification can also include an adaptive fine-tuning mechanism of parameters: after the calculated shaping scheme Plan BRSS is sent to the controller and loaded online, the perturbation sweep is retested to close the loop response. If the verification finds that the residual sensitivity ε(Ω d ) fails to meet the preset upper bound ε max , the system will automatically roll back the parameters and automatically fine-tune the center frequency and depth of the notch filter, and retest again until the constraint condition is met.

[0074] In embodiment three, the constraint generation on the perturbation signal can be divided into a basic implementation and a preferred implementation to embody the multi-level nature of the scheme.

[0075] As a basic implementation, the perturbation signal only needs to satisfy the mass constraint and the enthalpy constraint. The process of this double-constraint scheme is as follows: based on the initial signal template (such as a symmetric double pulse) and the online small signal model, an initial mass flow rate increment sequence δm init [k] is generated. An optimization problem with double constraints is constructed and solved, and the objective function is also to minimize ||δm con -δm init || 2 , but the constraint conditions only include the mass constraint Σ k∈W δm [k]•Δt=0 and the enthalpy constraint Σ k∈W h(T[k])•δm [k]•Δt=0. The δm con [k] solved is inversely mapped to the executable valve position perturbation signal through the inverse mapping model.

[0076] Through the double-constraint scheme, the injected perturbation achieves a balance in matter and energy to the system, i.e., within an injection period, there is no net increase or decrease in the working medium mass in the system, and there is no net increase or decrease in the total enthalpy value of the system. The beneficial effect is that after completing the detection task on the dynamic characteristics of the system, the perturbation does not cause cumulative drift of the steady-state working condition point (such as the average pressure, average temperature) of the system. However, in order to further suppress the more subtle transient baseline fluctuations that the perturbation may introduce, the present application further proposes a preferred triple-constraint scheme, i.e., adding a time first-moment constraint on the basis of the double constraint, as in the original text of embodiment three.

[0077] The process of inverse mapping to the perturbation signal in embodiment three is obtained through the following steps:

[0078] First, determine the upper limit of the perturbation injection. Before generating the perturbation signal, read the physical limit parameters of the actuator (such as the control valve), for example, the maximum allowed disturbance amplitude δu max and the slope upper limit |Δu / Δt| max . According to these parameters, calculate the perturbation energy level that does not trigger nonlinear effects (such as saturation, hysteresis), and generate an injection upper limit constraint package Act Limit This step makes the subsequently generated perturbation signal physically feasible, and the response relationship can be well described by the online small signal model.

[0079] Second, generate the initial template and map. Under the Act Limit constraint, based on the double-pulse time domain template generator, set the initial values of pulse width, interval and amplitude, and construct an unconstrained initial valve position perturbation template Utemplate The sequence of valve positions is mapped to an initial sequence of mass flow increments δm init [k] as initial values for the subsequent optimization solution.

[0080] Third step, perform constrained optimization solution. As before, a quadratic programming solver is employed to solve for the sequence of constrained mass flow increments δm init || 2 with the objective of minimizing ||δm con [k] as initial values for the subsequent optimization solution.

[0081] Fourth step, inverse mapping to complementary sequence. Apply the valve-flow inverse mapping approximation model to inverse solve the sequence of mass flow increments δm con [k] to a sequence of valve positions {δu + (k), δu - (k)} that can be executed by the actuators, and finally obtain the perturbation pair U pair The complementary sequence here has a clear physical meaning:

[0082] In the single-valve injection scenario, it can be a positive pulse δu + (k) and a compensatory negative pulse δu - (k) of equal total quantity (after model mapping) but opposite direction, so that the total disturbance effect is balanced.

[0083] In the dual-valve coordinated injection system, it can be an opening increment δu + (k) of the main excitation valve and an accurately calculated closing increment δu - (k) of the compensation valve.

[0084] Fifth step, constraint closure verification and fine tuning. The final generated perturbation pair U pair is mapped back to mass flow increments through the small-signal model and re-substituted into the three (or two) constraint integral equations to calculate the constraint residual R constraints . The constraint is verified to be closed precisely by comparing the residual with a pre-set zero threshold (e.g., 1e-6). If the residual is out of limit, the system can resort to a back-off strategy to fine-tune the amplitude or time duration of U template and re-execute the second to fifth steps until the constraint is satisfied within the tolerance range. This step guarantees the mathematical rigor of the final injection signal.

[0085] According to an aspect of the present application, after obtaining the system response data, the method can further include a data packaging process: in a complete perturbation injection cycle, the synchronized and time-aligned {P, m*, T, u} multi-channel signals are collected and packaged to form a data block D containing metadata tags inj . The tags can include injection timestamp, perturbation template type used, current system operating mode, etc., to facilitate subsequent batch processing and traceability analysis.

[0086] According to an aspect of the present application, before subsequent coherence evaluation or structured leakage estimation, a symmetric preprocessing process of the injection signal and the response signal can also be included: in order to eliminate errors introduced by signal processing itself, the collected injection signal sequence u and the response signal sequence m* need to be processed by the same band-pass filter. This processing process does not introduce additional phase distortion between the input and output signals, and a set of clean and phase-consistent band-pass signal pairs {u f [k],y f [k]} are obtained, which can be directly used for system identification.

[0087] Optionally, the steps of generating quality-enthalpy double-constrained paired perturbations and anti-resonance injection can also be:

[0088] Reading Baseline1 and valve group physical limits, combining with the actuator step test to obtain the allowed amplitude δu max and the slope upper limit |Δu / Δt| max , and calculating the perturbation energy level that does not trigger nonlinearity based on this, generating the injection upper limit Act Limit . Under the constraint of Act Limit , the complementary sequence {δu + (k), δu - (k)} is synthesized based on the double-pulse template of u, which is mapped to the mass flow increment δm*[k] through the online small signal model, and the mass and enthalpy integrals and the time first moment within the window are implemented with zero constraints, to obtain the original perturbation U pair that satisfies the quality-enthalpy double constraint.

[0089] Further, taking U pair and f1 of Struct1, a band-stop / band-pass filter is constructed to pre-shape the spectrum of the perturbation, so that the main energy falls in Ω d and is mutually exclusive with {f1, 2·f1, 3·f1}, and the anti-resonance sequence U shaped is output. U shaped is sent to the pre-execution check module, and the slope limit and beat alignment are performed according to |Δu / Δt| max and the minimum dwell time (dwell), to generate the schedule Schedule inj to avoid exciting the valve core stickiness and sub-harmonic.

[0090] Preferably, according to Schedule inj The perturbation injection is implemented under closed loop condition, and the {P, m*, T, u} are synchronously collected and time-aligned to form a tagged data block D inj , encapsulating a complete injection-response pair. The u of D inj and m* are band-pass filtered and coherence evaluated to eliminate segments with low SNR or insufficient coherence, obtaining a band-pass signal pair {u f [k], y f [k]} for estimation, along with a validity flag.

[0091] Embodiment Four, this embodiment details the online estimation of structured leakage based on system response data. It solves how to achieve accurate, continuous and robust online estimation of structured leakage under the complex conditions of multi-mode switching, closed loop feedback and measurement noise coexistence.

[0092] Before online estimation, in order to ensure the quality of data for estimation, the method of the present application can further comprise: before online estimation, evaluating the frequency domain coherence between the injected perturbation signal and the obtained system response data; and according to the coherence threshold, screening the valid data segment from the system response data, so as to use the valid data segment for online estimation.

[0093] Specifically, after a complete perturbation injection-response collection is implemented (as in Embodiment Three), the injected valve position perturbation sequence u[k] and the response quality flow sequence m*[k] (or pressure sequence P[k]) can be subjected to cross-spectrum analysis, and the frequency domain coherence function γ d (Ω 2 ) of the two in the target perturbation frequency band Ω d is calculated. The value range of this coherence function is [0, 1], and the closer the value is to 1, the higher the proportion of the component caused by the injected signal in the response signal, and the higher the signal-to-noise ratio (SNR). A coherence threshold can be preset, for example, 0.8. Only when the calculated γ 2 (Ω d ) is greater than the threshold, the current collected data segment is considered valid and is used for subsequent leakage estimation steps. If the coherence is insufficient, it may be due to excessive background noise or excessive system nonlinearity, and the data segment should be discarded and the amplitude or frequency band of the next round of perturbation injection should be considered for adjustment.

[0094] After the effective data segments are screened out, the step of online estimation of the structured leakage related to the system pressure comprises: identifying a current operation mode of the system based on the system response data, the operation mode comprising a charging mode, a discharging mode and a maintaining mode; and dividing the system response data into mode data windows corresponding to each operation mode according to the switching of the operation mode, for subsequent structured leakage estimation.

[0095] For example, the operation mode can be identified by monitoring the sign and magnitude of the mass flow m*[k]. When m*[k] continuously takes a large positive value, it is determined as the charging mode; when it continuously takes a large negative value, it is determined as the discharging mode; and when it fluctuates in a small range around zero, it is determined as the maintaining mode. Whenever a mode switching moment is detected, a mark is made on the data stream, and the long-period response data is naturally divided into a series of window sets {W c ,W d ,W h} with mode labels. Parameter estimation is performed within each window, and the subtle differences in leakage characteristics in different operation modes can be more accurately captured.

[0096] Within each mode data window, the structured leakage estimation comprises: constructing a pressure function vector for each pressure data point based on the system pressure and the ambient pressure within the mode data window; the pressure function vector comprises: a constant bias term, a term linearly related to the pressure difference between the system pressure and the ambient pressure, and a term related to the square root of the absolute value of the pressure difference.

[0097] The structured leakage model here is established based on the physical mechanism of fluid leakage. Specifically, for the system pressure P[k] and the ambient pressure P ext at each sampling moment k (usually atmospheric pressure), the constructed pressure function vector is Φ(P[k]) = [1, (P[k]-P ext ), sqrt(|P[k]-P ext |)] T . The first term 1 corresponds to a parameter for fitting a constant leakage or sensor bias independent of pressure; the second term (P[k]-P ext ) corresponds to a parameter for fitting a laminar flow leakage through a narrow slit, the flow of which is approximately linearly related to the pressure difference; and the third term sqrt(|P[k]-P ext |) corresponds to a parameter for fitting a turbulent flow leakage through an orifice, the flow of which is approximately proportional to the square root of the pressure difference. This structured model can more accurately describe the real leakage behavior than a single linear model, and improve the estimation accuracy.

[0098] Since the system is in a closed-loop control state, the perturbation input u fThe ordinary least squares (OLS) regression of [k] and the pressure function vector Φ(P[k]) as arguments will result in biased parameter estimates.

[0099] To solve this problem, the structured leakage estimation further comprises: generating a band-pass input signal based on the perturbation signal, and constructing a candidate instrumental variable set from the band-pass input signal; evaluating the statistical independence between each candidate instrumental variable in the candidate instrumental variable set and the pressure function vector; and selecting instrumental variables for the structured leakage estimation from the candidate instrumental variable set according to the evaluation result of the statistical independence.

[0100] It should be noted that a qualified instrumental variable Z needs to satisfy two conditions: 1) strong correlation with the endogenous explanatory variable (here mainly u f [k]) ; 2) no correlation with the model residual (which contains unmodeled dynamics related to Φ(P[k])).

[0101] Specifically, a candidate instrumental variable set is constructed based on the band-pass input signal u f [k], for example, it can include multiple different order delay versions of u f [k], such as {u f [k-1], u f [k-2],...}, and sequences generated by nonlinear transformations (such as square terms u f [k] 2 ) thereof.

[0102] Further, the statistical independence of each candidate instrumental variable and each column in the pressure function vector Φ(P) is evaluated by calculating the partial correlation coefficient or mutual information between them. The candidate variables with high coupling degree with Φ(P) are removed.

[0103] On this basis, the selected instrumental variables are made to have sufficiently strong correlation with u f [k] through F test or Cragg-Donald weak instrumental test. Finally, one or a group of instrumental variables that meet the conditions are obtained to form the instrumental variable matrix Z, which is used for subsequent two-stage least squares (2SLS) solution. The solution formula is: θ*=(Z T •U) -1 •Z T •Y; wherein: θ* is the parameter vector to be estimated, including the link gain and leakage model parameters; Z is the selected instrumental variable matrix; T is the transpose operator; U is the regression argument matrix, whose columns are composed of the band-pass input u f and the pressure function column vector Φ(P); -1 is the matrix inversion; Y is the observation vector, i.e. the band-passed mass flow output y f ; • is the matrix multiplication.

[0104] To ensure continuity and stability of the estimation results at mode switching, the step for subsequent structured leakage estimation further comprises: at the mode switching between the patterned data windows, starting a disturbance-free window to suspend the injection of the perturbation signal; defining a bridging window with the data surrounding the mode switching in adjacent patterned data windows; within the bridging window, correcting the leakage model parameters obtained from the previous patterned data window to achieve continuity of cross-mode estimation.

[0105] The step of correcting the leakage model parameters comprises: calculating a residual sequence within the bridging window using the leakage model parameters before correction; constructing an objective function based on the time integral of the residual sequence; adjusting the leakage model parameters by minimizing the objective function to obtain corrected leakage model parameters.

[0106] Specifically, when the system detects a mode switching event, the injection of the perturbation signal is immediately suspended, and a short disturbance-free window is entered to avoid the dramatic dynamic pollution of the statistics during the switching process. N b data points before and after the switching point are taken to form a bridging window. The parameters θ c estimated by the previous mode (e.g., charging mode) are used to calculate the model prediction residual e raw [k] within the bridging window. An objective function is constructed, such as the sum of squares of residuals J = Σ k∈Wb (e raw [k]) 2 A constrained optimization algorithm is used to make minor adjustments to the parameters θ c to minimize the objective function J, obtaining the corrected parameters θ c '. This set of corrected parameters θ c ' will be used as the initial value for parameter estimation of the next mode (e.g., holding mode), or directly as the smooth transition result at the mode boundary. This bridging consistency mechanism effectively avoids the estimation results jumping at the mode boundary, which is not consistent with the physical reality due to data windowing processing.

[0107] To quantify the reliability of the estimation results, structured leakage estimation further comprises: after obtaining the leakage model parameters and the corresponding residual sequence using instrumental variables; using a robust covariance estimation method based on heteroscedasticity and autocorrelation, calculating the confidence interval of the leakage model parameters based on the residual sequence.

[0108] In actual industrial data, the residual sequence of 2SLS regression often does not meet the ideal assumption of independent and identically distributed, often showing heteroscedasticity (i.e., the variance of the residual changes over time) and autocorrelation (i.e., the residuals at different times are correlated). The traditional covariance matrix calculation method will underestimate the uncertainty of the parameters.

[0109] Therefore, the present embodiment preferably adopts a covariance estimation method robust to heteroscedasticity and autocorrelation, such as the Newey-West estimator (a kind of HAC estimator), to calculate the covariance matrix of the parameter estimate θ*. Based on the robust covariance matrix, the confidence interval (e.g. 95% confidence interval) of each parameter can be calculated, which provides a key basis for the experimenter to judge the statistical significance and reliability of the leakage estimate.

[0110] According to an aspect of the present application, after the system response data is divided into patterned data windows, renumbering of the time axis within each window can also be included: for each divided patterned data window W c ,W d , or W h , the time axis index k within it is renumbered from 0. This is intended to stabilize the subsequent calculation of statistics that depend on the time index, avoiding numerical calculation problems that may be caused by the excessively large global time index value.

[0111] According to an aspect of the present application, the process of constructing a candidate instrument variable set from the band-pass input signal can also include various generation strategies: in addition to using the delayed version of the input signal, more diversified candidate instrument variables can be generated by applying a band-pass filter with a slightly shifted center frequency (i.e. bandwidth shift) or performing sign inversion and other transformations to the original band-pass input u f [k], to increase the probability of selecting strong and effective instrument variables from them.

[0112] According to an aspect of the present application, after the instrument variables used for structured leakage estimation are selected, a numerical condition optimization process for the final instrument variable matrix can also be included: in order to improve the numerical stability of the subsequent 2SLS solution, the matrix Z composed of the selected instrument variables can be further processed, such as by removing redundant columns (i.e. columns with high multicollinearity) or using orthogonal basis reconstruction (e.g. QR decomposition) to generate a final instrument matrix Z final with better numerical condition.

[0113] According to an aspect of the present application, before performing the two-stage least squares solution, standardization preprocessing of the data can also be included: in order to improve the numerical condition of matrix inversion and other operations, standardization processing (e.g. subtracting the mean and dividing by the standard deviation) or dimensionless normalization processing can be performed on each column of the regression independent variable matrix U and the observation vector Y, respectively.

[0114] According to an aspect of the present application, after obtaining the leakage model parameters and calculating the instantaneous leakage trajectory, a smoothing post-processing procedure can be further included: the original leakage estimation time series output by the 2SLS method can be further processed by a robust mean filter or a median filter to effectively suppress individual sharp pulses caused by measurement noise or model mismatch, and output a more smooth and physically more coherent leakage trajectory Leak seq .

[0115] According to an aspect of the present application, after obtaining the structured leakage estimation, the method can further include calculating and outputting a set of experimental effectiveness indicators: these indicators are used to quantitatively evaluate the quality of this estimation, which can specifically include the minimum detectable leakage threshold (i.e. the minimum leakage amount that the method can distinguish under the current signal-to-noise ratio threshold), the cross-mode deviation (i.e. the jump amplitude of the leakage estimation in different modes before bridge correction), and the variance of the estimated parameters, which are collectively summarized as an experimental effectiveness indicator report.

[0116] Optionally, the mode-aware structured leakage estimation and experimental effectiveness output (IV / 2SLS) can also be:

[0117] From {u f [k],y f [k]} and the sign and threshold events of the original m*[k], the working conditions are divided into window sets {W c ,W d ,W h} according to charging / discharging / holding, and the time axis is renumbered in each window to stabilize the statistics. The pressure trajectory P[k] is imported into each window and the pressure function vector Φ(P) = [1, (P-P ext ), sqrt(|P-P ext |)] is constructed, which is combined with the band-pass input u f to form the regression matrix U, and the band-pass output y f is assembled into the observation vector Y, obtaining the regression pair {U, Y} which can be used for parameter identification.

[0118] Further, according to Schedule inj and the time-delayed version of u f [k], the instrumental variable matrix Z is constructed, and through mutual information and orthogonality test, the columns strongly coupled with Φ(P) are removed to obtain a statistically more independent excitation subspace, so as to improve the estimation unbiasedness. U* is obtained by performing the first-stage projection on U using Z, and the second-stage least squares is completed by regressing Y with U* to output the link gain G* and the leakage parameter vector θ*, and the residual sequence is reserved for subsequent consistency constraints.

[0119] Further, {G*, θ*} is substituted back into {u f[k], Φ(P[k])} and generate leak estimate m* by sliding window averaging leak_hat and suppress spikes by robust mean / median filtering to produce coherent leak trajectory Leak seq . Monitor valve group switching and m* sign flipping, automatically freeze perturbation and open bridge window before and after switching, impose consistency regularization on residual integral to correct {G*, θ*} and obtain cross-mode aligned estimation Est aligned . Calculate minimum detectable leak threshold (minimum resolution at SNR threshold), cross-mode bias and estimation variance from Est aligned , aggregate into experiment equivalent indicator ExpEq Metrics and write back to report module.

[0120] Embodiment five, this embodiment solves the problem that the preset sensitivity notch and anti-resonance frequency band are invalid due to the drift of the physical characteristics of the system (especially the fluid main mode frequency) during long-term operation or working condition changes (such as fluid temperature, component changes).

[0121] To solve this problem, the method of the present application can further comprise: based on the system operation data, online estimating the fluid main mode frequency of the system to obtain an updated main mode frequency.

[0122] Specifically, this online estimation can be realized by using the disturbance-free window introduced in embodiment four. During mode switching, although the active perturbation injection has been suspended, the natural micro-fluctuation in the system or the transient response caused by the mode switching itself still contains the modal information of the system. At this time, the peak value of the main mode frequency f1 can be identified by performing passive spectral analysis (for example, calculating its power spectral density by using Welch method) on the pressure P[k] or mass flow m*[k] sequence in the disturbance-free window. In another alternative implementation, in order to obtain modal information with higher signal-to-noise ratio, a probe signal with extremely low energy and extremely short duration can be actively injected in the disturbance-free window, such as a short-time chirp signal or a low-amplitude white noise signal. Since the system is in a relatively quiet state at this time, the response of this probe signal can be clearly captured, and the updated main mode frequency f1 can be more accurately estimated. 1new This estimation process can be performed periodically at each mode switching to realize continuous tracking of the main mode frequency.

[0123] After obtaining the updated main modal frequency, the method of the present application further comprises: using the updated main modal frequency to update the sensitivity notch constructed by the spectral shaping step and the anti-resonance frequency band of the spectral pre-shaping step in linkage. In other words, the perturbation is automatically frozen, the statistics are restarted, and the main modal frequency is self-calibrated by short frequency sweep or passive spectrum analysis before and after mode switching, and the center of the sensitivity notch and the band-stop frequency band (perturbation frequency band) are updated in linkage to achieve the bidirectional enhancement of switching transient suppression and anti-resonance self-calibration.

[0124] In one aspect, the newly estimated main modal frequency f 1new is fed back to the spectral shaping module of Embodiment II. The system will update the spectral shaping module based on the new main modal frequency f 1new . d The perturbation frequency band Ω 1new is re-planned to maintain a safe distance from the new resonance region {α•Ω 1new , 2•Ω d ,...}. On this basis, the sensitivity notch is moved to the new perturbation frequency band Ω 1new by adjusting the controller parameters (such as the center frequency and depth of the notch filter) online, and the stability margin is rechecked. On the other hand, the main modal frequency f 1new is fed back to the spectral pre-shaping module of Embodiment III. The system will re-design the band-pass / band-stop filter for anti-resonance based on the new main modal frequency f 1new , so that the generated perturbation signal energy can be accurately injected into the updated sensitivity notch, and the energy components at the new resonance frequency and its harmonics can be effectively suppressed.

[0125] Through this online estimation-linkage update closed loop, the present application achieves the bidirectional enhancement effect of switching transient suppression-anti-resonance self-calibration. That is, mode switching provides a time window for online identification, and the results of online identification in turn enhance the safety and effectiveness of perturbation injection testing in various modes. Preferably, after each linkage update is completed, the system generates a parameter snapshot with a version number, recording the current f 1new , perturbation frequency band Ω and the corresponding controller and filter parameters, forming a traceable experimental version record, and improving the robustness and intelligence level of the entire experimental equivalent simulation method.

[0126] According to one aspect of the present application, when the leakage estimation adopts a recursive algorithm (such as recursive least squares RLS), the step of starting the disturbance-free window can also correspondingly include emptying and resetting the cumulative quantities of the algorithm: after detecting the mode switching boundary and triggering the disturbance-free window, the system will automatically empty the cumulative quantities in the recursive estimation algorithm, such as the information matrix and momentum vector. This prevents the statistical characteristics of the previous mode from polluting or dragging the initial estimation stage of the new mode, so that the algorithm can quickly converge at the beginning of the new mode.

[0127] Example 6. This example provides a concrete, reproducible numerical calculation procedure to illustrate the structured leakage estimation algorithm in Example 4, in particular the application of the two-stage least squares (2SLS) method based on instrumental variables.

[0128] Suppose a data set is collected within a patterned data window. For simplicity, only 5 sample points of data are shown.

[0129] The ambient pressure is set as P ext = 0.1 MPa.

[0130] The collected system pressure sequence P is: [10.10, 10.12, 10.08, 10.11, 10.09], (unit: MPa).

[0131] The perturbation input sequence u f after band-pass filtering is: [0.5, -0.5, 0.5, -0.5, 0.5], (unit: % valve opening).

[0132] The mass flow response sequence y f after band-pass filtering, (as the observation vector Y), is: [0.012, -0.008, 0.011, -0.009, 0.012], (unit: kg / s).

[0133] According to the structured model in Example 4, the pressure function vector Φ(P[k]) = [1, (P[k] - P ext ), sqrt(|P[k] - P ext |)] is constructed for each data point. T

[0134] k = 1: P[1] = 10.10 -> P - P ext = 10 -> Φ(P[1]) = [1, 10, 3.162] T .

[0135] k = 2: P[2] = 10.12 -> P - P ext = 10.02 -> Φ(P[2]) = [1, 10.02, 3.165] T .

[0136] k = 3: P[3] = 10.08 -> P - P ext = 9.98 -> Φ(P[3]) = [1, 9.98, 3.159] T .

[0137] k = 4: P[4] = 10.11 -> P - P ext = 10.01 -> Φ(P[4]) = [1, 10.01, 3.164]​T .

[0138] k=5: P[5]=10.09 -> P-P ext =9.99 -> Φ(P[5])=[1,9.99,3.161] T .

[0139] The linear model to be estimated is y f [k]=G•u f [k]+θ1•1+θ2•(P[k]-P ext )+θ3•sqrt(|P[k]-P ext |)+e[k].

[0140] Let θ*=[G*,θ 1_hat ,θ 2_hat ,θ 3_hat ] T be the vector of parameters to be estimated. G* is the link gain estimate, θ 1_hat is the constant leakage / bias term, θ 2_hat is the linear leakage coefficient, and θ 3_hat is the turbulence leakage coefficient.

[0141] The regression matrix of regressors U (size 5x4) is composed of the columns of u f and Φ(P): U=[[0.5,1,10,3.162],[-0.5,1,10.02,3.165],[0.5,1,9.98,3.159],[-0.5,1,10.01,3.164],[0.5,1,9.99,3.161]]; the observation vector Y (size 5x1) is: Y=[0.012,-0.008,0.011,-0.009,0.012] T .

[0142] Assume that u f is correlated with the residual e[k] (endogeneity) and that Φ(P) is uncorrelated with the residual (exogeneity). Instrumental variables need to be found for the endogenous variable u f .

[0143] Construct the instrumental variables, specifically, select the first order lag of u f , u f [k-1] as the instrumental variable. Assume that u f [0]=0, then Z u =[0,0.5,-0.5,0.5,-0.5] T . For the exogenous variables Φ(P), they are their own best instrumental variables. Thus, the complete instrumental variable matrix Z (size 5x4) is:

[0144] Z = [[0, 1, 10, 3.162], [0.5, 1, 10.02, 3.165], [-0.5, 1, 9.98, 3.159], [0.5, 1, 10.01, 3.164], [-0.5, 1, 9.99, 3.161]].

[0145] In the first stage, regress each column of U on all columns of Z to obtain the fitted values of U, U*. The formula is U* = Z • (Z T •Z) -1 •Z T •U. Compute Z T •Z (4x4 matrix), find its inverse (Z T •Z) -1 , and then compute Z T •U (4x4 matrix). Multiply them together to get U* (5x4 matrix). The numerical computation of matrix multiplication and inversion is a well-known technique and will not be discussed here.

[0146] In the second stage, use U* as the new regressors to perform ordinary least squares (OLS) regression on Y to obtain the final parameter estimates θ*. The formula is θ* = (U* T •U*) -1 •U* T •Y. Again, first compute U* T •U* (4x4 matrix), find its inverse, then compute U* T •Y (4x1 vector), and multiply them together to get the final parameter vector θ*.

[0147] After the above calculations, suppose the resulting parameter vector is θ* = [0.0201, 0.0002, -0.0001, 0.0005] T . This means:

[0148] The link gain estimate G* = 0.0201 ((kg / s) / (% valve position)). The constant leakage / bias term θ 1_hat = 0.0002 (kg / s).

[0149] The linear leakage coefficient θ 2_hat = -0.0001 ((kg / s) / MPa). The turbulent leakage coefficient θ 3_hat = 0.0005 ((kg / s) / MPa 0.5 ).

[0150] Therefore, the resulting structured leakage model is:

[0151] m* leak_hat (P) = 0.0002 - 0.0001 • (P - 0.1) + 0.0005 • sqrt(P - 0.1).

[0152] For example, when the system pressure is 10.1 MPa, the estimated instantaneous leakage is:

[0153] m* leak_hat = 0.0002 - 0.0001 • (10) + 0.0005 • sqrt(10) ≈ 0.0002 - 0.001 + 0.00158 ≈ 0.00078 kg / s.

[0154] When performing equivalent simulation experiments, the operating parameters (such as pressure, flow rate, temperature) of the laboratory bench often differ by orders of magnitude from the parameters of the target actual scene. In order to make the conclusions obtained in the laboratory effectively guide the actual engineering, a set of scale mapping criteria need to be defined in advance. In one application scenario of the present application, the mapping relationship can be defined as: P target (t) = λ P • P lab (t); m target (t) = λ m • m lab (t); T target (t) = λ T • T lab (t); where: P target (t) is the pressure trajectory of the target scene; P lab (t) is the pressure trajectory of the experimental bench; m target (t) is the mass flow trajectory of the target scene; m lab (t) is the mass flow trajectory of the experimental bench; T target (t) is the temperature trajectory of the target scene; T lab (t) is the temperature trajectory of the experimental bench; λ P is the pressure scale factor; λ m is the flow scale factor; λ T is the temperature scale factor; t is time; • is multiplication.

[0155] The scale factors λ P , λ m , λ T are constants determined in advance based on similarity criteria (such as Reynolds number similarity, Mach number similarity, etc. fluid mechanics criteria) and system energy balance relationships. Through this mapping relationship, the structured leakage model m* lab (P lab ) estimated by the method of the present application at the experimental bench pressure P leak_target can be converted to the equivalent leakage model m* target (P eq ) in the target actual scene. This makes it possible to evaluate the leakage characteristics of large-scale energy storage systems in a cost-controllable laboratory environment, improving research and testing efficiency.

[0156] One of the objectives of the method of the present application is to achieve an equivalent simulation of constant pressure gas storage experiments, in order to make the key physical quantities trajectories (such as pressure) of the experimental bench accurately reproduce the pre-set reference target trajectories. In order to quantify the fidelity of this reproduction, an experimental equivalence comprehensive index E eq .

[0157] In one particular embodiment, the index can be calculated by the following formula:

[0158] E eq = w P • mean(|P lab (t) - P ref (t)| / P ref ,max) + w m • mean(|m lab (t) - m ref (t)| / m ref ,max) + w T • mean(|T lab (t) - T ref (t)| / T ref ,max) + w S • mean(|dP lab / dt - dP ref / dt| / S ref );

[0159] wherein: E eq is the experimental equivalence comprehensive error; w P , w m , w T , w S are non-negative weights and sum to 1; mean is the time average operator; |·| is the absolute value; P lab (t), m lab (t), T lab (t) are the experimental bench trajectories; P ref (t), m ref (t), T ref (t) are the reference target trajectories generated from the aforementioned scale mapping; P ref ,max, m ref ,max, T ref ,max are the normalization factors of each reference quantity; dP lab / dt, dP ref / dt are the pressure time derivatives; S ref is the reference slope normalization constant; / is division; - is subtraction; + is addition; • is multiplication; t is time.

[0160] The comprehensive index Eeq The smaller the value is, the higher the equivalence of the experimental simulation is. It not only considers the static deviation of physical values such as pressure, flow rate, temperature, etc., but also considers the compliance of the dynamic response of the system through the pressure change rate term |dP lab / dt-dP ref / dt|, and the compliance of the dynamic response of the system. The accurate online leakage estimation realized by the present application plays an important role in reducing E eq . For example, the estimated leakage amount m* leak_hat can be input into the pressure controller of the system as a feedforward compensation term to offset the adverse effects of leakage on the system pressure, so that P lab (t) can more closely track P ref (t), obtain a lower E eq value, and obtain more reliable experimental equivalent simulation results.

[0161] According to one aspect of the present application, the sensitivity notch configuration method is specifically as follows, which is used to solve the technical problem of how to realize spectrum shaping in a closed-loop constant pressure control system to improve the signal-to-noise ratio of the perturbation signal response.

[0162] The process of determining the target perturbation frequency band specifically includes obtaining the open-loop transfer function G(s) and the controller transfer function C(s) of the system, and calculating the closed-loop sensitivity function S(jω)=1 / (1+G(jω)C(jω)). Wherein, j is the imaginary unit, ω is the angular frequency, and s is the Laplace operator. The fluid main mode frequency f1 and its harmonic frequencies 2f1, 3f1 of the system are determined by a sweep frequency experiment, and these resonance points are marked on the frequency axis. When the target perturbation frequency band Ω d needs to meet the constraint condition: |Ω d -n·Ω1|>0.2·Ω1, where n=1, 2, 3, and Ω1=2πf1 is the main mode angular frequency. For example, when f1=5Hz, Ω d can be selected within the range of [8, 12]Hz.

[0163] The specific implementation of configuring the sensitivity notch in the controller includes: connecting a notch filter N(s) in series with the original pressure controller C(s) to form a modified controller C'(s)=C(s)·N(s). The notch filter adopts a second-order band-stop structure: N(s)=(s 2 +2ζ1ω n s+ω n 2 ) / (s 2 +2ζ2ω n s+ω n 2 ); wherein, ω nζ1 is the center frequency of the notch, set as the center of the target perturbation frequency band; ζ2 is the numerator damping ratio, ranging from 0.01 to 0.05 to form a deep notch; ζ3 is the denominator damping ratio, ranging from 0.5 to 1.0 to ensure stability; s is the Laplace operator. The notch depth can be controlled by adjusting the ratio of ζ1 / ζ2, with a typical value of 0.02 / 0.7 yielding a sensitivity improvement of 20 dB.

[0164] Stability margin verification was performed using the Nyquist criterion. The open-loop transfer function L'(jω)=G(jω)C'(jω) after adding the notch filter was calculated, and its Nyquist plot was plotted. The curve was ensured not to encircle the point (-1,j0), and the shortest distance to this point was greater than 0.3, corresponding to a gain margin of 6dB. Simultaneously, the phase margin PM=180°+arg(L'(jω)C ... c ))>30°, where ω c Let L' be the crossover frequency, satisfying |L'(jω) c If the margin is insufficient, increase ζ2 or decrease the notch depth and redesign. In actual controller implementation, it can be discretized using a digital filter, and the sampling frequency should be greater than 10 times the notch center frequency.

[0165] According to one aspect of this application, the process of constructing and selecting instrumental variables, and the implementation of structured leakage estimation using the two-stage least squares method, are as follows:

[0166] A bandpass input signal is generated based on the perturbation signal, and a set of candidate instrumental variables is constructed from the bandpass input signal. The construction of candidate instrumental variables includes the following methods: The first type is time-delay instrumental variables, which take the bandpass input signal u... f Multiple delayed versions of [k] Z1={u f [k-1],u f [k-2],...,u f [k-10]}; The second type is the frequency shift tool variable, which applies a bandpass filter with a center frequency shift of ±0.5Hz to uf[k] to obtain Z2={u f_shift1 [k],u f_shift2 [k]}; The third category is nonlinear transformation instrumental variables, calculating Z3={u f 2 [k],|u f [k]|,sign(u f [k])·sqrt|u f [k]|}. The three categories are merged to form a candidate set Z. candidate .

[0167] The specific method for assessing the statistical independence between candidate instrumental variables and the stress function vector is as follows: for each candidate instrumental variable z iand each column φ of the pressure function vector Φ(P) j Calculate the partial correlation coefficient r ij =cov(z i ,φ j ) / sqrt(var(z i )·var(φ j )); where cov is the covariance operator and var is the variance operator. The independence threshold is set to 0.15. If |r ij If |>0.15, then z is considered to be i With φ j There is a strong correlation, z i Remove from the candidate set. Perform a weak instrumental test on the retained instrumental variables and calculate the F-statistic for the first-stage regression: F = (Ri)i 2 / (1-R 2 ))·((nk-1) / k); where R 2 is the coefficient of determination for the first-stage regression, where n is the sample size and k is the number of instrumental variables. An F > 10 is required to classify a variable as a strong instrumental variable.

[0168] The specific solution process for two-stage least squares is as follows: In the first stage, the endogenous variable u... f Regressing the instrumental variable matrix Z yields the fitted value u. f_hat =Z(Z'Z) -1 Z'u f In the second stage, the augmented regression matrix U=[u f_hat ,Φ(P)], where Φ(P)=[1,(PP ext ),sqrt|PP ext |]' represents the pressure function vector. Solve for the parameter vector θ=(U'U) -1 U'y yields the link gain G* and leakage coefficient [θ1,θ2,θ3]. For example, when the system pressure is 10 MPa and the ambient pressure is 0.1 MPa, the structured leakage flow rate is calculated as follows:

[0169] m leak =θ1+θ2·(10-0.1)+θ3·sqrt(10-0.1)=θ1+9.9θ2+3.146θ3(kg / s).

[0170] According to one aspect of this application, the process of triple-constraint optimization solution and bridging window correction algorithm is as follows, mainly used to solve the problems of numerical stability and cross-mode estimation continuity in constraint optimization.

[0171] The process involves generating a perturbation signal that satisfies the triple constraints of mass-enthalpy-time first moment. The specific solution process is as follows: Construct the Lagrangian function L:

[0172] L = ||δm* - δm init || 2 + λ1 · (Σδm[k]Δt) + λ2 · (Σh(T[k])δm*[k]Δt) + λ3 · (Σk · δm*[k]Δt);

[0173] where λ1, λ2, λ3 are Lagrange multipliers.

[0174] Convert the problem into KKT conditions for solving, form a linear equation group: [2I, A'; A, 0] [δm*; λ] = [2δm init ; 0]; where I is the unit matrix, A is the constraint matrix, and its three rows correspond to the three constraint conditions respectively. Use QR decomposition to solve the equation group, when the condition number of the matrix is greater than 100, introduce the regularization term ε = 1e-6, modify it to (2I + εI) to improve the numerical condition.

[0175] If the constraints are incompatible, resulting in no solution, then relax the constraints according to the priority, relax the time first moment constraint to |Σk·δm[k]Δt|<0.01, relax the enthalpy constraint to |Σh(T[k])δm*[k]Δt|<0.1, and the mass constraint is always strictly satisfied. The window length W is typically taken as 200-500 sampling points, corresponding to 2-5 seconds of injection time.

[0176] The specific algorithm for correcting the leakage model parameters in the bridging window: the bridging window size Nb is taken as 30 data points before and after the mode switching. Define the residual sequence e[k]=y[k]-G·u[k]-θ'Φ(P[k]), where θ is the parameter to be corrected. Construct the objective function J(Δθ)=Σ k∈Wb (e[k]-ΔθΦ(P[k])) 2 +ρ||Δθ|| 2 ; where Wb is the bridging window, Δθ is the parameter correction amount, and ρ=0.01 is the regularization coefficient. By taking the derivative and setting it to zero, the correction amount is obtained: Δθ=(Φ'Φ+ρI) -1 Φ'e; the corrected parameter θ new =θ old +Δθ. To ensure physical reasonableness, limit the correction amplitude |Δθ i / θ i_old |<0.2. The jamming window length is set to 50 sampling points (0.5 seconds) to ensure that the mode switching transient is fully attenuated.

[0177] Optionally, the determination process of the coherence threshold is as follows: through Monte Carlo simulation, evaluate the relationship between the coherence function and the estimation accuracy under different noise levels. When the expected estimation error is less than 5%, the coherence function γ 2 >0.75. Considering the engineering margin, set the threshold to 0.8. For different frequency bands, the threshold can be adjusted adaptively: γ2 threshold = 0.8 - 0.1(f / f1), where f is the current frequency, f1 is the dominant modal frequency, ensuring a slightly lower coherence when far from the resonance point.

[0178] The perturbation signal amplitude range is exemplified as follows: through step response experiments, the linear working interval of the valve is determined to be [20%, 80%] opening. The perturbation amplitude δu max is set to be 30% of the distance from the current steady-state working point to the linear boundary. For example, when the steady-state working point u o = 50%, δu max = min(50-20, 80-50) x 0.3 = 9%. Meanwhile, the pressure fluctuation constraint |ΔP / P o | < 0.5% must be satisfied, through the system gain G p is estimated to be δu max < 0.005P o / G p . The smaller value of the two constraints is taken as the final limit.

[0179] The real-time guarantee is achieved through algorithm optimization: recursive least squares (RLS) is used instead of batch 2SLS, and the information matrix update formula is P [k] = (P [k-1] - P [k-1] z [k] z’ [k] P [k-1] / (1+z’ [k] P [k-1] z [k] )) / λ f ; where P is the information matrix, z is the instrumental variable, and λ f = 0.995 is the forgetting factor. The parameter update is θ [k] = θ [k-1] + P [k] z [k] . The computational complexity is reduced from O(n 3 ) to O(p 2 ), where p is the number of parameters. At a sampling rate of 2 kHz, the time consumed by a single update is less than 0.1 ms, meeting the real-time requirement.

[0180] According to one aspect of the present application, the process of Newey-West covariance estimation is exemplified as follows, used to calculate the confidence interval of the leakage parameter.

[0181] The covariance evaluation method is robust to heteroscedasticity and autocorrelation, and the calculation formula of the Newey-West estimator is: Σ NW = (X’X) -1 · Ω · (X’X) -1Where X is the regression matrix and Ω is the long-term covariance matrix. The estimate of Ω is: Ω = Σ o +Σ j=1 L w j (Σ j +Σ j '); where Σ j =(1 / n)Σ n t=j+1 e t 2 x t x' t-j Let e ​​be the j-th order autocovariance. t For the residual, x t w is the regression variable. j =1-j / (L+1) represents the Bartlett kernel weights; L=floor(4(n / 100)) (2 / 9) ) represents the bandwidth parameter. For a sample of n=1000, L≈5.

[0182] Parameter θ i The 95% confidence interval is calculated as: CI i =θ i ±1.96·sqrt(Σ NW [i,i]); where Σ NW [i,i] represents the diagonal elements of the covariance matrix. For example, when θ²=0.001 and the standard error SE=0.0002, the confidence interval is [0.0006,0.0014]. If the confidence interval contains zero, it indicates that the leakage component is not statistically significant and can be removed from the model to simplify the structure.

[0183] This invention, from the design of high-fidelity perturbation signals to robust online estimation algorithms under multi-mode and closed-loop conditions, and then to an adaptive parameter calibration mechanism, ultimately achieves accurate, continuous, and reliable online identification of pressure-related structured leaks in constant-pressure gas storage systems, providing important technical support for the experimental verification of novel energy storage technologies.

[0184] This invention achieves optimal performance in a specific frequency band Ω by actively shaping the sensitivity function of the pressure closed-loop controller. dThe inner structure sensitivity notch is designed, and the corresponding main energy perturbation signal is also accurately fallen into the frequency band, and the problem that the excitation signal is controlled to be disappeared is solved. The principle is that the technology is not in competition with the constant voltage controller, but by locally changing the frequency response characteristic, a solution is provided for the detection signal of a specific frequency. The energy of the injected perturbation signal can bypass the strong suppression of the controller, and the system response can be reflected with very high signal-to-noise ratio, and the disturbance suppression capability of the system in all other frequency bands is not affected. In the constant voltage energy storage experiment, in the scene requiring high pressure stability, the method can effectively extract the dynamic characteristics submerged by weak leakage without sacrificing the control performance of the system, and provides a high-quality data basis for subsequent accurate and unbiased estimation.

[0185] The present application generates a perturbation signal by constructing and solving an optimization problem with mass, enthalpy and time first moment triple constraints, realizes online detection of zero disturbance of the energy storage system. Specifically, the zero integral constraints of mass and enthalpy ensure that the perturbation signal will not introduce net increase of working medium mass or heat to the system within an injection cycle, avoiding cumulative interference to the system steady-state pressure and temperature baseline. And further time first moment zero constraint, by ensuring the balanced distribution of disturbance on the time axis, suppresses the low-frequency drift of system state that may be caused by disturbance asymmetry. In the long-time constant pressure experiment, the purity of experimental data is realized, and the system pseudo-response introduced by the test behavior itself is avoided to be misjudged as leakage characteristics, and the confidence of online diagnosis is improved.

[0186] The present application solves the problem of biased estimation of closed-loop system by introducing a pattern-aware instrumental variable (IV) two-stage least squares (2SLS) estimation algorithm, and combining a structured leakage model based on physical mechanism. Through pattern recognition, the data is divided into windows, and a structured model containing laminar flow and turbulent flow terms is used to improve the fitting accuracy of the model. By designing a tool variable that is strongly correlated with the real disturbance, but statistically independent of the system closed-loop noise and state variables, the 2SLS algorithm can effectively break the variable correlation caused by system endogeneity, eliminate the contaminated part of the excitation signal, and only use the clean part for parameter regression. This makes the algorithm can accurately identify the real causal relationship from pressure change to leakage flow from complex closed-loop data, obtain unbiased and consistent leakage parameter estimation, and its result can truly reflect the physical health status of the system.

[0187] The application solves the discontinuity problem of cross-mode estimation by designing a disturbance-free window / bridge window and a parameter correction mechanism based on residual integral minimization. During the system mode switching process, the method suspends the perturbation injection to avoid data pollution, and uses the data on both sides of the switching point to form a bridge window. By fine-tuning the parameters of the previous mode, the prediction error in the bridge window is minimized, realizing the smooth transition of model parameters at the mode boundary. The final output of the leakage estimation is a continuous trajectory with time and working condition changes, which has physical meaning, rather than a series of jumps and broken fragments at the switching point without reasonable explanation. For energy storage systems that need long-term continuous monitoring of health status, the continuity of the estimation is the basis for realizing advanced diagnostic functions such as fault trend prediction and residual life assessment.

[0188] The application solves the potential problem of poor robustness of existing fixed parameter methods under variable working conditions by establishing a self-calibration coupling closed loop between online estimation of resonance frequency and spectral shaping parameters, and obtains adaptive ability to environmental and working condition changes. In the long-term operation of energy storage systems, changes in temperature, pressure and other factors will cause the drift of fluid sound speed and other physical parameters, causing changes in system resonance frequency f1. The application continuously identifies the latest value of f1 online, and links it back to the sensitivity notch design module and the perturbation signal anti-resonance filter design module at the front end, realizing the accurate matching of the excitation channel and safety protection with the current real dynamic characteristics of the system. This adaptive ability enables the entire high-precision identification scheme proposed by the application to maintain its effectiveness and safety without human intervention for a long time, enhancing its reliability and automation level in practical engineering applications.

[0189] The above describes the preferred embodiments of the application in detail, but the application is not limited to the specific details in the above embodiments. Within the technical concept of the application, various equivalent transformations of the technical solutions of the application can be made, and these equivalent transformations all belong to the protection scope of the application.

Claims

1. A multi-mode constant-pressure gas storage experimental equivalent simulation method, characterized in that, This includes the leakage characteristic assessment process, specifically: Determine the target perturbation frequency band; Read the perturbation signal and perform spectrum preshaping to make its main energy fall into the target perturbation frequency band. Inject the spectrum-preshaped perturbation signal into the system and obtain the system response data. Based on system response data, structured leakage related to system pressure is estimated online, and the estimated value of structured leakage is obtained.

2. The method according to claim 1, characterized in that, Generating perturbation signals includes: An initial perturbation signal is generated, and a pre-configured online small-signal model is applied to map it into a mass flow increment sequence; Obtain the real-time temperature sequence to calculate the corresponding enthalpy-flow increment sequence; Within a preset time window, zero integral constraints are simultaneously applied to both the mass flow rate increment sequence and the enthalpy flow rate increment sequence to determine and store the perturbation signals that satisfy the mass and enthalpy constraints.

3. The method according to claim 2, characterized in that, Within a preset time window, zero integral constraints are simultaneously applied to both the mass flow rate increment sequence and the enthalpy flow rate increment sequence, including: Construct a time-weighted mass flow increment sequence; Within a preset time window, a zero-integral constraint is applied to the time-weighted mass flow increment sequence to ensure that the perturbation signal satisfies the mass constraint, enthalpy constraint, and first-order time moment constraint.

4. The method according to claim 3, characterized in that, To ensure that the perturbation signal satisfies quality constraints, enthalpy constraints, and first-order time moment constraints, including: Generate the corresponding initial mass flow rate increment sequence based on the pre-configured signal template; A constrained optimization solution method is adopted to obtain the constrained mass flow rate increment sequence under the condition of satisfying triple constraints; By applying the valve-flow inverse mapping model, the constrained mass flow rate increment sequence is decomposed into a perturbation signal.

5. The method according to claim 1, characterized in that, Spectral preshaping of perturbation signals includes: Obtain the dominant fluid modal frequencies of the system; A filter is used to process the perturbation signal in order to suppress the energy components of the perturbation signal at the main mode frequency and its harmonic frequencies.

6. The method according to claim 1, characterized in that, Online estimation of structured leaks related to system stress, including: Based on system response data, the current operating mode of the system is identified, including inflation mode, discharge mode, and holding mode. Based on the switching of operating modes, the system response data is divided into patterned data windows corresponding to each operating mode for subsequent structured leakage estimation.

7. The method according to claim 6, characterized in that, Used for subsequent structured leakage estimation, including: At the mode switching point between patterned data windows, activate the do-not-disturb window to pause the injection of perturbation signals; Define the bridging window by using data around the mode switch in adjacent patterned data windows; Within the bridging window, the leakage model parameters obtained from the previous patterned data window are corrected to achieve continuity in cross-pattern estimation.

8. The method according to claim 7, characterized in that, The parameters of the leakage model were corrected, including: Using the leakage model parameters before correction, calculate the residual sequence within the bridging window; A target function is constructed based on the time integral of the residual sequence. The leakage model parameters are adjusted by minimizing the objective function to obtain the corrected leakage model parameters.

9. The method according to claim 6, characterized in that, Within each patterned data window, the structured leakage estimate includes: Based on the system pressure and environmental pressure within the patterned data window, a pressure function vector is constructed for each pressure data point. The pressure function vector includes: a constant bias term, a term linearly related to the pressure difference between the system pressure and the ambient pressure, and a term related to the square root of the absolute value of the pressure difference.

10. The method according to claim 9, characterized in that, Structured leakage estimation also includes: A bandpass input signal is generated based on the perturbation signal, and a set of candidate instrumental variables is constructed from the bandpass input signal; Evaluate the statistical independence between each candidate instrumental variable in the candidate instrumental variable set and the stress function vector; Based on the assessment results of statistical independence, instrumental variables for structured leakage estimation are selected from the set of candidate instrumental variables.

Images