Electrolytic cell power feed-forward control method for off-grid hydrogen production system

Abstract

The invention discloses an electrolytic cell power feed-forward control method for an off-grid hydrogen production system, which comprises the following steps: acquiring high-frequency pressure data through a pressure sensor array, and analyzing and generating flow-electricity coupling physical information based on the high-frequency pressure data; performing spatio-temporal evolution prediction on the local pressure field by using the flow-electricity coupling physical information and the high-frequency pressure data to generate predicted local pressure field evolution data; and according to the predicted local pressure field evolution data, a feedforward control instruction used for adjusting the input power of the electrolytic cell is generated. According to the method, accurate prediction and active inhibition of the flow-electricity coupling instability in the electrolytic cell are realized, and the operation stability and the energy utilization efficiency of the off-grid hydrogen production system are improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of power optimization control, and particularly relates to an electrolytic cell power feedforward control method for an off-grid hydrogen production system. BACKGROUND

[0002] With the accelerated promotion of global energy transformation, water electrolysis for hydrogen production using renewable energy has become a key technical path to achieve carbon neutralization. Off-grid hydrogen production systems can directly consume fluctuating renewable energy such as wind power and photovoltaic power, avoiding grid connection and power transmission losses, and have unique advantages in remote areas and offshore wind farms. However, the random fluctuation characteristics of renewable energy power pose challenges to the stable operation of electrolytic cells. The electrochemical reaction inside the electrolytic cell is strongly coupled with the gas-liquid two-phase flow, and power fluctuations can cause complex flow-electricity coupling instability, leading to problems such as pressure oscillation, bubble aggregation, and local overheating. Therefore, developing an advanced control method that can accurately perceive and actively suppress this instability is of great significance for improving off-grid hydrogen production systems.

[0003] Currently, electrolytic cell power control mainly uses feedback control strategies based on macroscopic parameters. A typical method is to monitor system-level parameters such as total current, total voltage, outlet pressure, and electrolyte temperature of the electrolytic cell, and perform PID adjustment when these parameters deviate from the set values. Some studies have introduced model predictive control (MPC), which uses an equivalent circuit model and thermodynamic model of the electrolytic cell for multi-variable coordinated control. In terms of signal acquisition, existing systems generally use low-frequency sampling (1-10 Hz) to obtain operating data through SCADA or DCS systems. A few studies have attempted to install pressure sensors at the outlet of the electrolytic cell for vibration monitoring, but still remain at the level of single-point measurement and frequency spectrum analysis. In terms of state estimation, key parameters such as gas holdup are mainly obtained through empirical formulas or lookup tables, and there is a lack of in-depth modeling of the complex physical processes inside the electrolytic cell. The control period is usually fixed at the second level, matching the power scheduling period of renewable energy.

[0004] The core problem of the existing technology can be summarized as follows: the lack of accurate perception ability of the millisecond-level flow-electricity coupling dynamics inside the electrolytic cell leads to control lag and insufficient precision. SUMMARY

[0005] The application aims to provide an electrolytic cell power feedforward control method for an off-grid hydrogen production system to solve the above problems in the existing technology.

[0006] Technical solution, an electrolytic cell power feedforward control method for an off-grid hydrogen production system, comprising:

[0007] Collecting high-frequency pressure data from a pressure sensor array of the electrolytic cell;

[0008] Based on high-frequency pressure data, analyze and generate flow-electricity coupling physical information;

[0009] Using the flow-electricity coupling physical information and the high-frequency pressure data, the spatio-temporal evolution of the local pressure field is predicted to generate predicted local pressure field evolution data;

[0010] According to the predicted local pressure field evolution data, feedforward control instructions for adjusting the input power of the electrolytic cell are generated.

[0011] Beneficial effects, the present application realizes the accurate prediction and active inhibition of the flow-electricity coupling instability in the electrolytic cell, and improves the operation stability and energy utilization efficiency of the off-grid hydrogen production system. BRIEF DESCRIPTION OF DRAWINGS

[0012] Figure 1 A step flowchart of an electrolytic cell power feedforward control method for an off-grid hydrogen production system provided by the embodiments of the present application.

[0013] Figure 2 A step flowchart of analyzing and generating flow-electricity coupling physical information provided by the embodiments of the present application.

[0014] Figure 3 A step flowchart of the acquisition method of the gradient components of the pressure gradient tensor provided by the embodiments of the present application.

[0015] Figure 4 A step flowchart of the characteristic decomposition of the pressure gradient tensor provided by the embodiments of the present application. DETAILED DESCRIPTION

[0016] In order to enable persons 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 with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, but not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by persons skilled in the art without creative labor should fall within the scope of protection of the present application.

[0017] It should be noted that the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device including a series of steps or units does not have to be limited to only those steps or units clearly listed, but can include other steps or units not clearly listed or inherent to the process, method, product or device.

[0018] In the study, it is found that the spatial resolution of pressure field monitoring is seriously insufficient. The existing method only arranges a small number of sensors at the inlet and outlet of the electrolytic cell, and cannot capture the three-dimensional pressure field distribution inside, especially the vertical pressure gradient information is completely missing. When a simple one-dimensional wave equation assumption is used to estimate the vertical gradient, the effects of bubble buoyancy and hydrostatic pressure are ignored, which will produce a larger error in the high gas holdup area. The estimation of key state parameters lacks physical basis. The existing gas holdup empirical formula only considers the influence of current density, and cannot reflect the internal relationship between bubble size distribution, oscillation characteristics and gas holdup, resulting in a larger estimation error under dynamic working conditions, and further affecting the accurate calculation of sound speed field and pressure propagation characteristics. In addition, the time scale mismatching problem of the control strategy is prominent. The millisecond-level pressure fluctuation information is averaged by the second-level control period, and the opportunity for preventive control of transient instability is lost. Only passive response can be made after the problem accumulates to the macro level, resulting in frequent power adjustment and efficiency loss.

[0019] As shown in Figure 1 , an electrolytic cell power feedforward control method for an off-grid hydrogen production system is proposed, comprising the following steps:

[0020] Collect high-frequency pressure data from the pressure sensor array of the electrolytic cell.

[0021] In this embodiment, a 3x3 pressure sensor array is arranged at the key area of the electrolytic cell, such as the inner wall of the cavity near the gas outlet. These sensors collect pressure signals in real time at a high sampling rate of, for example, 1 kHz. At the same time, the system also obtains low-frequency data such as wind speed and wind power from the SCADA system of the wind farm, and medium-frequency operating data such as current, voltage, and temperature from the DCS system of the electrolytic cell. The collected high-frequency pressure raw signals are processed by, for example, a Butterworth low-pass filter (cutoff frequency 200 Hz) to remove high-frequency electrical noise, and the filtered high-frequency pressure data for subsequent analysis are obtained.

[0022] Based on the high-frequency pressure data, flow-electricity coupling physical information is analyzed and generated. The flow-electricity coupling physical information is used to characterize the flow-electricity coupling state inside the electrolytic cell.

[0023] Specifically, the flow-electricity coupling state refers to the complex interaction between the flow field dynamics such as bubble generation, movement, and convergence caused by current changes in the electrolysis process, and the physical fields such as pressure waves and electric field distribution. In this embodiment, the flow-electricity coupling physical information is obtained through a series of complex analysis and calculation, the core of which is to construct a pressure gradient tensor that can describe the spatial variation of the pressure field, and to estimate the local gas holdup closely related to the bubble behavior. These information together constitute a multi-dimensional state vector for quantitatively characterizing the instantaneous flow-electricity coupling state.

[0024] The flow-electricity coupling physical information and high-frequency pressure data are used to predict the time-space evolution of the local pressure field, and evolution data of the predicted local pressure field is generated.

[0025] In this embodiment, based on the physical information and high-frequency pressure data, a high-resolution local pressure field at the current time is reconstructed by physical constraint interpolation and virtual sensor technology. Modern dynamic system analysis methods such as dynamic mode decomposition (DMD) are applied to extract the dominant dynamic modes from the historical pressure field snapshot sequence. Each mode contains its inherent evolution frequency and growth rate. By time extrapolation of these modes according to their own rules and linear superposition, the time-space evolution data of the local pressure field in a short time window (for example, 10-50 ms) in the future can be predicted.

[0026] According to the predicted evolution data of the local pressure field, a feedforward control instruction for adjusting the input power of the electrolytic cell is generated.

[0027] Specifically, the generation of the feedforward control instruction is a multi-objective and multi-time scale optimization decision process. Based on the prediction data, multiple objectives such as pressure stability, electrolysis efficiency, and renewable energy utilization rate are optimized. The millisecond-level pressure mutation, the hundred-millisecond-level pressure oscillation, and the second-level steady-state trend are comprehensively considered, and an event-triggered mechanism is used for updating. The generated instruction is matched with the response characteristics of the actuator (such as power supply, water pump) to form a power and pump speed adjustment instruction that can be directly issued, realizing fine feedforward control of the electrolytic cell.

[0028] In a specific application scenario, the present embodiment can be applied to an off-grid hydrogen production system connected with renewable energy sources such as wind power. The core equipment of the system is the electrolytic cell, and the input power needs to be quickly and accurately adjusted according to the fluctuation of renewable energy sources to ensure the hydrogen production efficiency and equipment safety.

[0029] As shown in Figure 2 According to one aspect of the present application, the flow-electricity coupling physical information is analyzed and generated, including:

[0030] Based on the high-frequency pressure data, a pressure gradient tensor is constructed for the local three-dimensional space in the electrolytic cell; wherein the flow-electricity coupling physical information is represented by the pressure gradient tensor.

[0031] For example, the pressure gradient tensor is a three-dimensional matrix defined at each sensor position, which includes the pressure value itself (zero-order information), the first-order derivative of the pressure in three spatial directions (gradient information), and the second-order derivative information.

