Micro-grid hydrogen power supply power accurate adjustment and overshoot suppression method

Abstract

The invention discloses a micro-grid hydrogen power supply power accurate adjustment and overshoot suppression method, which comprises the following steps: acquiring sound waves generated in an electrochemical process in a proton exchange membrane in a hydrogen power supply to obtain an acoustic response signal set; based on the acoustic response signal set, reconstructing the current density in the proton exchange membrane to obtain a spatial distribution diagram with high temporal-spatial resolution; generating a directional acoustic compensation signal according to the non-uniformity represented by the spatial distribution diagram, and applying the directional acoustic compensation signal to a target area of the proton exchange membrane; and synthesizing the data, and generating a final control instruction for adjusting the hydrogen power supply through a multi-time scale coordination control architecture. According to the invention, real-time perspective and targeted physical intervention of the internal state of the fuel cell are realized through a non-intrusive acoustic means, and the dynamic response speed, the power regulation precision and the operation stability of the hydrogen power supply are remarkably improved.

Description

TECHNICAL FIELD

[0001] The present application relates to power electronics, in particular, a method for precise regulation and overshoot suppression of hydrogen power in a microgrid. BACKGROUND

[0002] Under the background of building a new power system dominated by renewable energy, microgrids, as a key technology for efficient utilization of distributed energy and improving power supply reliability, have attracted much attention. Hydrogen energy, due to its high energy density and clean and pollution-free characteristics, is considered one of the most potential energy carriers. Proton exchange membrane fuel cells (PEMFC) as the core device for efficient conversion of hydrogen energy play an increasingly important role in microgrids, especially in terms of stabilizing renewable energy fluctuations, providing backup power and black start support, which shows great application value. Therefore, studying how to achieve precise and rapid regulation of hydrogen power in microgrids to ensure stable and reliable grid-connected operation under severe load changes is of great theoretical and engineering significance to improve the overall dynamic performance and power quality of microgrids.

[0003] Currently, research on hydrogen power control has made some progress. Existing technical solutions mainly focus on optimizing control algorithms. For example, on the basis of traditional proportional-integral-derivative (PID) control, intelligent algorithms such as fuzzy logic and neural networks are introduced to improve the controller's adaptability to system nonlinear characteristics. In addition, to improve the response speed of the system, some solutions use feedforward control strategies to predict load changes and adjust the supply of reactants (such as hydrogen and air flow) in advance. In terms of diagnosis, electrochemical impedance spectroscopy (EIS) is a commonly used analysis tool that can obtain average information about the internal electrochemical processes of the battery (such as ohmic resistance, charge transfer resistance, and mass transfer resistance) by applying a small-amplitude sinusoidal disturbance and measuring the impedance response of the system, providing a basis for system state assessment and fault diagnosis. These methods have improved the steady-state control accuracy and dynamic response characteristics of hydrogen power to some extent.

[0004] However, the prior art still faces deep technical bottlenecks in dealing with the rapid transient process of hydrogen power supply, especially the power overshoot problem under large power step response. These problems are not caused by the complexity of the control algorithm itself, but by the lack of perception of the internal physical state of the hydrogen power supply and the non-targeted control means. Specifically, the prior art mainly has the following deep-seated problems: First, there is a lack of real-time, high-resolution diagnostic capability for the physical root cause of power overshoot. The essence of power overshoot is that the current density inside the proton exchange membrane (PEM) produces a severe spatial non-uniformity on the millisecond time scale at the moment of power step-up. However, existing diagnostic tools such as electrochemical impedance spectroscopy (EIS) have a response speed of seconds and can only provide averaged information over the entire electrode area, making it completely invisible to such transient, localized internal state evolution processes, and making it impossible for the controller to obtain the root cause of the problem. Second, there is a spatial mismatch between control execution and physical problems. Existing control methods, whether adjusting the total flow of reactants or adjusting the output of the DC / DC converter, the object of action is the entire fuel cell, which is a globalized control input. However, it is trying to solve a localized non-uniformity problem. This kind of flood irrigation control method cannot differentiate between cold spots (under-reaction areas) and hot spots (over-reaction areas) in the membrane, and its inherent hysteresis and non-targeted nature can easily exacerbate system oscillation and overshoot. SUMMARY

[0005] The purpose of the invention is to solve the technical problems of severe power overshoot and slow dynamic response caused by the inability of the prior art to diagnose and regulate the spatial non-uniformity of the internal electrochemical process of the proton exchange membrane in real time.

[0006] The technical solution is a method for precise power regulation and overshoot suppression of a microgrid hydrogen power supply, comprising: Acquire acoustic response signal sets generated by the acoustic waves produced by the electrochemical process in the proton exchange membrane inside the hydrogen power supply; Based on the acoustic response signal sets, reconstruct the current density inside the proton exchange membrane to obtain a spatial distribution map; Analyze and generate an acoustic compensation signal according to the non-uniformity represented by the spatial distribution map, and apply it to the target area of the proton exchange membrane; Integrate the spatial distribution map, the acoustic compensation signal, and the externally input power demand signal and the feedback current output power to generate the final control instruction for regulating the hydrogen power supply.

[0007] The beneficial effects of the present application are that the real-time perspective and targeted physical intervention of the internal state of the fuel cell are achieved through non-invasive acoustic means, which significantly improves the dynamic response speed, power regulation accuracy, and operating stability of the hydrogen power supply. The related technical effects will be described in detail below. Attached Figure Description

[0008] Figure 1 This is a flowchart of the present invention.

[0009] Figure 2 This is a flowchart of the acoustic response signal set acquired by the present invention.

[0010] Figure 3 This is a flowchart of the process of reconstructing the spatial distribution map according to the present invention.

[0011] Figure 4 This is a flowchart illustrating how the present invention obtains the accurate time delay matrix. Detailed Implementation

[0012] Example 1: In this example, to achieve accurate diagnosis of the electrochemical processes within the proton exchange membrane, a physically reliable acoustic diagnostic system must first be constructed. For instance, this construction process overcomes the inherent limitations of idealized models in the art. Specifically, those skilled in the art might initially attempt to construct the diagnostic system based on a conventional electrostriction effect model and an idealized Nafion membrane sound velocity. For example, approximately 3 × 10⁻⁶ m² / m ... - ¹ 0 The design incorporated an electrostriction coefficient (β) of m / V and a sound velocity (c) of approximately 1200 m / s. However, through in-depth research, the inventors of this invention discovered that the actual electrostriction effect of the Nafion film is much weaker (β is approximately 10 m / V). -12 The sound velocity in the hydrated state is much lower (approximately 300-800 m / s) and exhibits significant anisotropy. This deviation from the original model results in extremely weak signals and low signal-to-noise ratios in the low-frequency range (e.g., 3-10 kHz), and the excessively long wavelength (λ=c / f) makes the physical diffraction limit much greater than 5 mm, which cannot meet the requirements for high-resolution imaging.

[0013] To fundamentally solve the above problems, this embodiment adopts a better system construction scheme. First, this embodiment identifies the electrochemical acoustic emission (AE) signal, which is characterized by suddenness and broadband nature, emitted by microbubbles formed from products (such as hydrogen and oxygen) during nucleation and collapse in the electrochemical reaction process. This signal is the main source of detectable sound waves within the membrane, and its signal intensity is much higher than that of the electrostriction effect.

[0014] Accordingly, in the configuration of system parameters, in order to effectively capture such high-frequency and wide-band AE signals, the core working frequency band of the diagnostic system is set to 40-60 kHz. Specifically, the piezoelectric sensor used has a wide-band response of 10 Hz-100 kHz, and the sampling rate of the data acquisition system is set to 200 kHz to meet the Nyquist criterion for distortionless acquisition of signals up to 100 kHz. In addition, the system is also configured with a low-noise amplifier with a preamplifier gain of 40 dB and a hardware anti-aliasing filter with a cutoff frequency of 80 kHz.

[0015] By modifying the diagnostic mechanism from the electrostrictive effect to the electrochemical acoustic emission and raising the working frequency band from low frequency to high frequency, the diagnostic system constructed in this embodiment not only solves the problem of weak signals, but also improves the physical resolution to the millimeter level by shortening the wavelength of the acoustic wave (λ=10 mm when c=500 m / s and f=50 kHz), laying a physical foundation for subsequent realization of enhanced resolution of 5 mm.

[0016] Embodiment 2 describes a method for accurate power regulation and overshoot suppression of a micro-grid hydrogen power supply.

[0017] Step 1, collect the acoustic waves generated by the electrochemical process in the proton exchange membrane inside the hydrogen power supply to obtain a set of acoustic response signals.

[0018] In this embodiment, the set of acoustic response signals, denoted as S _acoustic , refers to a set of time series signal data representing the acoustic activity inside the membrane synchronously collected by multiple sensors arranged on the proton exchange membrane.

[0019] Specifically, the implementation of this step is as follows: at the four corner points of the proton exchange membrane, one PVDF piezoelectric film sensor (e.g., thickness 28 μm, bandwidth 10 Hz-100 kHz) is arranged respectively. An external synchronous trigger signal (e.g., provided by a GPS clock, with a jitter of less than 10 ns) is used to ensure that the collection times of the four sensors are strictly aligned. Each sensor continuously records the original acoustic pressure signal s_raw_i(t) generated by the electrochemical acoustic emission effect at a sampling rate of 200 kHz, where i=1, 2, 3, 4. The collected original signals are sequentially subjected to 50 Hz notch filtering to remove power frequency interference, and then an effective frequency band is extracted through a band-pass filter (e.g., 20-80 kHz), and finally four preprocessed acoustic signals s_i(t) are obtained. The four signals {s1(t), s2(t), s3(t), s4(t)} collectively constitute the set of acoustic response signals S_acoustic.

[0020] In this embodiment, the purpose of this step is to non-invasively and high spatiotemporal resolution acquire the original physical information that can reflect the spatial distribution state of the electrochemical processes such as proton transport and product generation inside the membrane. High sampling rate and accurate synchronization are the basis for subsequent high-precision time delay calculation and acoustic tomography.

[0021] In some optional embodiments, the number of sensors can be more than four, for example, can be additionally arranged along the edge of the membrane or in the central region of the membrane to increase the density of acoustic ray coverage and thus improve the quality of the reconstructed image. The material of the sensor can also be PZT or other piezoelectric materials with wideband response.

[0022] Step 2, based on the set of acoustic response signals, reconstruct the current density inside the proton exchange membrane to obtain a spatial distribution map.

[0023] In this embodiment, the spatial distribution map, denoted as J(x, y, t), is a two-dimensional or three-dimensional digital image that can change with time, and the value of each pixel point on the image represents the local current density of the corresponding physical position of the proton exchange membrane at time t.

[0024] Specifically, the implementation of this step is a multi-stage processing flow. Based on the set of acoustic response signals S acoustic, the time delay of acoustic wave propagation is accurately extracted by calculating the cross-correlation function or cross-power spectrum between each pair of sensor signals, thereby constructing an accurate time delay matrix T delay. The accurate sound speed model c_final(x, y) considering anisotropy and water content gradient is applied, such as in embodiment 1. The time delay matrix T delay and the sound speed field c_final(x, y) are used as input, and an acoustic tomography algorithm (for example, algebraic reconstruction technique ART or filtered back projection FBP) is used for calculation to solve the inverse problem and obtain the initial current density spatial distribution map. Optionally, to further improve the detail recognition ability, a super-resolution reconstruction algorithm can also be applied to the initial distribution map to generate the final spatial distribution map J(x, y, t).

[0025] Since the root cause of power overshoot and other problems is the non-uniformity of current distribution in the membrane, this step can convert the non-directly observable electrochemical non-uniformity problem into a quantifiable and intuitive heat map, providing a decision basis for subsequent diagnostic analysis and targeted compensation.

[0026] In some optional embodiments, the acoustic tomography algorithm can be based on the synchronous iterative reconstruction technique (SIRT) in addition to ART or FBP. The super-resolution reconstruction algorithm can be based on compressed sensing or deep learning (such as SRCNN), or a combination of multiple algorithms.

[0027] Embodiment 3. Mainly describes the signal collection and multi-resolution spectral feature extraction process based on electrochemical acoustic emission. That is, the step of collecting the acoustic response signal set in Embodiment 2 is further elaborated from the aspects of the physical mechanism of the signal, the time-frequency domain processing method, and the fine extraction of the spectral features.

[0028] In this embodiment, one preferred implementation of Step 1 of Embodiment 2, i.e., the step of collecting the acoustic response signal set, includes: Step 1.1. Identify and capture the electrochemical acoustic emission signals emitted by the microbubbles formed by the products in the electrochemical process when they are nucleated and collapse, which have burst and broadband characteristics; and the captured electrochemical acoustic emission signals are taken as the main data components constituting the acoustic response signal set.

[0029] Among them, the electrochemical acoustic emission signal refers to the transient elastic wave generated due to the rapid release of local energy (e.g., the formation and rupture of microbubbles) driven by the electrochemical reaction. It is different from the bulk effect (i.e., electrostriction) generated by the material due to the change of electric field, and its signal characteristics are characterized by extremely short duration (e.g., less than 100 µs), concentrated energy, and very wide frequency spectrum (e.g., covering 20-40 kHz).

[0030] Specifically, to achieve the identification and capture of such signals, first set a dynamic detection threshold for the broadband acoustic signal s_i(t) collected by the sensor, for example, the threshold can be set to 3 times the standard deviation of the signal background noise. When the signal amplitude exceeds the threshold in a short time, it is considered that an acoustic emission (AE) event is captured. By analyzing the waveform characteristics of the event, such as its rise time, duration, and amplitude exponential decay characteristics, it can be distinguished from the stationary background noise or other interference signals.

[0031] In this embodiment, all the AE event signals identified and meeting the above-mentioned burst and broadband characteristics are collected as the most critical data components in the acoustic response signal set S_acoustic. The motivation of this step is that the inventors have found that, compared with the electrostriction effect, the electrochemical acoustic emission event has stronger correlation with the local current density (i.e., the reaction rate) and higher signal intensity. Therefore, shifting the focus of diagnosis from continuous and weak bulk effect signals to discrete and high-energy AE events can significantly improve the signal-to-noise ratio and accuracy of diagnosis.

[0032] Step 1.2, further, the captured electrochemical acoustic emission signal is transformed into time-frequency domain by using wavelet packet decomposition technology, and is separated into a group of time-frequency domain sub-band signals; and from the time-frequency domain sub-band signals, acoustic emission characteristic parameters representing each acoustic emission event are calculated and extracted, the parameters including the energy, duration and occurrence rate of the event; the acoustic emission characteristic parameters jointly constitute the structured input data for the subsequent reconstruction step.

[0033] Among them, the wavelet packet decomposition (WPD) can decompose the signal in full band, thereby providing high resolution in time domain and frequency domain at the same time, and is particularly suitable for analyzing non-stationary transient signals. The acoustic emission characteristic parameters refer to a group of numerical values extracted from the waveform of each AE event, which can quantify the physical characteristics of the event.

[0034] Specifically, the implementation of this step is: for each captured AE event signal, wavelet packet decomposition is applied for processing. For example, Daubechies-10 (db10) is selected as the wavelet basis function, and 5-layer decomposition is performed, so that the original signal is decomposed into 32 mutually orthogonal frequency bands. On this basis, the characteristic parameters of each AE event are calculated, and the exemplary calculation method is as follows: event energy E_AE=∫|s_event(t)| 2 dt, wherein s_event(t) is the signal waveform of the event, and the integral interval is the duration of the event; the event duration T_d can be defined as the time length from the first time the signal exceeds the threshold to the last time the signal is below the threshold; the event occurrence rate R_AE is defined as the number of AE events captured per unit time (for example, 10 ms). These calculated parameters {E_AE, T_d, R_AE} are structured and stored as input for the subsequent current density reconstruction algorithm.

[0035] Since the AE signal has typical non-stationary characteristics, the traditional Fourier transform will lose its time information. WPD can finely describe each AE event in time and frequency domains, so that the energy, duration and other key parameters of the event can be more accurately calculated, providing high-quality input data for establishing the quantitative relationship between acoustic characteristics and electrochemical reaction rate. As an optional solution, other time-frequency analysis methods such as short-time Fourier transform (STFT) or Hilbert-Huang transform (HHT) can also be used, but considering the calculation efficiency and adaptability to sudden signals, wavelet packet decomposition is the preferred solution of this embodiment.

[0036] Step 1.3, further, to extract accurate spectral features from the electrochemical acoustic emission signal, further comprises a multi-resolution analysis process: performing a wideband spectral analysis on the electrochemical acoustic emission signal to identify one or more candidate resonance frequency bands containing potential resonance peaks; only within the candidate resonance frequency bands, applying a spectral refinement technique to generate a fine power spectrum with ultra-high frequency resolution; and from the fine power spectrum, accurately locating the frequency and intensity of the resonance peaks to form a set of accurate resonance peak parameters.

[0037] wherein the spectral refinement technique refers to a class of algorithms that can perform spectral analysis within a specific frequency range of interest with much higher frequency resolution than the conventional Fourier transform, such as Chirp-Z transform (CZT).

[0038] Specifically, the implementation of this process is as follows: first, perform a standard FFT transform on the collected signal to obtain a wideband power spectrum covering the entire acquisition bandwidth (e.g. 0-100 kHz), which may have a relatively low frequency resolution (e.g. about 48.8 Hz). Through a peak search algorithm, identify the energy concentrated frequency band on the wideband spectrum, for example, find that 40-60 kHz is a candidate frequency band containing potential resonance peaks.

[0039] Instead of analyzing the entire bandwidth, only for the candidate frequency band of 40-60 kHz, the Chirp-Z transform is applied. Exemplarily, set the start frequency f_start=40 kHz, the stop frequency f_stop=60 kHz, and the calculation points M=2000, so that a fine power spectrum within this frequency band can be generated, with a frequency resolution as high as (60000-40000) / 2000=10 Hz.

[0040] On this fine power spectrum, the center frequency f_res (with an accuracy of ±1 Hz) and the peak intensity A_peak of the resonance peak can be accurately located by, for example, parabolic interpolation or Gaussian fitting. The motivation of this multi-resolution analysis process is to balance computational efficiency and analysis accuracy. The computational cost of high-resolution analysis of the full frequency band is high, and this two-step method uses the strategy of coarse positioning + refinement to concentrate computational resources on the most critical frequency band, so that accurate resonance peak parameters can be obtained under the premise of meeting real-time requirements, providing accurate frequency basis for subsequent acoustic compensation (such as resonance sound field application in embodiment 8).

[0041] Embodiment 4 describes the construction and iterative solving process of anisotropic multi-physical field coupled sound velocity model. This embodiment is used to elaborate in detail the acquisition process of the core premise required for the reconstruction of the spatial distribution map in embodiment 2, i.e. the acquisition process of the sound velocity benchmark.

[0042] Step 2.1, the step of reconstructing the spatial distribution map, is achieved by applying an anisotropic multi-physics coupled sound speed model; the model jointly models the anisotropic material properties of the proton exchange membrane, the water content gradient along its thickness direction, and the propagation direction and dispersion effect of the acoustic wave; thereby providing an accurate sound speed benchmark that dynamically changes with spatial position, direction, and frequency for the reconstruction process.

[0043] wherein the anisotropic multi-physics coupled sound speed model is a mathematical function, the sound speed value c calculated by which is not a constant, but a function of multiple physical variables, i.e. c = f(x, y, z, φ, T, f, θ,...), wherein (x, y, z) is the spatial position, φ is the propagation direction, T is the temperature, f is the frequency, θ is the water content, etc. Specifically, the preferred sound speed model c_coupled constructed in the embodiment can be expressed as the product or complex function combination of multiple physical effects.

[0044] An exemplary model is as follows: c_coupled(x, y, z, φ, f) = c_wet(θ(z)) * f_aniso(φ) * f_temp(T) * f_dispersion(f).

[0045] wherein: c_wet(θ(z)) is the relationship between sound speed and water content θ and membrane thickness position z, for example c_wet(θ) = 350*(1-θ)^2.3 + 650*θ^1.7, wherein 350 m / s and 650 m / s are the approximate sound speeds under dry and wet membranes respectively; at the same time, the water content θ itself is also a function of the thickness z, for example θ(z) = θ_surface*exp(-z / L_d), wherein L_d is the water diffusion characteristic length (such as 30 µm), which describes the water content gradient along the thickness direction. f_aniso(φ) is the anisotropy factor, φ is the angle between the sound wave propagation direction and the membrane plane, for example f_aniso(φ) = 1 + K_aniso*cos 2 (φ), wherein K_aniso is the anisotropy coefficient (such as 0.3), indicating that the sound speed of in-plane propagation is faster than that of thickness direction propagation. f_temp(T) is the temperature correction factor, for example f_temp(T) = 1 - K_T*(T-T_ref), wherein K_T is the temperature coefficient (such as 0.002 / °C). f_dispersion(f) is the dispersion effect correction factor, for example f_dispersion(f) = 1 + K_f*log(f / f_ref), wherein K_f is the dispersion coefficient (such as 0.06), indicating that the high-frequency sound wave propagation speed is slightly faster than the low-frequency.