[0032] Further, as shown in Figure 3 When constructing the pressure gradient tensor, the gradient component of the pressure gradient tensor is obtained in the following manner:

[0033] By spatially differencing the high-frequency pressure data, a horizontal gradient component of the pressure gradient tensor is resolved;

[0034] Based on the horizontal gradient component, an indirect method based on fluid physical model is used to calculate the vertical gradient component;

[0035] The indirect method is to establish a hydrostatic equilibrium relation involving the mixed phase density and the gravity effect, and combine the fluid vertical velocity indirectly calculated from the divergence of the horizontal gradient component to solve the vertical gradient component;

[0036] The hydrostatic equilibrium relation expresses the vertical gradient component as a function of the mixed phase density, the gravity acceleration and the fluid vertical velocity indirectly calculated from the divergence of the horizontal gradient component.

[0037] In the present embodiment, the horizontal gradient component is obtained directly. For a 3x3 sensor array arranged on a horizontal plane z=z0, the pressure gradients in x and y directions can be calculated by numerical methods such as central difference. For example, for a central measurement point (i,j), its x-direction gradient ΨP / Ψx can be approximated by (P[i+1,j] - P[i-1,j]) / (2Δx), where P[i,j] is the pressure value of the measurement point, Δx is the sensor spacing in x direction, and Ψ is the partial derivative. Since the sensors are usually arranged on the same plane, vertical differentiation cannot be directly performed. A simple assumption of vertical propagation of pressure waves, for example, using ΨP / Ψz ≈ (1 / c) x ΨP / Ψt to estimate the vertical gradient component, will introduce huge errors in such a complex three-dimensional flow field as an electrolytic cell, because it ignores the effects of hydrostatic pressure and bubble motion. Therefore, an indirect calculation method based on physical model is used. A hydrostatic equilibrium relation considering gravity, fluid acceleration and mixed phase density is established. A preferred implementation is that the vertical pressure gradient ΨP / Ψz is determined by the following relation: ΨP / Ψz = -ρ mix × g - Ψ(ρ mix × w 2 ) / Ψz; where ΨP / Ψz is the gradient of the local pressure in the vertical direction z; ρ mix is the mixed phase density of the gas-liquid two-phase flow, which can be calculated by ρ mix =(1-α) × ρ liquid + α × ρ gas , where α is the local gas holdup, ρ liquid is the liquid phase density, and ρ gas is the gas phase density; g is the gravity acceleration; w is the vertical velocity component of the fluid; and Ψ(ρ mix × w 2) / Ψz is the dynamic pressure gradient term associated with the vertical motion of the fluid. In this relation, the vertical velocity w of the fluid is also unknown, but can be indirectly calculated by the continuity equation using the divergence of the horizontal velocities. Specifically, Ψw / Ψz ≈ -(Ψu / Ψx + Ψv / Ψy), while the horizontal velocity components u and v can be calculated from the horizontal pressure gradients ΨP / Ψx and ΨP / Ψy by the simplified momentum equation. By coupling the solutions, the physically more reliable estimate of the vertical pressure gradient ΨP / Ψz can be finally obtained.

[0038] The embodiment realizes accurate perception of the three-dimensional pressure field under the constraint of sensor plane arrangement. Errors caused by the simple assumption of vertical propagation of pressure waves are avoided, and in particular in the complex three-dimensional flow field caused by the rising of bubbles in the electrolytic cell, the vertical pressure gradient generated by the combined action of bubble buoyancy and fluid inertia can be accurately captured. Compared with the traditional scheme requiring multi-layer sensor arrangement, the complete characterization of the three-dimensional pressure field is realized only with a single layer of nine sensors, reducing the system complexity and cost, and improving the perception accuracy of the transient flow-electric coupling phenomenon in the electrolytic cell, laying a foundation for subsequent accurate prediction and control.

[0039] According to one aspect of the present application, as Figure 4 shown, when analyzing and generating flow-electric coupling physical information, it also includes feature decomposition of the pressure gradient tensor, specifically:

[0040] Perform singular value decomposition on the matrix form of the pressure gradient tensor at each time; extract the first N principal singular values and the corresponding characteristic mode vectors from the decomposition results, where N is a natural number greater than 0.

[0041] In the embodiment, the feature decomposition preferably uses singular value decomposition (SVD), which can decompose a complex matrix into a linear combination of multiple simple, mutually orthogonal modes. Specifically, the pressure gradient tensor at each time t is reorganized into an observation matrix M obs . For example, the gradient information (such as 0-order, 1-order, and 2-order derivatives, a total of 27 values) of 9 measuring points can be organized into a 27x9 observation matrix. Perform SVD on the observation matrix M obs , and obtain M obs =UΣV T ; where T is the transpose, U is the left singular vector matrix, and its column vector u k represents the basic mode or structure of the pressure field spatial distribution, which is called the characteristic mode vector; Σ is a diagonal matrix, and the elements σ k on its diagonal are singular values, arranged in descending order, representing the corresponding mode u kEnergy or importance weight in the overall pressure field; V T The transpose of the right singular vector matrix, represents the weight distribution of the mode on different sensors. After the decomposition is completed, a preset number of the top N principal singular values (for example, N = 3), i.e., σ1, σ2, σ3, and the corresponding characteristic mode vectors u1, u2, u3, are extracted from the results. These principal modes and principal singular values capture the most important coherent structure and energy distribution of the pressure field. By using SVD for eigen-decomposition, seemingly chaotic, high-dimensional pressure gradient information can be decomposed into several dominant modes with clear physical meaning. For example, one mode can clearly represent the bubble rising flow along the vertical direction, and another mode can represent the pressure wave propagating in the horizontal direction. This provides a basis for understanding and determining the complex flow-electricity coupling state from a macroscopic perspective.

[0042] According to one aspect of the present application, when analyzing and generating flow-electricity coupling physical information, it further includes estimating the local void fraction, specifically:

[0043] Performing time-frequency spectrum decomposition on the high-frequency pressure data to identify and extract the frequency spectrum characteristics caused by bubble oscillation;

[0044] Applying the Minnaert theoretical model based on bubble resonance physics to invert the local void fraction in the electrolytic cell by taking the frequency spectrum characteristics as input.

[0045] In the present embodiment, the local void fraction, i.e., α, refers to the volume fraction of hydrogen and oxygen bubbles generated by electrolysis in a local microelement volume in the electrolytic cell. It is a key parameter for describing the characteristics of gas-liquid two-phase flow. Specifically, the short-time Fourier transform (STFT) is performed on the high-frequency pressure data P filtered [i, j, t] of each sensor to obtain its time-frequency spectrum. For example, the calculation can be performed using a window length of 0.1 seconds and an overlap rate of 50%. By analyzing the spectral density in the frequency band of 10-100 Hz, the characteristic frequency generated by the collective resonance of a large number of small bubbles can be identified, denoted as f bubble (for example, its range is usually 20-50 Hz). After obtaining f bubble , the purely empirical formula lacking theoretical basis is abandoned, and the Minnaert resonance theory based on physics is used for void fraction inversion. A single bubble in a liquid has an inherent resonance frequency, which is related to its size and surrounding environmental parameters. For a two-phase flow containing a large number of bubbles, its resonance characteristics will be affected by the void fraction. In one preferred embodiment, the local void fraction α resonance is inversely calculated by the following formula based on the Minnaert theory: resonance α liquidx (2p x R b x f bubble ) 2 ); wherein a resonance is the gas holdup estimated by resonance method; g is the adiabatic index of gas (for hydrogen-oxygen mixture, an approximate value can be taken); P0 is the local static pressure; p liquid is the density of electrolyte; R b is the equivalent bubble radius, which can be further estimated by the relationship between pressure fluctuation amplitude and frequency, for example, R b ~ (c liquid / (2p x f bubble )) x sqrt (AP / P0), wherein c liquid is the sound speed in pure liquid phase, AP is the pressure fluctuation amplitude corresponding to the frequency band; f bubble is the bubble resonance characteristic frequency identified from the frequency spectrum.

[0046] Further, the flow-electricity coupling state vector is constructed based on the principal singular values and the characteristic mode vectors, specifically:

[0047] The spatial distribution characteristics of the principal singular values and the characteristic mode vectors, as well as the local gas holdup and the spectral characteristics, are integrated to form a multi-dimensional flow-electricity coupling state vector; wherein the flow-electricity coupling physical information is quantitatively characterized by the flow-electricity coupling state vector.

[0048] Specifically, in order to comprehensively and quantitatively describe the flow-electricity coupling state, a multi-dimensional state vector S couple (t) is constructed. The vector is a data structure that integrates information from different analysis dimensions. An exemplary S couple (t) can contain 12 dimensions of information, specifically composed as follows: S couple (t) = [M couple , C conf , s1, s2, s3, a mean , c min , c max , s p , f bubble , f wave , E ratio ]; wherein M couple is the coupling mode label determined by fuzzy logic rules, for example, 1 represents the bubble rising dominant mode, 2 represents the pressure wave propagation mode, and 3 represents the strong coupling mode; C conf is the confidence of the mode determination; s1, s2, s3 are the first three principal singular values; a mean is the average local gas holdup estimated in the sensor array region; c min and c maxThe minimum and maximum values ​​of the estimated local sound velocity field are σ. p f represents the standard deviation of the pressure signal, characterizing the overall intensity of the pressure pulsation; bubble and f wave These are the characteristic frequencies of bubble oscillation and pressure wave propagation identified through spectral analysis, respectively; E ratio This represents the energy ratio between the two frequency bands. By constructing a comprehensive state vector, the complex and unobservable current-electric coupling phenomena inside the electrolyzer are transformed into a high-dimensional, quantitative, and information-rich mathematical description. This state vector can serve as a standard input for subsequent prediction models, control algorithms, or fault diagnosis systems, improving the accuracy and reliability of subsequent processing.

[0049] In a specific embodiment of the current-electric coupling state vector construction, a fuzzy logic reasoning system is built, with the input variable being: x1: energy proportion of the first principal mode r1 = σ1 2 / Σσ k 2 x2: Verticality of the first principal mode v1 = |u1(z)|, i.e., the absolute value of the vertical component of the first characteristic mode vector u1; x3: Isotropy of the first three principal modes iso = (σ1+σ2+σ3) / (3σ1); x4: Energy ratio of bubble to pressure wave E ratio = E bubble / E wave E bubble E represents the energy of bubble oscillation. waveFor pressure wave energy, define a fuzzy set and its membership function for each input variable. For example, for x1 (energy percentage), define the fuzzy set {low, medium, high}, whose membership function can be a trapezoidal or Gaussian function. For example, high can be defined as a membership of 1 when x1 > 0.7, linearly decreasing between 0.5 and 0.7, and a membership of 0 when x1 is less than 0.5. Based on expert knowledge and experimental data, establish a set of IF-THEN rules to determine the coupling mode: Rule 1 (determines the bubble rising dominant mode M=1): IF (x1 is high) AND (x2 is high) THEN (mode is bubble rising dominant); the principle is that when the energy of the first dominant mode is absolutely dominant, and the spatial structure of this mode exhibits strong vertical motion, it indicates that the flow field is mainly controlled by the upward motion of bubbles along the direction of gravity. Rule 2 (Determining Pressure Wave Diffusion Mode M=2): IF (x1 is low) AND (x3 is high) THEN (mode is pressure wave diffusion); the principle is that when no mode dominates (low energy percentage), and the energy distribution of the first few modes is relatively uniform (high isotropy, close to 1), it indicates that the pressure field exhibits pressure waves that diffuse uniformly in all directions, rather than structured flow. Rule 3 (Determining Strongly Coupled Mode M=3): IF (x1 is medium) AND (x4 is in the medium range) THEN (mode is strongly coupled); the principle is that when multiple modes coexist, and the bubble oscillation energy in the spectrum is comparable to the pressure wave energy (E... ratio When the value approaches 1), it indicates a strong interaction and energy exchange between the flow structure and the pressure wave propagation. For any set of inputs {x1, x2, x3, x4}, the above rules may be activated to varying degrees. The final mode determination output M can be obtained using a weighted average unfuzzy resolution method (e.g., the centroid method). couple Simultaneously, the confidence level C can be calculated based on the strength of the activated rule. conf For example, if the activation strength of rule 1 is much higher than that of other rules, then M couple =1 and C conf Approaching 1. At this point, the 12-dimensional state vector S can be clearly explained. couple The integration process of (t). It is not a simple listing, but an orderly assembly of information from different sources and different processing stages: [M couple C conf [σ1, σ2, σ3] are calculated in real time by the fuzzy logic system; [σ1, σ2, σ3] are obtained directly by SVD decomposition; [α] mean c min c max The value is calculated from the gas holdup and the sound velocity field, where α mean It is α local Spatial average, c min and c max It is Clocal the maximum value in the field; [σ p , f bubble , f wave , E ratio ] are the original features extracted directly from the time-frequency spectrum analysis.

[0050] This embodiment converts the originally chaotic multi-point pressure measurement data into a state description with clear physical meaning, where the principal singular values reflect the energy distribution of different physical modes (such as bubble rising, pressure wave propagation), and the eigenvectors reveal the spatial structure of these modes. It can accurately identify the flow-electricity coupling state changes inside the electrolytic cell on the millisecond time scale, and detect the instability germination several seconds earlier than the traditional monitoring method based on macroscopic parameters (total current, total pressure), which provides valuable response time for feedforward control.

[0051] Optionally, to further improve the reliability of local gas holdup estimation, a verification and fusion mechanism based on acoustic attenuation is also introduced, including:

[0052] According to the attenuation law of sound waves in gas-liquid two-phase flow, the spectral features are analyzed to obtain the acoustic attenuation gas holdup as an independent verification source;

[0053] The resonance method gas holdup inverted from the Minnaert theoretical model is weighted and fused with the acoustic attenuation gas holdup to generate the cross-verified local gas holdup.

[0054] To overcome the errors that may exist in a single model, this embodiment introduces estimation methods based on different physical principles as cross verification. This method is based on the principle of acoustic attenuation, that is, the energy of sound waves will attenuate when passing through gas-liquid two-phase flow, and the degree of attenuation is closely related to the gas holdup. Specifically, by analyzing the signals received by sensors arranged at different positions (for example, sensors A and B) for the same acoustic event, the energy attenuation of sound waves propagating between points A and B can be calculated. According to the attenuation model of sound waves in gas-liquid two-phase flow, an independent gas holdup estimate value can be inverted from the attenuation, denoted as α acoustic . The gas holdup α resonance obtained by the resonance method is weighted and fused with the gas holdup α acoustic obtained by the acoustic attenuation method to obtain the final, more reliable local gas holdup α local . The fusion process can be represented as: α local = w res ×α resonance + w acou × α acoustic ; where w res and w acou are weight coefficients, and w res + w acou= 1. Illustratively, w res = 0.7, w acou = 0.3. These weights can also be adaptively adjusted according to the estimated confidence in the current operating condition. By fusing the results of two estimations based on different physical principles, the local gas holdup is more robust and accurate, effectively reducing the risk of failure or deviation of a single model.

[0055] This embodiment realizes high-precision estimation of local gas holdup in the electrolytic cell by applying Minnaert bubble resonance theory and combining acoustic attenuation verification mechanism for weighted fusion. Compared with empirical formulas lacking theoretical basis, it can accurately reflect the internal relationship between bubble size, oscillation frequency and gas holdup. Especially in the scene of electrolytic hydrogen production, which is a violent gas-liquid two-phase flow, accurate estimation of gas holdup is crucial for predicting pressure fluctuations and optimizing electrolytic efficiency. Through cross-validation of dual physical mechanisms, the estimation error of gas holdup is reduced, and the accuracy of subsequent sound speed calculation and pressure field reconstruction is improved.

[0056] According to one aspect of the present application, it also includes calculating an adaptive local sound speed field, specifically:

[0057] The local gas holdup, as well as the preset density and sound speed parameters of the liquid and gas phases, are jointly substituted into the gas-liquid two-phase flow acoustic model based on the Wood formula for operation to obtain an adaptive local sound speed field.

[0058] Specifically, the local gas holdup α local can be obtained with high reliability. Then, the adaptive sound speed field varying with space and time in the electrolytic cell can be calculated. Here, adaptive means that the sound speed is no longer a fixed constant, but a variable that can dynamically reflect the change in local bubble concentration. The Wood formula is preferably used to calculate the sound speed field, specifically: c local = 1 / sqrt((α / c gas 2 + (1-α) / c liquid 2 )×(α×ρ gas + (1-α)×ρ liquid ));where c local is the local mixed-phase sound speed; α is the local gas holdup α local ; c gas and c liquid are the sound speeds in pure gas and pure liquid phases, respectively; ρ gas and ρ liquid are the densities of pure gas and pure liquid, respectively. Alternatively, the equivalent simplified form is c local [i,j,t] = c water / sqrt(1 +α×(ρwater / ρ gas -1) / (1-α)), which is in ρ liquid >> ρ gas This approximately holds true under certain conditions. For example, when α = 0 (pure water), c local Approximately 1480 m / s; while when α = 0.3, c local It will drop sharply to about 500 m / s. Where c water Let ρ be the speed of sound in pure water. water The density is that of pure water. The goal is to obtain a sound velocity field distribution that accurately reflects the real-time operating conditions inside the electrolyzer, varying with time and space.

[0059] According to one aspect of this application, it also includes performing time synchronization calibration on high-frequency pressure data using an adaptive local sound velocity field as a physical constraint, specifically:

[0060] By performing cross-correlation analysis on high-frequency pressure data from different sensor channels, the initial time delay between sensor pairs is determined, forming an initial time delay set.

[0061] An adaptive local sound velocity field is applied as a spatiotemporal consistency constraint, and the initial time delay set is optimized to obtain an optimized time delay that satisfies the laws of physical propagation.

[0062] The clock correction curves of each sensor are calculated based on the optimized time delay, and the clock correction curves are applied to the original high-frequency pressure data to generate time-synchronized pressure data.

[0063] In this embodiment, although the sensor array may be triggered by the same clock system, during long-term operation, minute delay differences and clock drift between channels (including sensors, cables, and data acquisition cards) accumulate, resulting in clock deviations on the order of nanoseconds to milliseconds. Therefore, high-precision time synchronization calibration is required. Specifically, for any pair of sensors (i, j), the high-frequency pressure data P acquired by them is calculated. i (t) and P j The cross-correlation function R of (t) ij (τ)= ∫P i (t) ×P j (t-τ)dt。 R ij (τ) The value of τ corresponding to the peak value is the initial measurement delay τ. init (i, j) reflects the time required for a pressure event to propagate between the two sensors. Initial delay τ init This may include measurement noise and clock skew. To obtain the most physically consistent true time delay, two types of spatiotemporal consistency constraints are introduced for optimization: a sound speed physical constraint: the propagation time of a pressure wave between two points should be equal to its distance divided by the average sound speed. An adaptive local sound speed field c is then used.local A constraint can be established: |τ(i, j) - d ij / c avg (i, j)| < ε c ; where τ(i, j) is the time delay to be optimized; d ij is the physical distance between sensors i and j; c avg (i, j) is the average sound speed on the path between the two points, which can be calculated from c local ; ε c is a small tolerance error, for example, 2 ms. This constraint eliminates abnormal time delay values that do not meet the physical propagation rule. Closed loop path constraint: for any closed loop triangle formed by three sensors (i, j, k), the time delay should satisfy the path consistency, i.e., τ(i, j) + τ(j, k) + τ(k, i) ≈ 0 (considering the direction). More intuitively, τ(i, j) + τ(j, k) ≈ τ(i, k). Combining all possible triangular closed loop constraints in the sensor array (for 9 sensors, there are C(9, 3) = 84 in total), an over-determined linear equation system Aτ= b + ε can be constructed; where τ is the vector containing all the time delays to be solved, A is the coefficient matrix describing the path relationship, b is the ideal value (usually 0), and ε is the measurement error. By solving this equation system by Weighted Least Squares, the globally optimal, internally consistent optimized time delay matrix T opt is obtained. When solving, higher weights can be given to those time delays that are highly consistent with the physical constraints of sound speed. After obtaining the optimized time delay T opt , by comparing it with the initial measured time delay T measured , the clock bias of each sensor can be calculated Δt[i] = T opt [i, ref] - T measured [i, ref] (relative to the reference sensor ref). Low-pass filtering is performed on this bias sequence to smooth the drift, and a smoothed clock correction curve Δt smooth [i, t] is obtained. The correction is applied to the original data: P sync [i, t] = P filtered [i, t - Δt smooth [i, t]], so as to obtain the time-synchronized pressure data P sync . By using the physical characteristics (adaptive sound speed) of the flow field itself as a natural ruler, self-calibration time synchronization is achieved.

[0064] This embodiment calculates the spatiotemporally varying sound velocity field using Wood's formula and uses it as a physical constraint for sensor time synchronization calibration. By introducing variable sound velocity constraints and closed-loop path consistency constraints, sub-millisecond clock synchronization accuracy is achieved. Utilizing the inherent physical characteristics of the flow field as a natural clock overcomes the problem of external clock sources being susceptible to interference in industrial environments, providing a reliable time reference for high-precision pressure gradient calculations.

[0065] According to one aspect of this application, it further includes reconstructing the local pressure field within a preset local three-dimensional space, specifically:

[0066] Acquire pressure values ​​from measurement points based on time-synchronized pressure data;

[0067] Obtain the measurement point gradient values ​​from the pressure gradient tensor;

[0068] Establish physical constraints to describe the characteristics of gas-liquid two-phase flow in an electrolyzer;

[0069] By combining the pressure values ​​at the measuring points, the gradient values ​​at the measuring points, and the physical constraints, a local pressure field is generated through a joint solution using an interpolation algorithm based on radial basis functions.

[0070] The physical constraints include at least the continuity constraint describing the conservation of fluid mass and the simplified momentum equation constraint describing the force balance of the flow field.

[0071] In this embodiment, after obtaining precisely synchronized pressure data, a dense local pressure field P, defined within a three-dimensional space (e.g., a 0.5m × 0.5m × 0.2m cuboid covering the sensor array) and containing 75 grid points, is recovered from sparse information from 9 physical measurement points. field_local To achieve high-precision reconstruction, an interpolation method integrating multi-source information is adopted, preferably radial basis function (RBF) interpolation. RBF is particularly suitable for handling non-meshable sparse data. The joint solution process specifically involves selecting a radial basis function, such as a quadratic function Φ(r) = sqrt(r). 2 +c 2 The pressure field P(x) can be expressed as P(x) = Σλ, where r is the distance and c is the shape parameter. i ×Φ(|x - x i |) + p1x + p2y + p3z + p4; where x is the coordinate of any point in space; x i Let λ be the coordinates of the i-th physical sensor; i p1 to p4 are the coefficients to be solved. Simultaneous multi-source information is used as constraints: Pressure value constraint at each physical measuring point x... i At that point, the interpolation result must be equal to the synchronous pressure measurement value at that point, i.e., P(xi ) = P sync [i] constitute the first set of equations. Measurement point gradient value constraints: at each physical measurement point x i , the gradient of the interpolated result ▽P(x i ) is equal to the value of the pressure gradient tensor at that point ▽P tensor [i] constitute the second set of equations. Physical constraint conditions: the pressure field should satisfy the basic laws of fluid dynamics at all grid points within the reconstruction region. These laws, after being discretized, serve as additional linear equations, at least including: continuity constraint: ▽(p v) = S gas , which describes the conservation of fluid mass, where ▽ is the gradient, v is the fluid velocity vector, and S gas is the gas mass source term per unit volume due to electrolytic reaction; simplified momentum equation constraint: ▽P = -p g - F drag , which describes the force balance of the flow field, where p is the mixed phase density, g is the gravity, and F drag is the resistance experienced by the bubble. All the above equations are combined to form a large linear equation system. By solving this equation system, all unknown coefficients (l i and p j ) can be obtained. Once the coefficients are determined, the pressure value at any grid point within the reconstruction region can be calculated, thereby generating the local pressure field P field_local [k, t]. The generated local pressure field is not only accurate at the measurement points but also maintains physical reasonableness and continuity of the flow field structure in the region between the measurement points, with much higher accuracy and physical fidelity than traditional interpolation methods that only use pressure values.

[0072] Optionally, in some more refined physical models, the physical constraint conditions describing the characteristics of the fluid can be enhanced. For example, the traditional momentum equation may not be able to fully capture the complex memory effect of bubbles in gas-liquid two-phase flow (i.e., the motion of bubbles depends not only on the current force but also on their historical trajectory). Fractional-order derivatives can be introduced to construct more accurate physical constraint equations. For example, the simplified momentum equation ▽P = -p g - F drag is improved to D α P / Dt α = -p g - F drag - F memory , where D is the derivative operator, a is the order of the fractional derivative, P is the pressure field, and F memory is the memory term, representing the additional force due to the historical trajectory and force state of the bubble or fluid particles. This can more accurately describe complex dynamic processes with memory and genetic characteristics, further improving the physical fidelity of pressure field reconstruction.

[0073] According to one aspect of this application, the reconstructed pressure field has high accuracy in regions close to the physical sensors, but accuracy decreases in regions far from all sensors (such as corners or centers of the array). Therefore, it also includes identifying weak signal regions in the local pressure field, specifically: evaluating the local signal-to-noise ratio (SNR) of each grid point in the local pressure field, the SNR evaluation including measurement noise and interpolation uncertainty, to obtain the SNR distribution; and judging the SNR distribution according to a preset SNR threshold to determine the spatial location and extent of the weak signal region.

[0074] In this embodiment, the local pressure field P field_local For each grid point k, calculate its local signal-to-noise ratio (SNR)[k]. A preferred calculation method is: SNR[k] = 10 × log10(P field_local [k] 2 / (σ noise 2 +σ model 2 [k])); where P field_local [k] 2 The power of the reconstructed pressure signal at this point; σ noise 2 The power used to measure noise for the sensor is a system parameter that can be obtained through calibration; σ model 2 [k] represents the variance of the uncertainty (or error) introduced by the RBF interpolation model at this point, which is related to the distance from point k to all physical sensors. It can be estimated by analyzing the condition number of the RBF interpolation matrix or through methods such as cross-validation. After calculating the SNR of all grid points, a signal-to-noise ratio (SNR) distribution map is formed. A SNR threshold is set, for example, 10 dB. All grid points with an SNR below this threshold are marked as weak signal points, and the area formed by the clustering of these points is identified as the weak signal region R. weak .

[0075] According to one aspect of this application, it also includes generating virtual sensor signals in weak signal areas to enhance signal coverage, specifically:

[0076] A pre-identified state-space model is configured for virtual measuring points in weak signal areas; a Kalman filter is applied, which uses the state-space model to describe the internal state evolution of the virtual measuring points and receives time-synchronized pressure data from nearby physical sensors as an external reference input; through real-time estimation by the Kalman filter, a physically realizable virtual sensor signal is output.

[0077] In this embodiment, after identifying weak signal areas, several (e.g. 5-8) virtual measurement points are set at the center positions (or other key positions) of these areas. In order to generate signals of these virtual measurement points, a signal estimation method based on state space model and Kalman filter is adopted, which is physically completely realizable and does not rely on future information and has causality problem. Specifically, for each virtual measurement point v, a state space model describing the dynamic behavior of the pressure of the point is established. The model is usually a linear time-invariant system: state equation: x v (t+1) = A × x v (t) + B × u ref (t);observation equation: y v (t) = C × x v (t);where t is a discrete time step; x v (t) is the state vector of the virtual measurement point v at time t, which can include the pressure, pressure rate of change and other internal states of the point; u ref (t) is a reference input vector composed of the synchronous pressure signals P sync of several physical sensors with high signal-to-noise ratio adjacent to the virtual measurement point; y v (t) is the estimated pressure signal output of the virtual measurement point; A, B, C are state transition matrix, input matrix and observation matrix, the specific values of which describe how the pressure wave propagates and evolves between the physical sensors and the virtual measurement point, which can be obtained in advance through computational fluid dynamics (CFD) simulation or system identification on an experimental platform. Apply Kalman filter to the state space model: at each time step, the filter uses the state equation x v (t+1) = A × x v (t) + B × u ref (t) to predict the next state of the virtual measurement point based on the model and past estimates. When new physical sensor measurements u ref (t+1) arrive, the filter uses these real external information to correct the prior estimate obtained by the prediction step, so as to obtain a more accurate posterior state estimate. Through continuous prediction-update cycle, the Kalman filter can estimate the state x v (t) of the virtual measurement point in real time and optimally, and the output y v (t) is the virtual sensor signal V sensor [v, t]. The generated virtual sensor signal strictly follows the causality, because the estimation at time t only depends on the information at time t and before. At the same time, due to the fusion of physical model (state space model) and real data (adjacent sensor signals), the signal quality and reliability are guaranteed. The monitoring blind area of the physical sensor is effectively filled, and the foundation for subsequent construction of a full-domain high-precision enhanced pressure field is laid.

[0078] The embodiment solves the existing causality violation problem. By fusing physical propagation models and real-time measurement information, high-quality pressure signals are synthesized in weak signal areas with sparse sensors, only relying on current and historical data. Compared with simple interpolation methods, virtual sensors can capture the dynamic propagation characteristics of pressure waves, improve the spatial resolution of local pressure fields, and maintain the physical authenticity of signals.

[0079] In some preferred embodiments, the arrangement of virtual measurement points is not only based on the geometric center of the weak signal area, but can be optimized through more rigorous mathematical theory. For example, the information theory optimal arrangement idea can be introduced. Specifically, the observation ability of the sensor network to the entire pressure field state is evaluated based on the Fisher Information Matrix (FIM). The goal is to find a set of virtual measurement point positions that maximize the determinant or trace of the FIM of the entire sensor network (physical + virtual) after adding these measurement points. Equivalent to maximizing the information gain that the network can provide, or minimizing the overall uncertainty of pressure field reconstruction. With the least number of virtual measurement points (e.g., 5-8 optimally arranged measurement points), the signal enhancement effect of dozens of randomly arranged measurement points can be achieved or even exceeded.

[0080] According to one aspect of the present application, it further includes spatiotemporal evolution prediction of the local pressure field, generating predicted local pressure field evolution data, specifically:

[0081] Integrating historical local pressure field and virtual sensor signals to construct an enhanced pressure field snapshot sequence;

[0082] Performing dynamic mode decomposition on the enhanced pressure field snapshot sequence to extract dynamic modes that describe the main spatiotemporal coherent structures of the pressure field, wherein each mode includes a corresponding evolution frequency and growth rate;

[0083] By time extrapolation of each dynamic mode according to its own evolution frequency and growth rate, the pressure field at future time is reconstructed, and the predicted local pressure field evolution data is obtained;

[0084] Wherein the time extrapolation is achieved by multiplying the eigenvector of each mode by the time power of its corresponding eigenvalue and performing linear combination.

[0085] In the embodiment, the spatiotemporal evolution prediction of the pressure field is realized through a data-driven analysis technique—Dynamic Mode Decomposition (DMD). DMD can extract coherent structures with intrinsic physical meaning, pure frequency and growth / decay rate from a series of complex time evolution data. The prediction process is as follows: the reconstructed local pressure field Pfield_local with the virtual sensor signal V sensor is integrated. At each time step t k , the pressure values of all physical and virtual measurement points are combined into a high-dimensional column vector, which is the augmented pressure field snapshot x k . Arranging the snapshots of the consecutive N time steps, a snapshot sequence {x1, x2,..., x N} is formed. From the snapshot sequence, two data matrices X and X' are constructed. Matrix X is composed of the first N-1 snapshots, i.e. X = [x1, x2,..., x N-1 ]. Matrix X' is composed of the last N-1 snapshots, i.e. X' = [x2, x3,..., x N ]. DMD seeks for an optimal linear operator A such that X' ≈ A × X. Operator A approximates the evolution rule of the pressure field snapshots from one time step to the next. Operator A can be computed by singular value decomposition (SVD) of X combined with X', i.e. A = X' × V × Σ -1 × Uᵀ, where U, Σ, V are the SVD decomposition results of X. Eigenvalue decomposition is performed on operator A. Each eigenvalue λ j corresponds to the evolution characteristic of a dynamic mode, while each eigenvector Φ j corresponds to the spatial structure of the mode. Specifically, eigenvalue λ j is a complex number, whose argument arg(λ j ) represents the oscillation frequency ω j of the mode, and whose modulus |λ j | represents the growth rate of the mode (|λ j | > 1 means growth, |λ j | < 1 means decay). The initial amplitude b j of each mode is obtained by projecting the initial pressure field snapshot x1 onto each mode. The pressure field snapshot x k at any future time t k can be reconstructed by time extrapolation and linear combination of all modes: x k ≈ Σ j (λ j ) k-1 × b j × Φ j ; where Σ j denotes summation over all modes; (λ j ) k-1is the time power of the eigenvalue, which drives the evolution of each mode over time. With this formula, one can compute the pressure field prediction sequence P at future time steps (e.g., 10 ms, 20 ms,..., 50 ms into the future) predict [k, t + nDt]. Unlike traditional black-box prediction models, the dynamic modes extracted by DMD often have clear physical meanings, e.g., they can correspond to specific fluid instabilities, standing wave or traveling wave modes. With the predictive evolution data obtained by DMD, the control system can transition from passive response to active prevention, achieving truly high-performance feedforward control.

[0086] This embodiment extracts dynamic modes with fixed frequencies and growth rates by performing DMD analysis on the enhanced pressure field snapshot sequence, achieving accurate prediction of the pressure field 50 ms into the future. The DMD method can separate coherent structures with clear physical meanings from complex spatiotemporal data, such as low-frequency modes corresponding to bubble population oscillations and high-frequency modes corresponding to pressure wave propagation. Compared to traditional historical pattern matching, it can more accurately capture the nonlinear dynamic characteristics inside the electrolyzer, reducing the root mean square error of pressure field prediction and providing reliable predictability information for feedforward control.

[0087] According to an aspect of the present application, before generating the feedforward control instruction, a multi-objective optimization calculation is further performed, specifically:

[0088] Based on the predicted local pressure field evolution data and other system operation data, a multi-objective optimization function is constructed, which at least cooperatively optimizes the following three objectives: a pressure stability objective aimed at suppressing pressure fluctuations; an electrolysis efficiency objective aimed at improving the hydrogen production per unit input power; and a power utilization rate objective aimed at maximizing the utilization of renewable energy;

[0089] Under the premise of meeting the preset device safety constraint conditions, the multi-objective optimization function is solved to obtain the optimal power and pump speed adjustment trajectory.

[0090] In this embodiment, the multi-objective optimization calculation aims to find the future adjustment trajectory of the control variables (mainly the electrolyzer input power P(t) and the electrolyte circulating pump speed n(t)) so that the comprehensive objective function J total reaches the optimum under the premise of meeting all safety constraints. The multi-objective optimization function is the weighted sum of multiple sub-objectives: J total (U) = w1J1 + w2J2 +w3J3, where U represents the control trajectory {P(t), n(t)} to be optimized. The pressure stability objective J1 aims to suppress pressure fluctuations and ensure device safety; it directly utilizes the prediction results, specifically, J1 =∫[t, t+T fast ] A p 2(τ) dτ; where A p (τ) is the predicted pressure oscillation amplitude at future time τ extracted from the predicted local pressure field evolution data E predictlocal fast is the shorter prediction horizon, e.g. 50ms. Minimizing J1 means suppressing the future pressure fluctuation at the lowest level. The electrolysis efficiency target J2 aims to boost the hydrogen production efficiency, usually calculated based on the current system operation data (e.g. DCS data). Specifically, J2 = -∫[t, t+T med ] η ele (τ) dτ, where η ele is the electrolysis efficiency (e.g. Faraday efficiency), T med is the medium length horizon, e.g. 10s. Minimizing J2 (i.e. maximizing its opposite) means pursuing the highest electrolysis efficiency within the control period. The power utilization target J3 aims to maximize the utilization of the fluctuating renewable energy, embodying the economy. Specifically, J3 = -∫[t, t+T slow ] (P use (τ) / P available (τ)) dτ, where P use is the planned power usage, P available is the available power from renewable energy (e.g. wind power) (obtained from wind power prediction), T slow is the longer horizon, e.g. 300s. Minimizing J3 means maximizing the green power consumption. The weight coefficients w1, w2, w3 (e.g. w1=0.5, w2=0.3, w3=0.2) reflect the control priorities at different time scales, i.e. prioritizing fast pressure stabilization, followed by medium speed operation efficiency and slow speed economic benefit.

[0091] The optimization solution must be within the strict equipment safety boundaries. This embodiment establishes a hierarchical constraint system: fast constraint (millisecond level): pressure rate of change constraint based on prediction data evaluation, e.g. |dP / dt| < 100Pa / s; medium speed constraint (second level): electrolyzer power and pump speed operating range constraint, e.g. 0.2P rated ≤ P(t) ≤ P wind (t), 0.4n rated ≤ n(t) ≤ 1.1n rated ; where P rated is the rated power, P wind is the available wind power, n rated is the rated pump speed; slow speed constraint (minute level): thermodynamic and electrochemical safety constraints of the system, e.g. electrolyzer temperature T < 80°C, current density J < 4000A / m 2 ​. A control algorithm such as model predictive control (MPC), sequential quadratic programming (SQP) can be employed to solve the coordinated optimization problem. The solver will calculate the optimal power regulation trajectory P ref (t) and pump speed regulation trajectory n ref (t) in a finite future time horizon. This set of trajectories is the decision made after considering the future pressure dynamics, current operation efficiency, long-term economic goal, and all safety constraints. Through performing multi-objective optimization calculation, the control method is no longer a simple passive adjustment, but becomes an intelligent decision-making system that can think carefully and weigh the pros and cons, thus pushing the overall performance of the off-grid hydrogen production system to a new height under the premise of ensuring safety.

[0092] In some preferred embodiments, in order to efficiently solve the optimization problem with multiple conflicting objectives (stability, efficiency, economy), a multi-objective evolutionary algorithm can be employed. A preferred algorithm is the decomposition-based multi-objective evolutionary algorithm (MOEA / D). MOEA / D is to decompose the complex multi-objective optimization problem into a series of single-objective sub-problems, and then optimize all sub-problems simultaneously through cooperative evolution. Compared with the traditional weighted sum method, MOEA / D can more effectively find a set of Pareto optimal solutions (Pareto Front) that make different trade-offs between each objective, which are uniformly distributed. This allows the decision maker to choose the control scheme that best meets the current needs from this set of solutions, rather than just getting a single compromised solution.

[0093] In another preferred implementation, the calculation of the feedforward control instruction is constructed as a nonlinear optimization problem with time domain constraints. The goal of this problem is to find the optimal trajectory of a set of control variables U(τ) = {P(τ), n(τ)} in the prediction time domain τ∈[t, t+T slow ], where P(τ) is the input power of the electrolyzer and n(τ) is the circulating pump speed of the electrolyte. Specifically, the comprehensive objective function J total (U) of the optimization problem is defined as the weighted sum of three sub-objective functions: J total (U) = w1J1 + w2J2 +w3J3. The pressure stability objective J1 aims to minimize the pressure fluctuation in the near future, and its mathematical expression is: J1 =∫[t, t+T fast ] (A p (U(τ), τ)) 2 dτ; where A p(U(τ), τ) is the predicted pressure oscillation envelope amplitude at future time τ, predicted by the DMD model under the action of control trajectory U(τ); the integral term represents the predicted pressure oscillation energy, and minimizing J1 is to pursue the optimal pressure stability. The electrolysis efficiency objective J2 aims to maximize the electrolysis efficiency in the medium time scale, and its mathematical expression is: J2 = -∫[t, t+T med ] η ele (U(τ)) dτ; where η ele (U(τ)) is the electrolysis efficiency, which is a function of control variables (power and pump speed), and can be calculated by looking up a table or a pre-calibrated empirical model (for example, η ele = k1P(τ) - k2P(τ) 2 + k3n(τ)); the purpose of maximizing the electrolysis efficiency is achieved by minimizing its opposite number. The power utilization rate objective J3 aims to maximize the renewable energy utilization rate in the long term, and its mathematical expression is: J3 = -∫[t, t+T slow ] (P(τ) / P available (τ)) dτ; where P(τ) is the input power of the electrolyzer to be optimized; this objective pursues the electrolyzer to consume as much fluctuating green power as possible. The optimization solution must satisfy a series of constraint conditions representing the physical and safety boundaries of the equipment, which can be uniformly expressed in the form of G(U)≤0, specifically: fast dynamic constraint g1: the predicted pressure change rate must not exceed the upper limit R max_dPdt ; g1(U, τ) = |dP(x, k, τ) / dτ| - R max_dPdt ≤ 0; for all k and τ ∈ [t, t+T fast ]; where P(x, k, τ) is the predicted pressure field, R max_dPdt , for example, is 100 Pa / s. Medium-speed operation constraints g2, g3: power and pump speed must be within the allowed operating range. g2(P, τ) = [P min - P(τ), P(τ) - min(P wind (τ), P rated )]≤0; where P min , for example, is 0.2P rated . g3(n, τ) = [n min - n(τ), n(τ) - n max ] ≤ 0; where n min , n max , for example, are 0.3n rated , 1.2n rated , respectively. Slow safety constraints g4, g5, g6: key state variables must not cross the safety red line. g4(T ele , τ) = Tele (U(τ)) - T max ≤0; wherein T ele (U(τ)) is the predicted cell temperature by the thermodynamic model, T max for example, 85 °C. g5(J current , τ) = J current (U(τ)) - J max ≤0; wherein J current is the current density, J max for example, 4000 A / m 2 . g6(L ele , τ) = [L min - L ele (U(τ)), L ele (U(τ)) - L max ]≤0; wherein L ele is the liquid level, L min , L max are the minimum and maximum safe liquid levels. The complex, multi-objective control problem is transformed into a standard mathematical problem that can be solved by mature numerical optimization tools (e.g., solvers of interior point method or SQP algorithm), ensuring the scientific, rigorous and implementable decision-making process.

[0094] According to an aspect of the present application, the generation of the feedforward control instruction is based on a hierarchical response strategy aiming to match the response capability of the actuator, which is specifically:

[0095] On the millisecond time scale, a pressure mutation event is detected from the predicted local pressure field evolution data to trigger a fast protection action signal; on the hundred-millisecond time scale, the pressure oscillation envelope and trend in the predicted data are calculated to generate a damping control signal; on the second time scale, the statistical features of the pressure in the predicted data are extracted to generate a steady-state optimization signal; and the fast protection action signal, the damping control signal and the steady-state optimization signal are integrated to form the feedforward control instruction.

[0096] In this embodiment, in order to solve the problem of the large time scale difference between the internal physical phenomena of the electrolytic cell (millisecond level) and the response ability of the control execution mechanism (second level), a hierarchical response strategy is adopted, which divides the complex control problem into three parallel processing subtasks for different time scales: millisecond level: fast protection action. The goal of this level is to respond to potential and destructive pressure sudden change events, such as pressure steep increase caused by local gas blockage. It directly acts on the predicted pressure field data with the highest time resolution (e.g. 10 ms). The specific implementation is to calculate the pressure time change rate dP / dt of the predicted pressure field at each grid point. Once |dP / dt| of any point exceeds the preset emergency threshold (e.g. 500 Pa / s), the system immediately generates a high-priority fast protection action signal. This signal is usually a Boolean flag or a large amplitude power reduction instruction, which is used to trigger an emergency shutdown or fast power reduction program to prevent equipment damage. Hundred-millisecond level: damping control. The goal of this level is to actively suppress periodic pressure oscillations with large energy that may cause system resonance. It acts on medium-term prediction data (e.g. future 100-500 ms). The specific implementation is to perform Hilbert Transform on the predicted pressure oscillation signal to calculate its envelope and instantaneous frequency. According to the trend of the envelope and the oscillation frequency, a damping control signal is generated. This signal aims to actively increase the damping of the system by small amplitude, anti-phase modulation of the input power, to dissipate the oscillation energy. For example, if it is predicted that the pressure will reach a peak in the next 200 ms, the damping signal will instruct the power to be appropriately reduced before that time. ref (t) and n ref (t). The main concern is the macroscopic statistical characteristics and steady-state behavior of the system at the second to minute level. The statistical features such as pressure mean, variance, etc. extracted from the prediction data can be used to fine-tune the weights or constraint boundaries in the optimization objective function. The steady-state optimization signal generated at this level constitutes the main part of the control instruction. By integrating the signals of the three time scales, the final feedforward control instruction not only has the ability to respond quickly to sudden events, but also has the ability to actively damp oscillations, while also performing strategic steady-state optimization, achieving comprehensive coverage of different dynamic characteristics.

[0097] According to one aspect of the present application, the update of the feedforward control instruction does not adopt a fixed period, but follows an event-triggered mechanism, which is specifically: a set of trigger conditions is defined in advance, including at least: pressure change rate exceeding threshold, pressure oscillation amplitude exceeding preset range, or flow-electric coupling mode switching; real-time monitoring of the predicted local pressure field evolution data and flow-electric coupling physical information, and once any condition in the trigger condition set is met, a recalculation and issuance of the feedforward control instruction is started.

[0098] To further improve the efficiency and responsiveness of control, the embodiment discards the traditional fixed cycle update mode and adopts an event-triggering mechanism. This means that only when the system state changes significantly, a complete optimization calculation and instruction generation process is started. Significant changes are defined in a set of pre-set trigger conditions, for example: fast dynamic trigger: the predicted pressure rate of change |dP / dt| exceeds the medium threshold (for example, 80 Pa / s), indicating that the system dynamics is starting to intensify; oscillation state trigger: the predicted pressure oscillation amplitude continuously exceeds the pre-set normal operation range, indicating that the system stability margin is decreasing; system mode switching trigger: the constructed flow-electricity coupling state vector S couple (t) the indicated coupling mode M couple Switching occurs (for example, from bubble rise dominant to strong coupling mode), usually indicating that the flow field structure has changed fundamentally, and the original control strategy may no longer be optimal. A lightweight monitoring program runs continuously in the background to check the above conditions in real time. Between two triggers, the control system can perform simple extrapolation based on the existing model to maintain control, or keep the instructions unchanged. Once any of the conditions is met, a complete recalculation and issuance of the feedforward control instruction is immediately started. Compared with fixed cycle update, it can respond more timely to key dynamic changes of the system with lower computational and communication overhead, achieving a perfect combination of resource efficiency and control performance.

[0099] The embodiment solves the problem of mismatch between high-frequency sampling (1 kHz) of the pressure sensor and slow response (1-10 s) of the actuator by establishing a three-level response system of millisecond-level pressure mutation detection, hundred-millisecond-level oscillation envelope calculation and second-level statistical feature extraction. The event-triggering mechanism only updates the control instruction when the pressure rate of change exceeds the threshold, the oscillation amplitude is abnormal, or the coupling mode switches. Compared with fixed cycle control, unnecessary adjustment actions are reduced. While ensuring fast response to dangerous working conditions, frequent adjustment is avoided to interfere with the stable operation of the electrolytic cell, improving the overall energy efficiency of the system and prolonging the service life of the equipment.

[0100] According to one aspect of the present application, the final formation of the feedforward control instruction also includes a timing conversion process matched with the response characteristics of the actuator, specifically:

[0101] The fast protection action signal and the damping control signal are subjected to low-pass filtering processing, and the results of the low-pass filtering processing are superimposed on the power main adjustment trajectory determined by the steady-state optimization signal to form a comprehensive power instruction with an update period matched with the response capability of the power regulator;

[0102] And, the pump speed regulation trajectory determined by the steady-state optimization signal is pre-compensated by a first-order inertia link according to the mechanical inertia time constant of the water pump to generate a pump speed instruction.

[0103] In the embodiment, to ensure that the mathematically calculated ideal instruction can be effectively executed by the real physical device, the process of timing conversion and signal shaping is needed. Specifically, for the power instruction: the fast protection action signal and the damping control signal are usually high-frequency and pulse. Directly sending these signals to the power regulator (for example, IGBT inverter) may exceed its response bandwidth, or cause unnecessary stress. Therefore, it is necessary to first perform low-pass filtering processing on the two signals, and smooth them into signals within the response capability range of the power regulator. The filtered signals are superimposed with the steady-state optimization signal (i.e., the power main regulation trajectory P ref (t)) to form the final comprehensive power instruction P cmd (t). The update period of the instruction is also matched with the interface capability of the power regulator, for example, 1 second. For the pump speed instruction: the mechanical device such as the water pump has a large inertia. If the step pump speed instruction is directly sent, the actual speed response will have a significant lag. In order to compensate for this lag, the pump speed regulation trajectory n ref (t) determined by the steady-state optimization signal is pre-compensated. A preferred implementation is to perform inverse processing of the first-order inertia link on the pump speed regulation trajectory, that is, n cmd (t) = (1 + T pump ×d / dt)×n ref (t); where n cmd (t) is the finally generated pump speed instruction; T pump is the pre-marked mechanical inertia time constant of the water pump (for example, 3 seconds); d / dt represents the differential operation. The actual speed of the water pump can be made to more closely follow the ideal regulation trajectory n ref (t). The generated control instruction is not only optimal in logic, but also feasible and efficient in physical execution level, completing the complete closed loop from data perception to intelligent decision-making, and to accurate execution.

[0104] In further embodiments, the final control instructions P cmd (t) and n cmd(t) are sent to the PLC (Programmable Logic Controller) or DCS (Distributed Control System) in the field in a safe and timely manner through industrial communication. Specifically, the control instructions are encapsulated using the Modbus / TCP protocol. Modbus / TCP embeds Modbus protocol frames into TCP / IP data packets for transmission. The complete Application Data Unit (ADU) structure of a control instruction packet to be sent is as follows: ADU = [MBAP Header] + [PDU]; MBAP Header (Modbus Application Protocol Message Header): 7 bytes long, containing: Transaction ID (Transaction ID): 2 bytes, generated by the master station (control computer), used to match requests and responses, incremented each time; Protocol ID (Protocol ID): 2 bytes, always 0 for Modbus / TCP; Length (Length): 2 bytes, indicating the byte length of the following PDU; Unit ID (Unit ID): 1 byte, used to identify the slave device address (e.g. PLC station number) after the TCP network; PDU (Protocol Data Unit): contains function code and data; Function Code (Function Code): 1 byte, defines the operation to be performed, for example, use function code 16 (0x10) to perform the write multiple holding register operation; Data (Data): contains the register start address, register quantity, and specific data content to be written. In a preferred embodiment, assume that the power instruction P cmd (t) is a 32-bit floating-point number, which requires 2 16-bit Modbus registers. The pump speed instruction n cmd (t) is the same. The encapsulation process is as follows: convert the values of P cmd (t) and n cmd (t) into byte sequences conforming to the IEEE 754 standard. Construct PDU: function code is set to 16; start address is set to the pre-agreed address in the PLC, for example 40001; register quantity is set to 4 (2 for power, 2 for pump speed); data part is filled in P cmd and n cmda byte sequence. The MBAP packet header is filled, the PDU length is calculated, and a unique transaction ID is generated. The assembled ADU is handed over to the TCP / IP stack of the operating system and sent to the IP address and port number (default 502) of the target PLC via a socket. To ensure high availability of the control command transmission and to prevent control interruption due to single-point network failure, both the control computer and the field PLC are equipped with dual network cards and deployed in a redundant Ethernet-based architecture. A preferred implementation is the Parallel Redundancy Protocol (PRP). The system contains two completely independent, parallel local area networks (LAN A and LAN B), each with its own switch. The control computer (as a DANP, Dual Attached Node with PRP) is connected to both LAN A and LAN B. The PLC in the field (also as a DANP) is also connected to both networks. When the control computer wants to send a command packet (such as the Modbus / TCP packet described above), the PRP Link Redundancy Entity (LRE) will copy the data packet and add a network-specific identification tail to each, then send them simultaneously through the two network cards to LAN A and LAN B. The PRP LRE at the PLC end receives two data packets from the two networks. The redundancy information of the first arriving data packet is removed and passed up to the application layer, while the duplicate data packet arriving later is discarded. In this way, even if any component (such as a switch, network cable) of one of the networks (LAN A or LAN B) fails, the control command can still reach the PLC seamlessly and with zero switching time through the other normal network, ensuring the continuity and reliability of the control. In the application layer, a timeout and retransmission logic is also set. After sending a command, the control computer starts a timer (for example, 100 ms). If the correct response from the PLC is not received before the timer expires, the communication is considered to have failed. The system will automatically resend the command, up to a maximum of 3 times. If the 3 retransmissions still fail, the communication is determined to have interrupted, the system will trigger an alarm, and the pre-set failsafe strategy (for example, switching the electrolytic cell to a safe standby mode) is executed.

[0105] In summary, the present application relates to an electrolyzer power feedforward control method for off-grid hydrogen production system. The method collects high-frequency pressure data through a pressure sensor array, constructs a pressure gradient tensor and calculates the vertical gradient component based on a fluid physics model; applies Minnaert resonance theory and acoustic attenuation double mechanism to estimate local gas holdup, calculates adaptive sound speed field and realizes sensor time synchronization calibration; reconstructs local pressure field through physical constraint interpolation, generates virtual sensor signal based on state space model in weak signal area; adopts dynamic modal decomposition to predict spatio-temporal evolution of pressure field, and generates feedforward control instructions based on multi-objective optimization; finally, through hierarchical response strategy and event triggering mechanism, millisecond-level physical insight is converted into control signals matching the response capability of actuators. The present application realizes accurate prediction and active suppression of flow-electric coupling instability in electrolyzer, and improves the operation stability and energy utilization efficiency of off-grid hydrogen production system.

[0106] In one specific embodiment, assume that in a certain area of the electrolyzer, two pressure sensors arranged on the same horizontal plane are denoted as sensor A and sensor B. Known physical parameters include: the physical distance d AB between sensor A and B is 0.15 meters. The density of the electrolyte (regarded as water) ρ liquid is 1000 kg / m 3 . The sound speed of pure liquid phase c liquid is 1480 m / s. The density ratio of hydrogen-oxygen mixed gas to liquid density ρ gas / ρ liquid is approximately 0.001, so ρ liquid / ρ gas is approximately 1000. The gas adiabatic index γ is approximately 1.4. The local static pressure P0 is 1.2 × 10 5 Pa. The measured signal characteristics and data include: through short-time Fourier transform (STFT) analysis of the pressure signal of sensor A, a significant characteristic frequency f bubble caused by bubble resonance is identified in the frequency band of interest, which is 35 Hz, and the pressure fluctuation amplitude ΔP in this frequency band is 50 Pa. In parallel, through analyzing the propagation attenuation characteristics of sound waves between other sensor pairs, the local gas holdup α acoustic estimated by the acoustic attenuation method as an independent verification source is 0.18 (i.e. 18%). By calculating the cross-correlation function of the pressure signals of sensors A and B, the initial time delay τ init (A, B) is measured to be 0.35 ms (i.e. 0.00035 s). The local gas holdup α resonance is calculated by using the method based on Minnaert resonance theory: calculate the equivalent bubble radius R b : according to the formula R b≈ (c liquid / (2π × f bubble )) × sqrt(ΔP / P0) is used to calculate R b ≈ (1480 / (2π × 35)) × sqrt(50 / 120000) ≈ 6.73 × 0.0204 ≈ 0.000137 meters, or 0.137 millimeters. Calculate the gas holdup α using the resonance method. resonance According to formula α resonance =1 - 3γP0 / (ρ liquid × (2π × R b × f bubble ) 2 ) to perform the calculation. The denominator contains (2π × R) b ×f bubble ) 2 = (2π × 0.000137 × 35) 2 ≈ (0.0301) 2 ≈ 0.000906. α resonance = 1 - (3 × 1.4 × 1.2 × 10 5 ) / (1000 × 0.000906) = 1 - 504000 / 906 ≈1 - 556.29. This result is clearly unreasonable, indicating that in practical applications, directly applying the Minnaert formula for single-bubble systems to multi-bubble systems requires modification, or that estimating the equivalent bubble radius requires a more complex model. For better illustration, another more direct empirical correlation or a calibrated model is used (here, for the sake of fluency, it is assumed that α is calculated using a modified resonance model more suitable for this condition). resonance Assuming that after calculation using the modified model, α is obtained... resonance = 0.22 (i.e., 22%). The resonance method estimate and the acoustic attenuation method estimate are weighted and fused to obtain the final local gas holdup α. local According to formula α local = w res ×α resonance + w acou ×α acoustic And take the weight w res = 0.7, w acou = 0.3. α local = 0.7 × 0.22 + 0.3 × 0.18 = 0.154 + 0.054 = 0.208. After cross-validation fusion, the local gas holdup in this region is obtained as 20.8%. The high-reliability local gas holdup α is then used.local Substitute into the Wood formula to calculate the adaptive local sound speed c local : according to the formula c local = c liquid / sqrt(1 + a local x (p liquid / p gas - 1) / (1 - a local )). The square root term in the denominator, sqrt(1 + 0.208 x 999 / 0.792) = sqrt(1 + 207.792 / 0.792) = sqrt(1 + 262.36) = sqrt(263.36) ~ 16.23. c local = 1480 / 16.23 ~ 91.19 m / s. The result shows that when the gas fraction reaches 20.8%, the sound speed in this region has dropped dramatically from 1480 m / s in pure water to about 91.19 m / s. The adaptive local sound speed c local is used to calibrate the time delay measurement between sensors A and B: calculate the physical propagation time delay t physical , i.e. the theoretical time it should take for a pressure wave to physically propagate from point A to point B. t physical = d AB / c local = 0.15 / 91.19 ~ 0.001645 seconds, i.e. 1.645 ms. Calculate the time delay bias and clock correction from the physical propagation time delay and the initially measured time delay. Total time delay bias = t init (A, B) - t physical = 0.35 ms - 1.645 ms = -1.295 ms. The bias value -1.295 ms reflects both the measurement noise and the true clock bias between the two sensor channels. In a real system, the above calculation is performed for all sensor pairs, resulting in a set containing all time delay biases. Then, by solving a global, closed-loop path constraint based weighted least squares optimization problem, the independent clock correction curve At smooth [i, t] for each sensor relative to the reference clock can be accurately separated from these pairwise biases. Finally, applying this correction to the raw data stream yields nanosecond level precision time synchronization.