[0046] Since the proton exchange membrane is a highly heterogeneous, anisotropic composite material, its acoustic properties are influenced by multiple factors. Any simplified, idealized sound velocity model (such as assuming a constant sound velocity) will bring huge model errors to acoustic tomography, leading to serious distortion of the reconstructed current distribution. The multi-physics coupling model constructed in this embodiment can provide a highly realistic sound velocity benchmark by jointly modeling these key physical effects, which is the guarantee for accurate reconstruction.

[0047] Step 2.2, further, the process of providing an accurate sound velocity benchmark, specifically a solving process achieved through iterative optimization, which includes: In an initial or updated current sound velocity field, a set of theoretically calculated time delays is calculated by ray tracing; the theoretically calculated time delays are compared with an actually measured accurate time delay matrix to obtain a time delay error quantifying the deviation between the model and reality; And the time delay error is input to apply an iterative optimization algorithm to calculate the update amount of the current sound velocity field; the above steps are repeatedly executed until the time delay error converges within a preset threshold, and the converged current sound velocity field is the final accurate sound velocity benchmark.

[0048] Among them, ray tracing is a method for calculating the propagation path and time of waves in a non-uniform medium, and its basic principle is Fermat's principle or Snell's law.

[0049] Specifically, the implementation of the iterative optimization process is as follows: the first step (initialization): divide the membrane area into a grid (for example, 12x12), and assign an initial sound velocity value to each grid based on the modified model of embodiment 1 to form an initial sound velocity field c 0 (x,y). The second step (forward calculation): in the current sound velocity field c^k(x,y) (k is the number of iterations), for each pair of sensors (i,j), apply the layered ray tracing algorithm to calculate the shortest propagation path of the sound ray from i to j, and integrate along the path to obtain the theoretically calculated time delay τ_calc_ij. The third step (error calculation): compare the theoretically calculated time delay τ_calc of all sensor pairs with an actually measured accurate time delay matrix T_delay (its acquisition method is as in embodiment 6), and calculate a total time delay error, for example, the error function can be defined as E=Σ(τ_calc_ij-τ_ij_refined) 2 . The fourth step (inverse update): input the time delay error E to apply an iterative optimization algorithm to calculate the update amount Δc k (x,y) of the current sound velocity field c k .

[0050] Preferably, the Levenberg-Marquardt (LM) algorithm can be used, and its update formula is ck+1 = c k - η * (J T * J + λ * I) -1 * J T e; where c is the sound speed field vector; η is the learning rate; J is the Jacobian matrix, representing the rate of change of the theoretical time delay with respect to the sound speed field; λ is the regularization parameter; I is the identity matrix; and e is the time delay error vector. Step 5 (Loop and Convergence): Repeat steps 2 to 4 until the total time delay error E converges within a pre-set threshold (e.g., E < 10 -6 ) or the change in the sound speed field between two consecutive iterations is less than a certain minimum value. The converged sound speed field c^k(x, y) is considered as the final, accurate sound speed reference c_final(x, y) that matches the actual measured data.

[0051] Since the sound speed field and the propagation path of the acoustic wave are interdependent unknowns, they cannot be solved directly. This iterative process, through a closed-loop feedback of assuming the sound speed field -> calculating the path and time delay -> comparing the error -> revising the sound speed field, continuously lets the model approach the observed data, and eventually finds a self-consistent sound speed distribution that best explains the measured time delay data. As an optional solution, the algorithm used for optimization can also be the conjugate gradient method or genetic algorithm, etc. Especially when the error function presents multi-extremum characteristics, global optimization algorithms such as genetic algorithms may perform better.

[0052] Embodiment 5 describes a super-resolution reconstruction process based on heterogeneous algorithm fusion. This embodiment provides a detailed and preferred description of how to break through the physical diffraction limit and achieve higher spatial resolution in the step of reconstructing the spatial distribution map in Embodiment 2.

[0053] Step 2.3, the step of reconstructing the spatial distribution map, further comprises: processing the set of acoustic response signals by an acoustic tomography method to first obtain an initial distribution map limited by the physical diffraction limit; and applying a super-resolution reconstruction algorithm to the initial distribution map to further improve the detail recognition ability based on the initial distribution map, thereby generating the final spatial distribution map.

[0054] where the physical diffraction limit refers to the minimum detail size that can be resolved by any imaging system based on wave theory, which is about half of the wavelength (λ / 2) used. Super-resolution reconstruction algorithms refer to a class of computational methods that can recover high-frequency information beyond the physical diffraction limit from low-resolution observations, thus generating higher-resolution images. Specifically, after applying the sound velocity model and acoustic tomography algorithm (e.g., FBP) in application example 4, an initial current density spatial distribution map J_initial(x,y) can be obtained. Limited by the preferred acoustic working frequency (e.g., 50 kHz) and the sound velocity in the film (e.g., 500 m / s) of the present application, the corresponding physical diffraction limit is about 10 mm / 2 = 5 mm. Considering the noise and model errors in the actual system, the effective spatial resolution of the initial distribution map can only be 8-10 mm.

[0055] This resolution is still insufficient for identifying and locating smaller scale non-uniform cold or hot spots. To solve this problem, the present embodiment further applies one or more super-resolution reconstruction algorithms for post-processing after obtaining the initial distribution map J_initial(x,y) as input. By introducing prior knowledge about the signal or image (e.g., sparsity, structure), super-resolution algorithms can mathematically recover the high-frequency spatial information lost in the physical acquisition process, effectively improving the detail recognition ability of the final spatial distribution map without changing the hardware system, reaching or approaching the target resolution of 5 mm.

[0056] Step 2.4, further, the process of applying super-resolution reconstruction algorithms specifically includes: using at least two heterogeneous super-resolution reconstruction algorithms based on different physical or mathematical principles to process the initial distribution map in parallel, respectively, to generate multiple candidate distribution maps; and weighting and fusing the multiple candidate distribution maps to collect the advantages of different algorithms to form the final spatial distribution map. Wherein, heterogeneous super-resolution reconstruction algorithms refer to algorithms based on different theoretical foundations, for example, compressed sensing type algorithms based on signal sparsity assumption and deep learning type algorithms based on data driven. Specifically, a preferred implementation is to use the following three heterogeneous algorithms in parallel to process the initial distribution map J_initial(x,y): Compressed sensing (CS) based algorithm: this algorithm assumes that the current density map is sparse in a certain transform domain (e.g., discrete cosine transform DCT domain or wavelet domain). A high-resolution image is reconstructed by solving an L1 norm minimization problem, which is mathematically expressed as: min||s||1subject to||y-ΦΨs||2<ε; where s is the sparse coefficient to be solved; y is the low-resolution observation value (i.e., J_initial); Ψ is the sparse basis matrix; Φ is the measurement matrix; ε is the noise tolerance.

[0057] MUSIC-based algorithm: This algorithm is suitable for identifying point-like or small-area clustered sound sources. By performing eigen-decomposition on the covariance matrix of sensor signals, the signal subspace and noise subspace are separated, and the spatial spectrum function P_MUSIC(θ) = 1 / [a H (θ)*E_n*En H *a(θ)] is constructed; where a(θ) is the direction vector; E n is the noise subspace. By searching for the peak value of this spectrum function, the position of the sound source (i.e., the high current density center) can be located with high accuracy.

[0058] Deep learning-based algorithm: For example, a super-resolution convolutional neural network (SRCNN) is used. This network learns the complex nonlinear mapping from low-resolution features to high-resolution images through offline training on a large number of low-resolution-high-resolution image pairs. In actual application, J_initial(x, y) is directly input into the trained SRCNN, and a high-resolution candidate image can be quickly generated.

[0059] After obtaining the candidate distribution maps J_CS, J_MUSIC, and J_DL generated by the above three algorithms respectively, they are weighted and fused to form the final super-resolution spatial distribution map J_SR. The fusion formula can be expressed as: J_SR = w1*J_CS + w2*J_MUSIC + w3*J_DL; where w1, w2, w3 are weight coefficients, and w1 + w2 + w3 = 1. Since different algorithms have their own advantages and application scenarios. For example, the CS algorithm performs well in image overall structure recovery, the MUSIC algorithm is extremely accurate in local hotspot positioning, and the SRCNN is good at restoring texture and details. By fusing their outputs, we can learn from the best and get a final image that is more robust and more accurate than any single algorithm result. As an optional solution, the weight coefficients w_k can be fixed or adaptive. For example, the weight can be dynamically adjusted according to a certain feature (such as sparsity or smoothness) of the initial distribution map, and a larger weight is assigned to the algorithm that is expected to perform better in this scenario.

[0060] Embodiment 6 describes a high-precision time delay matrix acquisition method. This embodiment describes in detail how to obtain a high-precision time delay matrix as a key input of the reconstruction process in Embodiment 5 before the reconstruction process begins.

[0061] Step 2.5, the step of reconstructing the spatial distribution map, the implementation basis of which is to first obtain an accurate time delay matrix by performing the following operations: calculating the cross-power spectrum between any two sensor signals in the acoustic response signal set, and performing linear fitting on the phase information of the spectrum to analyze a set of high-resolution initial time delay values; performing physical consistency checking on the set of initial time delay values, the checking including symmetry checking and triangle inequality checking to identify and mark abnormal time delay data that does not conform to physical laws; and correcting or removing the abnormal time delay data to form an accurate time delay matrix for subsequent acoustic tomographic imaging calculation.

[0062] wherein the cross-power spectrum G_ij(ω) is the Fourier transform of the cross-correlation function of the two signals s_i(t) and s_j(t), which contains the phase difference information of the two signals at each frequency. The physical consistency checking is a process of verifying the effectiveness of the measured data according to basic physical laws (such as reciprocity and path relationship of sound wave propagation). Specifically, the implementation of the process is as follows: Initial time delay value analysis: for any pair of sensor signals s_i(t) and s_j(t), first calculate their Fourier transforms S_i(ω) and S_j(ω), respectively. Then calculate their cross-power spectrum G_ij(ω)=S_i(ω)*S_j*(ω); wherein S_j*(ω) is the complex conjugate of S_j(ω). Extract the phase spectrum φ_ij(ω)=angle[G_ij(ω)] of G_ij(ω). After phase unwrapping to eliminate 2π jumps, perform weighted linear least squares fitting on the phase spectrum within the effective frequency band of the signal (for example, 40-60 kHz), and the fitting equation is φ_ij(ω)=2πωτ_ij+φ0; wherein ω is the angular frequency; τ_ij is the initial time delay value to be solved; φ0 is the initial phase. By solving the slope of the linear equation, a high-resolution τ_ij can be obtained. The weight of the fitting can be taken as the modulus |G_ij(ω)| of the cross-power spectrum to enhance the contribution of the high frequency band with high signal-to-noise ratio.