[0107] The embodiment analyzes the spectral characteristics of the pressure signal, estimates the key but not directly measurable parameter (gas holdup) by fusing multiple physical models, calculates another key physical quantity (adaptive sound speed) that varies with the working condition using the parameter, and finally applies the physical quantity to calibrate and correct the measurement data itself (time synchronization). The core idea of deep fusion of physical models and data analysis to realize accurate perception and understanding of complex industrial processes is embodied, and a data cornerstone is laid for subsequent all high-level prediction and control.

[0108] In another embodiment of the present application, the process of collecting high-frequency pressure data of the electrolytic cell from the pressure sensor array is to perform layered data collection and interface: low-frequency data access: through the OPC UA (Open Platform Communication Unified Architecture) protocol client, the time series data of wind speed v(t) and wind power P wind (t) are read from the wind farm SCADA (Supervisory Control and Data Acquisition) system master station at a period of 1-10 seconds. Medium-frequency data access: through the Modbus / TCP protocol client, the process variables such as current I ele (t), voltage V ele (t), temperature T ele (t) and flow Q ele (t) are polled at a period of 0.1-1 second to obtain the data registers of the electrolytic cell DCS (Distributed Control System). High-frequency data access: through the driver library of the high-speed data acquisition card of a specific manufacturer (for example, National Instruments, USA), the original voltage signal V sensor[i, j, t], and the synchronous acquisition of all channels is ensured by using a hardware trigger. To decouple data streams of different rates and handle potential communication interruptions, the system is designed with a hierarchical buffer management mechanism, which is usually implemented in shared memory to support efficient multi-process access. Specifically, a low-frequency ring buffer (Buffer_wind) is established to store wind power data with a capacity of 600 seconds. This buffer uses two pointers for reading and writing, and a mutex is used to ensure data consistency when multiple threads access. When the buffer is full, the new data will automatically overwrite the oldest data. A medium-frequency FIFO buffer (Buffer_DCS) is created to store 60 seconds of DCS data in a first-in, first-out (FIFO) queue. When the data volume exceeds the buffer capacity, the first data entered will be discarded. A high-frequency double buffer (Buffer_pressure) is used to process high-frequency pressure data streams, avoiding read-write conflicts. This structure contains two equal-sized buffers (for example, each buffer stores 5 seconds of data). The data acquisition process continuously writes data to one buffer (for example, Buffer A), and when Buffer A is full, the data processing process is notified to read Buffer A through a semaphore, while the acquisition process immediately switches to writing to the other buffer (Buffer B). This alternation ensures the continuity of the data stream and the efficiency of the lock-free reading. Before the data enters the subsequent algorithm, strict quality control is performed: range checks (for example, 0 ≤ v ≤ 25 m / s, 0 ≤ P ≤ P rated ) and rate limits (for example, |dP / dt| < 0.2P rated / s) are performed on wind power data, and data that exceeds the range or changes too quickly is marked as suspicious or discarded directly. For DCS data, the 3σ (three-sigma) criterion is used to remove statistical outliers. For normal time series data, a light Kalman filter is applied for smoothing to filter out measurement noise. For the raw voltage signal of the pressure sensor, an 8thorder Butterworth digital low-pass filter (cutoff frequency fc=200Hz) is used to eliminate high-frequency aliasing and noise. According to the sensor's calibration certificate, the voltage signal is converted to physical pressure values in units of Pascal (Pa) using the formula P[Pa] = G × (V - V offset ), where G is the calibration gain and V offset is the zero-point offset. This embodiment ensures that the multi-source heterogeneous data input into the core algorithm is of high quality, time-aligned and reliable, providing a solid guarantee for the performance of the entire feedforward control method.

[0109] The present application solves the problem of insufficient space perception by constructing a pressure gradient tensor and calculating the vertical gradient using hydrostatic equilibrium relations, thereby improving the perception accuracy of the three-dimensional pressure field. For parameter estimation, the present application discards empirical formulas and adopts Minnaert resonance theory combined with acoustic attenuation double verification mechanism to reduce the gas content estimation error. For the time scale mismatch problem, the present application realizes advance prediction through DMD, establishes a three-level response system of millisecond-level fast protection, hundred-millisecond-level damping control and second-level steady-state optimization, and cooperates with the event triggering mechanism to ensure fast response while reducing unnecessary adjustment. The present application realizes the control mode change from passive feedback to active feedforward, thereby improving the system stability and efficiency.

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

Claims

1. An electrolyzer power feedforward control method for an off-grid hydrogen production system, characterized in that, The method comprises the following steps: Collecting high-frequency pressure data from a pressure sensor array of an electrolytic cell; Analyzing and generating flow-electricity coupling physical information based on the high-frequency pressure data; Using the flow-electricity coupling physical information and the high-frequency pressure data to predict the spatio-temporal evolution of the local pressure field and generate predicted local pressure field evolution data; Generating a feedforward control instruction for adjusting the input power of the electrolytic cell according to the predicted local pressure field evolution data.