[0063] Physical consistency checking includes: Symmetry checking: in theory, the propagation time from i to j should be equal to the propagation time from j to i, i.e. τ_ij=τ_ji. Calculate the asymmetry A_ij=|τ_ij-τ_ji| / max(|τ_ij|,|τ_ji|), if A_ij exceeds a small threshold (for example, 0.05), it is considered that the measurement pair is abnormal, which can be corrected to the average value (τ_ij+τ_ji) / 2 or removed.

[0064] Triangle inequality check: for any three sensors i, j, k, the propagation time delay between them should satisfy the triangle inequality |τ_ij|≤|τ_ik|+|τ_kj|. By traversing all possible sensor triplets, if a violation of this inequality is found, one or more of the time delay values in the triplet is marked as abnormal data.

[0065] Abnormal data handling: for all time delay data marked as abnormal, one can either correct or discard them. Correction can be done by a local or global optimization algorithm that minimizes the adjustment to the original measurements while satisfying all physical consistency constraints. Discarding is to mark the abnormal data as invalid and not use them in the subsequent tomographic reconstruction.

[0066] Since the accuracy of tomographic reconstruction is highly dependent on the accuracy of the input time delay data, the traditional method based on cross-correlation function peak positioning is limited by the sampling rate. The method based on mutual power spectrum phase fitting used in this embodiment utilizes the phase information of the entire frequency band and can achieve sub-sampling accuracy time delay estimation. More importantly, the subsequent physical consistency check step can effectively identify and handle gross errors caused by noise, multipath effects, etc., which is equivalent to adding a quality inspection process for the data entering the reconstruction algorithm, improving the stability and reliability of the final reconstruction result.

[0067] Embodiment 7, describes an interface mass transfer enhancement compensation method based on acoustic micro-flow.

[0068] Step 3.1, the step of generating and applying a directional acoustic compensation field, one implementation includes: identifying low current density regions from the spatial distribution map; and for the low current density regions, generating and focusing a beam of sound working in the medium frequency band; the acoustic beam induces acoustic micro-flow at the interface between the catalyst and the proton exchange membrane in the region to enhance the local mass transfer process in the region.

[0069] Among them, acoustic micro-flow is a stable fluid vortex driven by acoustic energy dissipation in the acoustic field. In the context of the modification scheme of the present embodiment, the medium frequency band is preferably 45-55 kHz.

[0070] Specifically, this step first uses a region recognition algorithm (such as the DBSCAN algorithm of embodiment 11) to analyze the J_SR(x, y, t) spatial distribution map generated in embodiment 5, automatically identifying and locating all regions with current density significantly lower than the average value, i.e. low current density region R_low. Subsequently, the acoustic actuator array of the control system (i.e. the piezoelectric sensor used to collect signals in embodiment 2, which is now working in transmission mode) generates a high-frequency acoustic beam focused on the center of R_low region.

[0071] The study shows that a common cause of local low current density is water flooding at the interface between the catalyst and the proton exchange membrane, i.e. too much liquid water forms a too thick mass transfer boundary layer, which hinders the diffusion of reactants (such as oxygen) to the active sites of the catalyst. The focused acoustic beam applied in this embodiment, whose energy is absorbed and scattered at the interface, will induce strong acoustic-induced microflows. This microflow, like countless micro-stirrers, can effectively disturb and thin the stagnant mass transfer boundary layer, accelerate the supply of reactants and the discharge of products, thereby directly strengthening the local mass transfer process in this region and improving its current density.

[0072] Step 3.2, further, the step of generating and focusing an acoustic beam working in the medium frequency band, specifically includes: for focusing, pre-record the acoustic response of the low current density region to obtain the channel information of the medium, and determine a set of phase delay parameters for adaptive focusing based on the time reversal of the channel information; for generating, calculate and generate a pulse width modulation control sequence according to the compensation needs of the low current density region to dynamically adjust the output power of the acoustic beam; and combine the phase delay parameters with the pulse width modulation control sequence to drive the acoustic actuator array to generate the final medium frequency acoustic beam.

[0073] Wherein, time reversal is to record the wave emitted by a source at a certain place, and then emit the time-reversed recorded signal back from the recording place, and the wave front will automatically converge to the original source point. Pulse width modulation (PWM) is a technology that controls the average power of a signal by changing the width of its pulses.

[0074] Specifically, the implementation of focusing and generating is as follows: Focusing (phase determination): a preferred implementation is to first emit a broadband probe pulse from one actuator (e.g. actuator 1) in the actuator array. Other actuators (2, 3, 4) record the response signals h_i1(t) (i=2, 3, 4) scattered and propagated through the target region R_low, which represents the channel impulse response from 1 to i through R_low. Then, time-reverse these response signals, i.e. h_i1(-t). When all actuators emit these time-reversed signals as excitation signals at the same time, according to the time reversal principle, the acoustic waves will be adaptively focused in the R_low region. The relative phase extracted from these h_i1(-t) signals, i.e. constitutes the required phase delay parameters {φ1, φ2, φ3, φ4}.

[0075] Generation (power regulation): a proper PWM duty cycle D is calculated for the compensation demand (e.g. the difference of its current density value from the mean value) in the R_low region. For example, an adaptive regulation strategy can be adopted: D(x, y) = D_base + D_var * (1 - J(x, y) / J_mean); where D_base is the base duty cycle (e.g. 15%); D_var is the variable range (e.g. 35%); J(x, y) and J_mean are the local and mean current density, respectively. The duty cycle is used to modulate a carrier wave with a frequency of 100 Hz, generating the PWM control sequence. Finally, a high-frequency sinusoidal carrier wave (e.g. 50 kHz) with a phase delay φ_i is modulated with the PWM control sequence to obtain the final signal driving the actuators. Time reversal is a very robust adaptive focusing method without the need of knowing the specific sound speed distribution of the medium, and is particularly suitable for PEM, a highly inhomogeneous and time-varying complex medium. With PWM for power regulation, the energy is allocated on demand, ensuring the compensation effect while avoiding unnecessary energy consumption and heat generation, improving the overall energy efficiency of the system.

[0076] Step 3.3, further, the step of generating and focusing the mid-frequency acoustic beam, further comprises a compensation effect prediction step: based on the parameters of the mid-frequency acoustic beam, the induced interfacial flow velocity distribution at the interface of the catalyst and the proton exchange membrane is calculated; and the interfacial flow velocity distribution is taken as input, and a mass transfer enhancement prediction model is applied to quantitatively calculate a predicted current density enhancement distribution; the predicted current density enhancement distribution is used to assist the generation of the final control instruction.

[0077] Specifically, the implementation of the prediction step is as follows: First, based on the determined acoustic beam parameters (frequency, amplitude and phase of each actuator), the acoustic pressure field p(x, y) and the acoustic pressure gradient grad p formed in the target region are calculated through an acoustic simulation model.

[0078] Then, a theoretical model of acoustic streaming (e.g. Rayleigh flow theory) is applied to calculate the induced interfacial flow velocity distribution v_interface(x, y), an example formula is: v_acoustic(x, y) = -1 / (2ρω)*Re[p**grad p]; where ρ is the medium density; ω is the angular frequency; p* is the acoustic pressure complex conjugate.

[0079] After the interface velocity distribution is obtained, a mass transfer enhancement prediction model is applied to calculate the change of current density. For example, the thinning effect of the mass transfer boundary layer thickness caused by acoustic streaming can be calculated first: δ_new(x, y) = δ0 / (1 + Sh_acoustic); where δ0is the initial boundary layer thickness; Sh_acoustic is the acoustic Sherwood number, which is positively correlated with v_interface.

[0080] Finally, the predicted current density enhancement distribution can be given by the following formula: ΔJ_predicted(x, y) = J0*(δ0 / δ_new-1); where J0is the reference current density under mass transfer limited conditions.

[0081] The method provides a feedforward information. Before the actual application of the acoustic field, the system can estimate the positive effect brought by the operation. ΔJ_predicted(x, y) can be used in a higher level control loop (such as the controller of embodiment 10), for example, the system can adjust the supply amount of hydrogen or air in advance according to the predicted current increase, so as to realize the collaborative control among multiple actuators, avoid the lag and instability caused by single adjustment, and finally achieve faster and more stable power response.

[0082] Embodiment 8 describes a proton conductivity regulation compensation method based on acoustic elastance effect, which is another optional implementation of embodiment 7 and is the second acoustic compensation path in the present application.

[0083] Another implementation of step 3.4, generating and applying a directional acoustic compensation field, includes: determining and obtaining a mechanical resonance frequency of the proton exchange membrane in advance by frequency scanning; identifying the low current density area from the spatial distribution map; and applying a resonance acoustic field to the low current density area with the mechanical resonance frequency as the working frequency; the resonance acoustic field in the area uses the acoustic elastance effect to dynamically modulate the proton transition barrier on the proton transport path, thereby directly regulating the proton conductivity of the area.

[0084] The mechanical resonance frequency is the frequency point at which the vibration amplitude of the membrane reaches the maximum when the external excitation frequency is equal to a certain inherent vibration frequency of the proton exchange membrane. The acoustic elastance effect refers to the alternating stress field generated by the propagation of acoustic waves in a medium, which can cause changes in physical properties such as elastic constant of the medium. In this embodiment, it specifically refers to the modulation effect of the stress field on the proton transition barrier.

[0085] Specifically, the implementation of this step is: through a frequency scanning program, the proton exchange membrane is excited in the range of 7-10 kHz (this is an example range, the actual frequency depends on the specific size and material properties of the membrane), and its acoustic response amplitude is monitored. Its mechanical resonance frequency f_membrane, such as 8.5 kHz, can be accurately determined by methods such as Lorentz fitting. As in Example 7, the low current density region R_low is identified. The acoustic actuator array is driven to generate and apply a resonant acoustic field with the resonance frequency f_membrane as the operating frequency.

[0086] The present method provides a different compensation mechanism from Example 7. Example 7 indirectly improves the current by improving mass transfer, while the present embodiment directly regulates the proton mobility μ_H+ to achieve compensation. At the resonance frequency, a small driving energy can excite a high-intensity mechanical stress field in the membrane. According to the acoustoelastic effect, this stress field can periodically change the energy barrier that protons need to overcome when jumping in the hydration channel, which is equivalent to providing assistance for proton transport, directly improving the proton conductivity and current density in this region. This method may have a better compensation effect for the problem of intrinsic conductivity decline caused by membrane aging, structural densification, etc.

[0087] Step 3.5, further, the step of applying a resonant acoustic field further comprises: constructing a target function aimed at maximizing the compensation effect in the low current density region; and using a gradient optimization algorithm to solve the target function to obtain an optimized spatial phase distribution; the finally generated resonant acoustic field carries the optimized spatial phase distribution, thereby achieving selective compensation of the low current density region.