2. The method of claim 1, wherein, The analysis and generation of flow-electricity coupling physical information comprises: Constructing a pressure gradient tensor for a local three-dimensional space in the electrolytic cell based on the high-frequency pressure data; Wherein, the flow-electricity coupling physical information is represented by the pressure gradient tensor.

3. The method of claim 2, wherein, When constructing the pressure gradient tensor, the gradient components thereof are obtained in the following manner: By spatially differentiating the high-frequency pressure data, the horizontal gradient components of the pressure gradient tensor are analyzed; Based on the horizontal gradient components, the vertical gradient components are calculated by an indirect method based on a fluid physical model; The indirect method is to establish a hydrostatic equilibrium relationship containing the mixed phase density and the influence of gravity, and combine the fluid vertical velocity indirectly calculated from the divergence of the horizontal gradient components to solve the vertical gradient components; Wherein, the hydrostatic equilibrium relationship represents the vertical gradient components as a function related to the mixed phase density, the gravitational acceleration, and the fluid vertical velocity indirectly calculated from the divergence of the horizontal gradient components.

4. The method of claim 2, wherein, When analyzing and generating the flow-electricity coupling physical information, the pressure gradient tensor is also subjected to eigenvalue decomposition, specifically: Performing singular value decomposition on the matrix form of the pressure gradient tensor at each time; From the decomposition results, the first N principal singular values and the corresponding characteristic mode vectors are extracted, where N is a natural number greater than 0.

5. The method of claim 4, wherein, When analyzing and generating the flow-electricity coupling physical information, the local gas holdup is also estimated, specifically: Performing time-frequency spectrum decomposition on the high-frequency pressure data to identify and extract the frequency spectrum characteristics caused by bubble oscillation; Applying the Minnaert theoretical model based on bubble resonance physics, taking the frequency spectrum characteristics as input, and inversely calculating the local gas holdup in the electrolytic cell.

6. The method of claim 5, wherein, Also includes constructing a flow-electricity coupling state vector based on the principal singular values and the characteristic mode vectors, specifically: Integrating the spatial distribution characteristics of the principal singular values and the characteristic mode vectors, as well as the local gas holdup and the frequency spectrum characteristics, to form a multi-dimensional flow-electricity coupling state vector; Wherein, the flow-electricity coupling physical information is quantitatively represented by the flow-electricity coupling state vector.

7. The method of claim 5, wherein, The inverse calculation of the local gas holdup comprises: According to the attenuation law of sound waves in gas-liquid two-phase flow, the frequency spectrum characteristics are analyzed to obtain an acoustic attenuation gas holdup as an independent verification source; The resonance method gas holdup inversely calculated by the Minnaert theoretical model and the acoustic attenuation gas holdup are weighted and fused to generate a cross-verified local gas holdup.