[0088] Wherein, the spatial phase distribution φ(x,y) is a function defined on the membrane plane, which specifies a specific phase relationship for the sound waves emitted from each actuator at each point in space, and can sculpt a synthesized acoustic field with a specific spatial form.

[0089] Specifically, one preferred implementation is to construct and maximize the following objective function J_target: J_target = ∫∫_R_low μ_enhanced(x, y) dxdy - λ1 * ∫∫_R_high μ_enhanced(x, y) dxdy; where μ_enhanced(x, y) is the protonic conductivity enhancement factor predicted from the sound field intensity; R_low and R_high are the low-current region that needs to be enhanced and the high-current region that needs to be avoided, respectively; λ1 is a weight coefficient (e.g., 0.5) to balance the two goals of target region enhancement and non-target region suppression. To solve this objective function, a gradient optimization algorithm can be used, such as the conjugate gradient method. By calculating the gradient of J_target with respect to the spatial phase distribution φ(x, y), and iteratively updating φ(x, y) in the gradient direction until J_target converges to a maximum value. The final optimal spatial phase distribution φ_opt(x, y) is loaded onto the signals driving each actuator.

[0090] The method can achieve selective and targeted compensation. If a uniform resonant sound field is simply applied, although it can enhance the performance of the low-current region, it will also enhance the high-current region, which may cause local overheating or exacerbate the Matthew effect. By optimizing the spatial phase distribution, the sound wave energy (i.e., stress field intensity) can be redistributed in space, achieving constructive interference (energy concentration) in the R_low region and destructive interference (energy cancellation) in the R_high region, thereby precisely heating only the cold spots that need to be compensated for, avoiding damage to hot spots.

[0091] Step 3.6, further, the step of using acoustic elastance effect to dynamically modulate proton transition barrier is achieved through a refined physical modeling process, which includes: based on the parameters of the resonant sound field, first calculate the acoustic-induced strain field induced by it in the proton exchange membrane; and combined with a pre-set effective acoustic elastic coefficient, the specific change value of the proton transition barrier is quantitatively calculated through the acoustic-induced strain field.

[0092] wherein, the acoustic-induced strain field ε_acoustic(x, y) is the periodic mechanical deformation field caused by the alternating pressure of the sound wave in the membrane. The proton transition barrier E_a is the activation energy required for the proton to jump from one water molecule cluster to the next in the Grotthuss transport mechanism.

[0093] Specifically, the modeling process is implemented as follows: based on the known resonant acoustic field parameters (e.g. the spatial distribution of acoustic pressure amplitude P(x, y)), the acoustic-induced stress field σ_acoustic(x, y) is calculated using the acoustic constitutive equation. Through the compliance tensor of the membrane (which can be determined by Young's modulus and Poisson's ratio), the acoustic-induced strain field ε_acoustic(x, y) is calculated from the stress field.

[0094] Applying the acoustic-elastic theory, the change of proton transition barrier ΔE_barrier is calculated. An exemplary calculation formula is: ΔE_barrier(x, y) = γ_eff * σ_acoustic(x, y) * V_m; where γ_eff is the effective acoustic-elastic coefficient, which is a material constant representing the strength of stress modulation on the barrier (e.g. 0.15 eV / GPa for high-frequency acoustic waves); V_m is the molar volume of the proton exchange membrane.

[0095] The refined modeling process provides a quantitative mathematical description of the physical mechanism of acoustic compensation, directly linking the macroscopic acoustic field parameters to the microscopic proton transport behavior. Through this model, the control system not only knows that the acoustic field can act, but can accurately predict how much the proton transition barrier will be reduced and the proton mobility will be improved under a given acoustic pressure. This quantitative prediction capability is the basis for precise closed-loop control, making acoustic compensation from an empirical adjustment method to a modeled, accurately measurable control input.

[0096] Embodiment 9, describes the control parameter adaptive optimization and knowledge transfer process based on reinforcement learning.

[0097] Step 4.1, the step of generating the final control instruction, further includes an adaptive process of optimizing its internal control parameters online: periodically collecting and calculating a set of performance evaluation indicators representing the overall operation effect of the method; using the performance evaluation indicators as input, applying a reinforcement learning algorithm to generate a set of updated control parameters; and using the updated control parameters to adjust the generation logic of the final control instruction.

[0098] The core idea of reinforcement learning (RL) is to let an agent learn how to maximize cumulative reward through trial and error in the interaction with the environment, and thus find the optimal behavior policy. The control parameters refer to the tunable coefficients in the controller, such as the proportional gain K_p and integral gain K_i of a PI controller, and the strength and frequency of acoustic compensation, etc. Specifically, the adaptive process is a loop that runs continuously. After a complete power regulation event (e.g., a power step response), the system calculates a set of performance evaluation indicators E, for example, the indicator set can include the final non-uniformity η_final, the power overshoot P_overshoot, the regulation time t_settling, and the steady-state fluctuation σ_ss. This indicator set E is input into a RL agent as part of the state information of the environment feedback.

[0099] The agent (e.g., a deep neural network) outputs an action according to the current state and its internal policy, which is a set of updated control parameters {K_p_new, K_i_new,...}. These new parameters will replace the old parameters for the next power regulation task.

[0100] The motivation of this adaptive process is that the hydrogen power supply is a complex nonlinear system, whose characteristics will drift with the changes in working conditions such as aging, temperature, humidity, etc. Any fixed, pre-tuned control parameters cannot remain optimal under all conditions. The RL online optimization process introduced in this embodiment makes the controller have a living characteristic, which can continuously and autonomously fine-tune its control strategy according to the performance of each task (i.e., performance indicators), so that the system performance always remains near-optimal in the long run.

[0101] Step 4.2, further, the process of applying the reinforcement learning algorithm specifically includes: according to the performance evaluation indicators, a scalar reward value is calculated through a pre-set comprehensive reward function, which aims to balance multiple optimization objectives such as power overshoot, regulation time, and running stability; the current system state, the performed action, and the calculated scalar reward value are collectively encapsulated into a state-action-reward triple, and stored in an experience replay buffer; and historical data is sampled from the experience replay buffer to form a training batch, which is used to train and update the internal model of the reinforcement learning algorithm.

[0102] ExperienceReplayBuffer is a FIFO queue storing historical interaction records (state-action-reward-next state quadruples), the RL agent does not directly learn from the current interaction data, but randomly samples a batch of historical data from the buffer to train, in order to break the data correlation and improve learning efficiency. Specifically, a preferred implementation is to use the Deep Deterministic Policy Gradient (DDPG) algorithm. In each learning cycle, the system performs the following operations: Calculate reward: According to the performance indicator E obtained in the last step, a scalar reward value r_t is calculated through a comprehensive reward function r=f(E). The specific form of this function will be described in the next step.

[0103] Package and store experience: Package the key information of this interaction into an experience tuple (s_t, a_t, r_t, s_{t+1}). Among them, s_t is the initial state of this task (for example, the initial non-uniformity η_initial and the power step amplitude ΔP); a_t is the control parameter used in this task (for example, {K_p, K_i,...}); r_t is the calculated reward value; s_{t+1} is the state at the end of the task. Store this tuple in an experience replay buffer with a fixed size (for example, capacity 1000).

[0104] Sampling and training: When the amount of data in the buffer is sufficient, a small batch of data (for example, 32 historical experiences) is randomly sampled from it. Use this small batch of data to train the internal models of the DDPG algorithm, namely the Actor network (policy network) and the Critic network (value network). This process provides a stable and efficient learning mechanism for the RL algorithm. The comprehensive reward function converts the multi-objective optimization problem into a single objective (maximize cumulative reward) problem, providing a clear direction for learning. The experience replay mechanism disrupts the order of experiences, so that each training sample is approximately independent and identically distributed, which greatly stabilizes the training process of deep neural network-based RL algorithms and avoids learning shocks or non-convergence problems caused by high correlation of data.

[0105] Step 4.3, further, the structure of the comprehensive reward function when calculating the scalar reward value includes: a linearly weighted combination of multiple optimization objectives in the performance evaluation indicator to form a basic reward part; and a nonlinear penalty term whose mechanism is: when one or more key indicators in the performance evaluation indicator exceed the preset safety threshold, the value of this penalty term will grow exponentially with the degree of exceeding; the final scalar reward value is determined by the basic reward part and the nonlinear penalty term.

[0106] Specifically, one preferred comprehensive reward function r_total can be designed as follows: r_total = r_base + r_penalty. Here, the base reward part r_base is a linear weighted sum of multiple performance indicators, aiming to guide the agent to make trade-off optimization among multiple objectives within a safe range.

[0107] For example: r_base = -w1 * η_final - w2 * (P_overshoot / P_rated) - w3 * (t_settling / t_target) - w4 * σ_ss; where w1, w2, w3, w4 are the weight coefficients of each term (e.g., 0.3, 0.3, 0.2, 0.2, respectively); P_rated and t_target are the rated power and target regulation time for normalization. The negative signs here indicate that we want to minimize these indicators.

[0108] The nonlinear penalty term r_penalty is a safety guardian. For example, when the power overshoot P_overshoot exceeds a preset safety threshold P_threshold (e.g., 15% of the rated power), this penalty term is activated. Its mathematical form can be: r_penalty = -exp((P_overshoot - P_threshold) / K_exp); where K_exp is a coefficient that controls the exponential growth rate.

[0109] This reward function combines conventional optimization and risk aversion. In the normal working range, r_penalty is zero, and the agent's behavior is dominated by r_base, which will strive to find an optimal balance point between response speed, stability, and uniformity. But once its reckless exploration behavior causes the system to touch the safety red line (e.g., a huge power overshoot occurs), the exponential penalty term will give a huge negative reward, and this sharp pain will quickly and powerfully correct the agent's strategy, making it stay away from the danger zone in subsequent exploration. This design ensures that the RL-driven adaptive process not only pursues excellence but also strictly adheres to the safety bottom line.

[0110] Step 4.4, further, the adaptive process further includes a step of knowledge transfer and quick invocation: the working condition features representing the current operating conditions, including the initial non-uniformity and power step amplitude, are constructed as a multi-dimensional index; and the updated control parameters generated by the reinforcement learning algorithm are stored in the positions of the multi-dimensional parameter lookup table corresponding to the multi-dimensional index; so that when similar working condition features are encountered in the future, the optimal control parameters can be directly retrieved and invoked from the multi-dimensional parameter lookup table and applied to the generation of the final control command.

[0111] ​Wherein, a multi-dimensional parameter lookup table (LUT) takes the discretized working condition feature combination as the address (index), and the corresponding optimal control parameter set as the storage content.

[0112] Specifically, the implementation of this step is to discretize the continuous working condition features. For example, the initial non-uniformity η_initial can be divided into 5 levels (such as [0-0.2], [0.2-0.4],...), the power step amplitude ΔP can be divided into 4 levels (such as [0-250W], [250-500W],...), and the ambient temperature T can be divided into 3 levels. These discretized levels together form a three-dimensional index [index_η, index_ΔP, index_T]. When the RL agent finds a set of optimal control parameters {K_p*, K_i*,...} through learning under a specific working condition, this set of parameters is stored in the corresponding cell of the three-dimensional LUT.

[0113] In future operation, when the system faces a new task, it will first obtain the current working condition features and convert them into the corresponding multi-dimensional index, and then directly query the LUT. If the corresponding cell of the index has stored historical optimal parameters, the system will directly call these parameters to execute the task without waiting for RL to perform new learning. If the corresponding cell is empty, the system will start the exploration and learning process of RL.

[0114] Knowledge transfer and quick call solve the problem of slow response that may exist in pure RL learning, and structure the accumulation and sedimentation of the results of the learning process. It builds a long-term memory system (LUT), which makes the ability of the RL agent to draw inferences from specific instances concrete. When encountering new problems, the experience notes (LUT) are consulted first, and if there are records, they are directly used, and if there are no records, the exploration and learning are used to fill in the blank, and the new knowledge is stored in the notes. This mechanism combines the fast response of the lookup table method and the adaptive innovation ability of reinforcement learning, and realizes the continuous accumulation and efficient reuse of knowledge.

[0115] Embodiment 10, describes a time scale coordinated control architecture and a dynamic overshoot suppression method.

[0116] Step 4.5, the step of generating the final control instruction, is implemented through a multi-time scale coordinated control architecture, which includes: a fast feedforward loop running at a millisecond time scale, which predicts the short-term change trend of power by analyzing the gradient of the spatial distribution map, and generates an active feedforward compensation component accordingly; and a slow feedback loop running at a hundred-millisecond time scale, which generates a corrective feedback correction component through a proportional-integral controller according to the steady-state deviation between the power demand signal and the current output power; the final control instruction is integrated from the feedforward compensation component and the feedback correction component.

[0117] In the multi-time scale coordinated control architecture, complex control tasks are divided into control loops of different rates, and each loop focuses on handling disturbances and errors within a specific dynamic range. Specifically, the implementation of the architecture is as follows: Fast feedforward loop: This loop runs at a very high frequency (e.g., period less than 10 ms, preferably 1 ms), and its core function is prediction and pre-emptive compensation. It does not rely on actual power output error, but directly uses the current density spatial distribution map J_SR(x, y, t) generated in embodiment 5. By calculating the gradient grad J of this map, the area where the current is about to change dramatically can be identified. Based on this gradient, a prediction model can estimate the possible change trend ΔP_pred of power in the near future (e.g., the next 10 ms). An exemplary prediction model is: ΔP_pred(t+T)=∫∫κ*grad J_comp*v_drift*dxdy; where κ is an empirical proportionality coefficient (e.g., 0.8); v_drift is the equivalent drift speed of protons in the membrane. Based on this prediction value, a feedforward compensation component C_ff=K_ff*ΔP_pred is generated; where K_ff is the feedforward gain.

[0118] Slow feedback loop: This loop runs at a lower frequency (e.g., period of 50-100 ms), and its core function is to correct and eliminate steady-state errors. It uses a traditional proportional-integral (PI) control strategy. First, calculate the steady-state deviation e_ss between the power demand signal P_demand(t) and the actual output power P_output(t). Then, generate a feedback correction component C_fb=K_p*e_ss+K_i*∫e_ss*dt through a PI controller; where K_p and K_i are the proportional and integral gains, respectively.

[0119] Instruction integration: The final control instruction is integrated from the two components: C_final=C_ff+C_fb.

[0120] The multi-time scale architecture combines the rapidity of feedforward control and the accuracy of feedback control. The internal dynamics of the hydrogen power supply (on the order of milliseconds) are much faster than its external power response (on the order of hundreds of milliseconds). Slow feedback loops cannot cope with fast internal state changes, inevitably leading to lag and overshoot. The fast feedforward loop, by directly seeing the internal state (current gradient), can actively compensate for power errors before they actually occur, effectively suppressing the dramatic fluctuations in the transient process. The slow feedback loop is responsible for handling the final steady-state errors caused by model prediction inaccuracies or external slowly varying disturbances, ensuring the accuracy of the system over a long period of time.

[0121] Step 4.6, further, the control architecture further comprises an overshoot dynamic suppression process: real-time monitoring of the current output power rise rate, when it exceeds the preset threshold, triggering an overshoot warning; after triggering the overshoot warning, performing overshoot cause diagnosis to distinguish whether the overshoot is caused by the non-uniformity of the spatial distribution map or by the integral saturation of the control algorithm, and outputting an overshoot cause type; and according to the overshoot cause type, executing the corresponding targeted suppression strategy, the strategy including: when judged as caused by spatial distribution non-uniformity, enhancing the compensation strength of the directional acoustic compensation field; when judged as caused by integral saturation, performing anti-saturation processing on the integral term of the controller. Wherein, the overshoot cause diagnosis is a multi-source information-based logical reasoning process, aiming to determine the root cause of the rapid rise in power and the possible triggering of overshoot. Specifically, the implementation of this suppression process is as follows: Triggering the warning: real-time calculation of the rate of change of the current output power dP_output / dt. When detecting that dP_output / dt exceeds a preset threshold (for example, 500 W / s), the system immediately triggers an overshoot warning flag.

[0122] Cause diagnosis: after the warning is triggered, the system immediately queries two key state quantities: one is the non-uniformity index η from the internal physical state (calculated by the method of embodiment 11); the other is the power deviation integral value ∫ΔP*dt from the controller state. The diagnostic logic can be a simple decision tree: if the value of η is very high at this moment (for example, η>0.6), it is determined that the overshoot is caused by spatial distribution non-uniformity; if the absolute value of ∫ΔP*dt is very large at this moment (for example, |∫ΔP*dt|>100 W*s), it is determined that the controller integral saturation causes the overshoot.

[0123] Targeted suppression: If it is determined that it is caused by spatial distribution non-uniformity, the system will immediately execute a physical layer suppression strategy, such as dynamically increasing the intensity (or power) of the currently applied acoustic compensation signal (such as embodiment 7 or seven) by 30%.

[0124] If the controller integral saturation is determined to be the cause, the system will immediately execute the suppression strategy at the algorithm level, i.e. Anti-Windup processing for the PI controller in the slow feedback loop, for example, the value of the integral term ∫e_ss*dt is forced to be limited within a pre-set reasonable range (e.g. [-50, 50] W*s), to prevent its excessive accumulation.

[0125] This process changes the overshoot suppression from a blind braking behavior to a precise intervention according to the symptoms. Different causes of overshoot require different solutions. Physical problems (non-uniformity) need physical means (enhanced acoustic compensation) to solve, while algorithmic problems (integral saturation) need algorithmic means (anti-saturation) to deal with. This diagnosis-based targeted suppression strategy avoids ineffective or incorrect adjustment actions, making the overshoot suppression process more efficient and stable.

[0126] Step 4.7, further, the step of performing overshoot cause diagnosis is implemented based on first constructing and analyzing a multi-dimensional overshoot feature vector, the construction process of which includes: calculating the power tracking error representing the current deviation, the power rate of change representing the dynamic trend, and the power deviation integral representing the historical cumulative effect from the current output power and power demand signal; and extracting the spatial non-uniformity index representing the physical source from the spatial distribution map; combining the above calculated and extracted multiple features to form a multi-dimensional overshoot feature vector.

[0127] Specifically, the multi-dimensional overshoot feature vector F_OS is constructed as follows: Power tracking error e_p(t): e_p(t)=P_output(t)-P_demand(t); it represents the size and direction of the deviation between the current output and the expectation.

[0128] Power rate of change dP / dt: calculated by differencing the time series of P_output(t), preferably using a five-point central difference formula to improve accuracy. It represents the instantaneous speed of power change.

[0129] Power acceleration d 2 P / dt 2 : calculated by differencing the dP / dt sequence again. It represents the trend of power change, whether it is accelerating or decelerating.

[0130] Power deviation integral ∫e_p*dt: cumulative integration of tracking error e_p(t). It represents the historical cumulative effect of the error, which is a direct measure of integral saturation.

[0131] Spatial non-uniformity index η: calculated from the latest current density spatial profile J_SR(x,y,t), η = std(J_SR) / mean(J_SR). It characterizes the severity of the physical root cause of the overshoot.

[0132] The final feature vector is F_OS = [e_p, dP / dt, d 2 P / dt 2 ,∫e_p*dt,η]. This vector is computed and updated in real-time. This step provides rich and comprehensive decision basis for the cause diagnosis in the previous step. A single index (e.g. looking at dP / dt only) is not enough to distinguish complex fault modes. This multi-dimensional feature vector, which combines information from both electrical control domain (power, error, integral) and internal physical domain (non-uniformity), constitutes a complete snapshot of the current system overshoot risk. Based on this information-rich vector, the diagnostic logic (whether it is a simple decision tree, or a more complex machine learning classifier such as fuzzy logic or support vector machine) can make much more accurate and robust judgments than those based on a single index.

[0133] Embodiment 11, describes a method for accurate identification and quantification of non-uniformity regions based on density clustering and multi-dimensional features. This embodiment describes in detail the processing of identifying low current density regions and quantifying non-uniformity in the present invention.

[0134] Step 2.6, a preferred non-uniformity evaluation and region identification method, includes: Extracting multi-dimensional statistical features, i.e. after obtaining the super-resolution current density profile J_SR(x,y,t) in Embodiment 5, not only the mean and standard deviation are calculated, but further higher-order statistical moments are calculated to describe the distribution pattern. Specifically, the skewness is calculated, which is used to characterize the degree of asymmetry of the distribution; and the kurtosis is calculated, which is used to characterize the sharpness or flatness of the distribution. These higher-order statistical quantities can more finely reflect the pattern of non-uniformity.