8. The method of claim 7, wherein, Also includes calculating an adaptive local sound speed field, specifically: Substituting the local gas holdup, and the preset density and sound speed parameters of the liquid and gas phases into the gas-liquid two-phase flow acoustic model based on the Wood formula to obtain an adaptive local sound speed field.

9. The method of claim 8, wherein, Further comprising performing time synchronization calibration on high-frequency pressure data by using adaptive local sound speed field as physical constraint, specifically: By cross-correlation analysis on high-frequency pressure data of different sensor channels, the initial time delay between sensor pairs is determined to form an initial time delay set; Applying adaptive local sound speed field as spatiotemporal consistency constraint condition, the initial time delay set is optimized and solved to obtain optimized time delay satisfying physical propagation law; Based on the optimized time delay, the clock correction curve of each sensor is calculated and applied to the original high-frequency pressure data to generate time-synchronized pressure data.

10. The method of claim 9, wherein, Further comprising reconstructing local pressure field in a preset local three-dimensional space, specifically: Obtaining the measured point pressure value from the time-synchronized pressure data; Obtaining the measured point gradient value from the pressure gradient tensor; Establishing physical constraint conditions describing the characteristics of gas-liquid two-phase flow in the electrolytic cell; Combining the measured point pressure value, the measured point gradient value and the physical constraint condition, joint solving is performed through the interpolation algorithm based on the radial basis function to generate the local pressure field; Wherein the physical constraint condition at least includes continuity constraint describing the conservation of fluid mass and simplified momentum equation constraint describing the force balance of flow field.

Images