[0135] Identifying abnormal regions using density-based clustering algorithm: To overcome the shortcomings of traditional threshold segmentation method, which is sensitive to noise and prone to produce broken small regions, the embodiment preferably uses a density-based spatial clustering algorithm (DBSCAN) to automatically identify high and low current density regions. The DBSCAN algorithm is an unsupervised learning algorithm that identifies clusters based on the local density of sample points, without the need to specify the number of clusters in advance. Specifically, each pixel point of J_SR(x, y, t) is regarded as a point p_i = [x_i, y_i, J_i] in a three-dimensional space, where J_i is the current density value of the point. Two core parameters of the DBSCAN algorithm are set: neighborhood radius ε (for example, set to 7.5 mm corresponding to the physical distance) and minimum point number MinPts (for example, 4). When the algorithm runs, all pixel points are automatically divided into core points, boundary points and noise points, and naturally form several clusters. Then, the average current density of each cluster is calculated and compared with the global average current density. The cluster with an average value significantly higher than the global average is identified as a high current region R_high, and the cluster with an average value significantly lower than the global average is identified as a low current region R_low. Since the physical cold or hot spots usually exhibit spatially continuous and abnormally dense regions, the core idea of the DBSCAN algorithm is highly consistent with this physical reality, so the boundary of the region identified by the DBSCAN algorithm is smoother and more complete, and it is not sensitive to noise, and the result is much better than simple threshold segmentation.

[0136] Constructing a comprehensive non-uniformity feature vector, including: integrating the information of multiple dimensions obtained in the above steps to construct a comprehensive non-uniformity feature vector U for high-level control decision. An exemplary U vector can include: U = [η_weighted, Kurtosis, Skew, A_high, A_low, grad_max, θ_principal]; where η_weighted is the weighted non-uniformity considering the weight of the central region; A_high and A_low are the area ratios of the high and low current regions; grad_max is the maximum value of the current gradient; θ_principal is the main direction of non-uniformity obtained by inertia tensor analysis. This comprehensive feature vector U provides a more comprehensive and detailed description of the internal state of the system for the reinforcement learning agent in embodiment 9 or the overshoot diagnosis logic in embodiment 10.

[0137] Embodiment 12, a method for describing the coherent synthesis and adaptive power allocation of the mixed mode compensation signal.

[0138] Step 3.7, one preferred compensation signal generation method, includes hybrid mode signal synthesis and adaptive power allocation. The hybrid mode refers to a compensation state that activates both interface mass transfer enhancement (S3a mode) and proton mobility modulation (S3b mode) simultaneously. Specifically, the implementation of this method is as follows: Hybrid mode signal synthesis: when the control system decides to enable the hybrid mode according to the diagnostic results (e.g., the value of non-uniformity η is in an intermediate zone, or its eigenvector U exhibits a composite feature of multiple causes), the compensation signals of the two modes need to be weighted and synthesized. The final driving signal U_combined_i(t) of the ith actuator can be expressed as: U_combined_i(t) = w_a*U_flow_i(t) + w_b*U_field_i(t); where U_flow_i(t) is the signal of S3a mode (as in Example 7); U_field_i(t) is the signal of S3b mode (as in Example 8); w_a and w_b are dynamic weighting coefficients, and w_a + w_b = 1. These two weight coefficients can be adaptively adjusted according to the non-uniformity η, and an exemplary adjustment strategy is: w_a = exp(-k*η) and w_b = 1-w_a; where k is the adjustment factor (e.g., k = 2). This strategy means that when the non-uniformity is mild, the S3a mode with lower energy consumption is mainly used; as the non-uniformity intensifies, the weight of the S3b mode gradually increases. To avoid intermodulation distortion that may occur when different frequency signals are synthesized, the system will monitor the frequency spectrum of the synthesized signal in real time, and if significant energy appears at the intermodulation frequency points such as 2f1-f2, the modulation depth will be appropriately reduced.

[0139] Adaptive power allocation: after signal synthesis, in order to achieve the best compensation effect under the limited total power budget (e.g., P_max = 1W), this embodiment adopts an adaptive power allocation strategy based on efficiency. First, the system needs to evaluate or query from historical data the compensation efficiency ε_i of each actuator i for the current compensation task, which is defined as the amount of target area current density improvement ΔJ_i / P_i that the actuator can bring with the consumption of unit power. Subsequently, the system solves an optimization problem: maximize the total current density improvement Σ(ε_i*P_i) under the condition of satisfying the total power constraint ΣP_i ≤ P_max. This optimization problem can be solved by the Lagrange multiplier method to obtain the optimal power allocation P_i_opt for each actuator channel. Finally, the signal amplitude of each channel is scaled according to this optimal allocation.

[0140] This method combines two compensation means with different physical mechanisms, which can more effectively deal with non-uniformity problems with complex causes. More importantly, the adaptive power allocation strategy based on efficiency is a precise energy management scheme. It can intelligently allocate the limited acoustic power budget to the currently most efficient actuator channel, thereby obtaining a compensation effect far superior to simple proportional scaling under the same energy consumption.

[0141] Embodiment 13 describes a method for tracking and dynamically updating the time evolution of the spatial distribution of current density. This embodiment aims to add a time dimension to the reconstruction of the spatial distribution map in Embodiment 5, providing a preferred implementation that is more suitable for transient process monitoring.

[0142] Step 2.7, a preferred method for reconstructing the spatial distribution map, further includes tracking and dynamically updating the time evolution of the map. Specifically, the implementation of this method is as follows: The system continuously generates the current density distribution map J_complete(x, y, t) with a fixed update period (e.g., Δt = 10 ms). To filter out random noise in single-frame images and observe long-term trends, the system uses an exponential weighted moving average (EWMA) filter to smooth the time series. The smoothed distribution map J_smooth(x, y, t) is calculated by the following formula: J_smooth(x, y, t) = α * J_complete(x, y, t) + (1 - α) * J_smooth(x, y, t - Δt); where α is the smoothing coefficient (e.g., 0.3), which determines the weight of new and old data.

[0143] After obtaining the smoothed time series, the system analyzes the dynamic change trend of the current density by calculating the time derivative ΞJ / Ξt. A simple calculation method is: ΞJ / Ξt = (J_smooth(t) - J_smooth(t - Δt)) / Δt. This derivative map ΞJ / Ξt can clearly show which areas in the membrane have rapidly rising or falling current density, and Ξ represents the partial derivative.

[0144] Based on the results of dynamic trend analysis, the system implements an adaptive update strategy. Specifically, the system sets a change rate threshold (e.g., 50 A / m 2 / s). For areas where |ΞJ / Ξt| exceeds this threshold, the system marks them as dynamic regions R_dynamic. Subsequently, the control system temporarily increases the update frequency for diagnosing and reconstructing these dynamic regions from the regular 10 ms to a higher 1 ms. For other static regions, the regular update frequency is maintained.

[0145] Since the power overshoot during the startup of hydrogen power is a typical transient problem, its internal state changes dramatically on the millisecond scale. The traditional monitoring method with fixed frequency updates may miss the key transient changes, leading to control lag. The time evolution tracking and dynamic update strategy proposed in this embodiment enables the system to have an attention mechanism. It can continuously monitor itself and intelligently focus valuable computing resources on the most dramatic changes in the key areas, capturing their dynamics with a higher frame rate. This strategy ensures that the system has stronger timeliness and pertinence when dealing with fast transient processes, and is another key guarantee for effectively suppressing overshoot.

[0146] Embodiment 14: A complete diagnostic-decision-compensation-control numerical calculation case This embodiment aims to quantitatively demonstrate the overall process of the present application through a specific, reproducible numerical case, thereby supporting all claims.

[0147] Suppose in a power step response scenario, the power demand P_demand is instantaneously increased from 200W to 1000W. Step A: Diagnosis process The sensors (1-4) located at the four corners of the 100mmx100mm membrane collect signals, and an accurate time delay matrix T_delay (unit: µs) is calculated by the method of embodiment 6, an example is as follows: T_delay=[[0,141.4,200.0,100.0],[141.4,0,100.0,200.0],[200.0,100.0,0,141.4],[100.0,200.0,141.4,0]].

[0148] Apply the sound speed model of embodiment 4 (assuming the average sound speed finally converges to c=707m / s), and apply the tomographic imaging and super-resolution algorithm of embodiment 5 to reconstruct the current density distribution J_SR(x,y,t). Suppose the reconstruction result shows that the current density in most areas is about 800A / m 2 , but there is an area near the coordinates (20mm, 80mm) with a current density of only 400A / m 2 .

[0149] Apply the method of embodiment 11 to calculate the global non-uniformity index η=0.45. At the same time, through the DBSCAN algorithm, the area centered at coordinates (20mm, 80mm) is accurately identified as a low current area R_low.

[0150] Step B: Decision and compensation calculation Since 0.3<η=0.45<0.5, the control system decides to enable the S3a mode (interface mass transfer enhancement) of embodiment 7.

[0151] The target focus is at (x_c, y_c) = (20, 80). Assume the four actuators are at (0, 0), (100, 0), (100, 100), and (0, 100) respectively. The working frequency f = 50 kHz. Calculate the distance d_i from each actuator to the focus, then the propagation time t_i = d_i / c, and finally the phase delay φ_i = 2πf*(t_i - t_min).

[0152] According to the current value of the R_low region 400 A / m 2 and the global mean 800 A / m 2 Apply the PWM duty cycle formula of Example 7 to calculate the required duty cycle D for this region.

[0153] Step C: Control command generation Apply the method of Example 10 to calculate the predicted power change ΔP_pred based on the gradient information in the reconstructed J_SR map, and then obtain the feedforward compensation component C_ff. Assume the calculation result is C_ff = 50 W.

[0154] The current power deviation e_ss = 1000 W - 200 W = 800 W. Assume the PI controller parameters K_p = 0.2 and K_i = 0.5, then in the first control period, the feedback correction component C_fb is mainly contributed by the proportional term, about 0.2 * 800 W = 160 W.

[0155] The initial control command C_initial = C_ff + C_fb = 50 W + 160 W = 210 W. This command is sent to the power regulation unit of the hydrogen power supply.

[0156] Step D: The system starts to respond by applying the acoustic compensation field (for local current enhancement of R_low) and the above electrical control command in parallel. Since the acoustic compensation directly acts on the source of non-uniformity, the η value quickly decreases. At the same time, due to the presence of the feedforward component, the control system can more quickly track the power demand. Finally, it is expected that the power overshoot P_overshoot can be controlled within 2% of 1000 W (i.e. < 20 W), the regulation time t_settling is shortened from the original 8 seconds to within 3 seconds, and the non-uniformity η_final is reduced to below 0.15. This numerical case clearly demonstrates how the various technical aspects of the invention work together to achieve precise power regulation and overshoot suppression.

[0157] Example 15, this embodiment aims to provide some specific, optional or optimized implementation details for the specific steps in the foregoing embodiments, which can be arbitrarily and selectively combined into the foregoing embodiments by those skilled in the art.

[0158] In the steps of power signal acquisition and processing in Embodiments 2, 9, etc., a preferred pre-processing method is to apply a Kalman filter to the raw power demand signal P_demand(t) and the current output power P_output(t) of the hydrogen power source respectively before calculating the real-time power deviation ΔP(t). Specifically, the variation process of the power signal is modeled as a linear dynamic system, and the standard Kalman filter recursion equation is applied to perform a prediction-update cycle for each measured value, thereby obtaining an optimal estimate of the true power value. The motivation for this is that the actual acquired power signal often contains high-frequency noise introduced by the sensor, power electronic converter or communication link. These noises can interfere with the accurate calculation of subsequent indicators such as deviation, rate of change, etc., and may even cause the controller to malfunction. As an optimal linear state estimator, the Kalman filter can effectively extract the true state of the signal from the noisy measurement sequence, obtaining a smoother and more reliable power curve, thereby providing a higher quality input signal source for the entire control system.

[0159] In the spectrum analysis process of Embodiment 3, in order to suppress spectral leakage, a window function needs to be applied to the signal segment. A preferred enhancement scheme is to use an adaptive window function selection strategy. Specifically, before performing the Fast Fourier Transform (FFT), the system first performs a brief stationarity discrimination on the current signal segment s_i(t) to be analyzed. For example, the ratio of the standard deviation to the mean of the short-time energy can be calculated. If the ratio is greater than a threshold value (e.g. 0.5), it is determined that the signal segment is a non-stationary signal containing AE events, and the Hamming window with lower side lobe peak value is preferred to better capture the main spectrum of transient events. If the ratio is less than the threshold value, it is determined that the signal is relatively stationary background noise, and the Blackman window with wider main lobe and faster side lobe attenuation is preferred to achieve higher side lobe suppression and clearer observation of potential weak resonances. Since different window functions have different trade-offs between main lobe width (affecting frequency resolution) and side lobe suppression (affecting dynamic range). Dynamically selecting the most suitable window function according to the real-time characteristics of the signal can achieve near-optimal spectrum analysis results in different scenarios.

[0160] In the S3b mode of embodiment 8, which compensates for the mechanical resonance frequency, a standing wave can be formed in the membrane when the acoustic field is applied at a single, precise resonance frequency f_membrane for a long time. The standing wave can cause the energy to be fixed in space, resulting in a wave crest with the strongest vibration and a wave node with the weakest vibration, which can make it difficult for the cold spots located in the wave node area to be effectively compensated. To solve this problem, a preferred enhancement scheme is to apply a weak, low-frequency frequency modulation or amplitude modulation to the generated resonance acoustic field signal. For example, when generating the driving signal u_res_i(t)=A_res_i*sin(2πf_membrane*t+φ_avg_i), further multiply it by a low-frequency modulation term [1+m*sin(2πf_mod*t)]; where m is the modulation depth (e.g. 0.1); f_mod is the modulation frequency (e.g. 10 Hz).

[0161] By introducing such a low-frequency perturbation, the frequency or amplitude of the excitation signal is slightly and slowly periodically changed around f_membrane. This change can continuously and slightly change the interference conditions of the acoustic field, effectively breaking the fixed standing wave pattern, and slowly moving the wave crest and wave node of the standing wave in space. The effect is that the acoustic field energy can be more evenly distributed in time average over the entire target compensation area, ensuring that all cold spots can be effectively affected by the acoustic field.

[0162] In the adaptive optimization process of embodiment 9 based on reinforcement learning, when the multi-dimensional parameter lookup table (LUT) has accumulated a large amount of data through long-term learning, a preferred knowledge utilization enhancement scheme is to extract human-readable control rules from it. Specifically, a decision tree algorithm (such as C4.5 or CART algorithm) can be applied, taking the working condition characteristics (η_initial, ΔP, etc.) in the LUT as input features and the corresponding optimal control parameters (K_p*, K_i*, etc.) as output labels to train a decision tree. After training, the path from the root node to each leaf node of the decision tree corresponds to a control rule in the form of IF-THEN. For example, the following rule can be obtained: IFη_initial>0.6 AND ΔPisHigh THEN set K_ff>0.4 AND Acoustic Intensity>8V. Although LUT can achieve the storage and quick call of knowledge, its content itself is a black box, which is not conducive to the understanding, verification and manual intervention of operation and maintenance personnel. By distilling the complex mapping relationship contained in the LUT into a set of clear logical rules, not only can the operation and maintenance personnel be provided with an interpretable system behavior diagnosis basis, but also the optimal strategy learned by the machine can be solidified as part of the expert system, which can play a role even in the simplified version of the system without running the RL algorithm, thereby enhancing the explainability, portability and reliability of the entire technical solution.

[0163] According to one aspect of the present application, in order to solve the problem of insufficient internal state diagnosis capability of proton exchange membrane in the prior art, the present application collects broadband acoustic emission signals generated in the electrochemical reaction process, which have high signal-to-noise ratio, and combines a multi-physical field coupled sound velocity model considering material anisotropy, water content gradient and acoustic dispersion effect, and applies acoustic tomography and super-resolution reconstruction algorithm. The technical path inverses and reconstructs the original acoustic time delay data into a dynamic current density distribution map with millisecond-level refresh rate and millimeter-level spatial resolution. This converts the transient and localized internal electrochemical inhomogeneity from a physical state that cannot be directly observed into an engineering data that can be quantitatively analyzed, thereby providing accurate and high-dimensional information input for subsequent control decisions.

[0164] On this basis, in order to solve the spatial mismatch problem between traditional control means and localized physical problems, the present application proposes a compensation strategy of targeted physical intervention. The strategy uses the high-resolution current density map reconstructed above to accurately identify the low-performance area through density clustering algorithm. Then, an acoustic compensating field is generated by driving the acoustic actuator array, which is spatially focused and energy-adjustable. The acoustic field applies precise physical intervention to the identified area by inducing interface micro-flow to strengthen local mass transfer process or using acoustic elastomeric effect to directly modulate proton mobility through specific physical mechanisms. This method replaces the global control method of traditional reaction total flow or total output voltage with a spatially selective physical adjustment, so that the control energy can directly act on the physical root of the problem, thereby avoiding the inherent response lag and non-targetedness of global adjustment.

[0165] By combining this high-resolution diagnostic capability with targeted physical interventions and integrating it into a multi-timescale coordinated control architecture, this invention achieves significant technical results. Test data shows that this method can reduce the power overshoot of hydrogen power sources from 15% to less than 2%, shorten the response time from 8 seconds to 3 seconds, and improve steady-state power fluctuations from ±5% to ±1%, fundamentally enhancing the dynamic response performance and power quality of hydrogen power sources in microgrids.

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

Claims

1. A method for accurate power regulation and overshoot suppression of a microgrid hydrogen power source, characterized in that, The method comprises the following steps: Collecting acoustic response signals generated by the electrochemical process in the proton exchange membrane inside the hydrogen power supply; Based on the acoustic response signal set, reconstructing the current density in the proton exchange membrane to obtain a spatial distribution map; Analyzing and generating acoustic compensation signals according to the non-uniformity represented by the spatial distribution map, and applying the acoustic compensation signals to the target area of the proton exchange membrane; Integrating the spatial distribution map, the acoustic compensation signal, and the externally input power demand signal and the feedback current output power to generate the final control instruction for adjusting the hydrogen power supply.

2. The method according to claim 1, characterized in that, The acoustic response signal set is collected, including: Identifying and capturing the electrochemical acoustic emission signals emitted by the micro-bubbles formed by the products in the electrochemical process during nucleation and collapse, which have burst and broadband characteristics; The captured electrochemical acoustic emission signals are used as the main data components of the acoustic response signal set.

3. The method according to claim 1, characterized in that, The spatial distribution map is reconstructed by applying an anisotropic multi-physical field coupling sound velocity model, including: The model jointly models the anisotropic material properties of the proton exchange membrane, the water content gradient along the thickness direction of the proton exchange membrane, and the propagation direction and dispersion effect of the acoustic wave; Providing accurate sound velocity benchmarks that dynamically change with spatial position, direction, and frequency for the reconstruction process.

4. The method according to claim 1, characterized by, The spatial distribution map also includes: Processing the acoustic response signal set by acoustic tomography to obtain an initial distribution map limited by the physical diffraction limit; And applying a super-resolution reconstruction algorithm to the initial distribution map to improve the detail recognition ability based on the initial distribution map to generate the final spatial distribution map.

5. The method according to claim 2, characterized in that, Further including: Transforming the captured electrochemical acoustic emission signals into the time-frequency domain and separating them into time-frequency domain sub-band signals; From the time-frequency domain sub-band signals, calculate and extract acoustic emission feature parameters representing each acoustic emission event, including the energy, duration, and occurrence rate of the event; The acoustic emission feature parameters are used as structured input data for the reconstruction step.

6. The method of claim 1, characterized in that, Reconstructing the spatial distribution map includes performing the following operations to obtain an accurate time delay matrix: Calculate the cross-power spectrum between any two sensor signals in the acoustic response signal set, and linearly fit the phase information of the spectrum to obtain high-resolution initial time delay values; Physically consistent verification is performed on the set of initial time delay values, including symmetry verification and triangle inequality verification, to identify and mark abnormal time delay data that do not conform to physical laws; Correct or remove the abnormal time delay data to form an accurate time delay matrix for acoustic tomography calculation.

7. The method of claim 1, characterized in that, One implementation of the step of generating and applying a directional acoustic compensation field includes: Identifying low current density areas from the spatial distribution map; And for the low current density areas, generate and focus acoustic beams working in the medium frequency band; The acoustic beams induce acoustic-induced micro-flow at the interface between the catalyst and the proton exchange membrane in the area, strengthening the local mass transfer process in the area.

8. The method according to claim 7, characterized in that, The step of generating and focusing acoustic beams working in the medium frequency band includes: To achieve focusing, pre-record the acoustic response of the low current density area to obtain the channel information of the medium, and determine a set of phase delay parameters for adaptive focusing based on the time reversal of the channel information; To achieve the generation, according to the compensation demand of the low current density area, a pulse width modulation control sequence is calculated and generated to dynamically adjust the output power of the sound beam; And combine the phase delay parameter with the pulse width modulation control sequence to drive the acoustic actuator array to generate the final intermediate frequency sound beam.

9. The method of claim 1, characterized in that, Generating and applying a directional acoustic compensation field can also be used for: Pre-measuring and obtaining the mechanical resonance frequency of the proton exchange membrane through frequency scanning; Identifying the low current density area from the spatial distribution map; Using the mechanical resonance frequency as the working frequency, a resonance sound field is applied to the low current density area; The resonance sound field uses the acoustic elastance effect to dynamically modulate the proton transition barrier on the proton transport path in the area, directly regulating the proton conductivity of the area.

10. The method according to claim 9, characterized in that, The step of applying the resonance sound field further includes: Constructing a target function for the purpose of maximizing the compensation effect in the low current density area; Solving the target function to obtain an optimized spatial phase distribution; the generated resonance sound field carries the optimized spatial phase distribution and selectively compensates for the low current density area.

Images