Control methods for electrolytic cells to participate in power system frequency regulation ancillary services

Abstract

This invention discloses a control method for electrolyzers participating in power system frequency regulation ancillary services, belonging to the interdisciplinary field of power system control and electrochemical technology. The method includes: collecting grid AGC commands, electrolyzer operating status, and high-frequency impedance signals to generate a synchronization dataset; based on this dataset, extracting millisecond-level electrochemical transient features by analyzing the impedance spectrum and performing interface state estimation; combining the interface state estimate, transient feature vector, and AGC power demand to perform transient sensing multi-objective optimization and determine the optimal power allocation scheme; and based on the optimal scheme and the interface state estimate, performing multi-timescale coordinated control and outputting a power control signal for the electrolyzer. This invention can deeply sense and adapt to the real-time internal state of the electrolyzer, fully releasing its frequency regulation potential while ensuring safety, and improving the performance of frequency regulation services and the stability of equipment operation.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of power system control and electrochemical technology, and in particular, it is a control method for electrolyzers to participate in power system frequency regulation auxiliary services. Background Technology

[0002] The increasing scale of grid-connected renewable energy sources such as wind and solar power poses a challenge to the frequency stability of the power system due to their volatility, necessitating rapid and flexible frequency regulation resources. Proton exchange membrane (PEM) electrolyzers, with their fast response, wide adjustment range, and flexible start-up and shutdown, are ideal units for providing frequency regulation ancillary services. Researching control methods for their participation in frequency regulation is of significant technical importance for improving grid stability and the coupling efficiency of the electricity-hydrogen system.

[0003] For electrolyzers participating in grid frequency regulation, existing technical solutions mainly revolve around two major directions: First, designing to track Automatic Generation Control (AGC) commands, such as using a proportional-integral (PI) controller for feedback adjustment based on power deviation. This involves adjusting the control based on the deviation between the target power issued by the grid and the actual output power of the electrolyzer. During the control process, to ensure the safe operation of the electrolyzer, a set of static, fixed operating constraints are typically set, such as maximum / minimum power limits and maximum power ramp-up rate limits. When the AGC command exceeds these preset static boundaries, the control system will perform amplitude or rate limiting, sacrificing some frequency regulation response to ensure equipment safety. Second, optimizing power response through simplified model predictive control (MPC). This involves predicting future power response using a simplified electrolyzer model to achieve better tracking performance. These methods have, to a certain extent, achieved effective tracking of AGC commands by the electrolyzer, laying the foundation for electrolyzers to participate in grid frequency regulation.

[0004] However, existing control strategies treat the electrolyzer as a black box, ignoring the internal electrochemical dynamics, leading to a mismatch between control and equipment status. Because the internal state (e.g., electrode bubble coverage, ion concentration gradient, etc.) cannot be perceived, control is somewhat blind and requires conservative static constraints, making it difficult to balance safety and rapid response requirements. Therefore, further research is needed to address the problems of existing technologies. Summary of the Invention

[0005] Purpose of the invention: To provide a control method for electrolytic cells to participate in frequency regulation ancillary services of power systems.

[0006] Technical Solution: Collect AGC frequency regulation commands and power demands from the power grid, electrolyzer operating status data, and broadband impedance spectrum data to generate a synchronous dataset; based on the synchronous dataset, extract millisecond-level electrochemical transient features and perform interface state estimation to obtain interface state estimates and transient feature vectors; combine the interface state estimates, transient feature vectors, and AGC power demands to perform transient sensing multi-objective optimization to determine the optimal power allocation scheme; based on the optimal power allocation scheme and interface state estimates, perform multi-timescale coordinated control and output power control signals for the electrolyzer.

[0007] Beneficial effects: This invention solves the problem of existing strategies treating the electrolyzer as a black box and ignoring its internal electrochemical dynamics. By sensing the internal state in real time, it eliminates blind control and matches control with the equipment state. Using the internal state dynamics as a safety constraint instead of static constraints, it balances rapid adjustment and safety, fully releasing the electrolyzer's frequency regulation potential while ensuring safety, thus improving the efficiency and reliability of the electrolyzer's participation in frequency regulation and enhancing its support for power grid stability. Attached Figure Description

[0008] Figure 1 A flowchart illustrating a control method for an electrolytic cell participating in frequency regulation ancillary services of a power system, as provided in an embodiment of this application.

[0009] Figure 2 This is a flowchart illustrating the execution of transient perception multi-target optimization provided in an embodiment of this application.

[0010] Figure 3 The flowchart provided in this application embodiment guides the particle swarm to explore and find a better solution space.

[0011] Figure 4 This is a flowchart illustrating intelligent relocation of specific particles in a particle swarm, provided as an embodiment of this application.

[0012] Figure 5 A flowchart illustrating the control parameters for dynamically adjusting particle swarm optimization provided in this application embodiment. Detailed Implementation

[0013] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0014] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0015] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0016] Existing technical solutions generally simplify the electrolyzer into a black-box model with fixed response characteristics. Their control strategies primarily rely on external power grid commands, severely neglecting the rapidly changing and complex electrochemical dynamics within the electrolyzer. This leads to a deep-seated mismatch between the control strategy and the physical state of the equipment. This mismatch results in an irreconcilable contradiction between pursuing rapid response and ensuring operational safety, specifically manifested in the following two aspects:

[0017] First, existing control strategies are inherently blind due to their inability to perceive internal states, resulting in compromised safety and efficiency during frequency modulation. The actual adjustable capability of an electrolyzer is not static but is severely constrained by its internal millisecond-level electrochemical interface states (such as the bubble coverage on the electrode surface and the concentration gradient of ions within the diffusion layer). For example, if a large number of bubbles have already accumulated on the electrode surface when the current density suddenly increases, forcibly increasing the power will lead to a sharp increase in ohmic losses and overpotential, causing voltage overshoot or even localized overheating, endangering the lifespan of the electrolyzer's membrane electrodes. Existing blind control methods are completely unaware of this and will still execute power adjustments according to external commands, pushing the equipment to the edge of potential risks. This disconnect between control and state forces the equipment to maintain a large safety margin at all times, preventing it from fully utilizing its performance.

[0018] Secondly, to adapt to all possible unknown operating conditions, existing control strategies must adopt a one-size-fits-all static and conservative design, suppressing the electrolyzer's rapid response potential and reducing the value of its frequency regulation services. Because it is impossible to assess the electrolyzer's internal health status in real time and accurately, the controller can only set global, static operating constraints (such as a very conservative power ramp-up rate) based on the worst-case scenario. This means that even if the electrolyzer is in excellent internal condition at a certain moment and is fully capable of providing a faster and greater power response, it will be limited by this static constraint, resulting in the waste of its valuable rapid adjustment potential. This conservative strategy causes the electrolyzer to operate at levels far below its actual adjustable capabilities most of the time, significantly reducing the quality of its frequency regulation services and weakening its competitiveness in the ancillary services market.

[0019] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0020] Example 1: A control method for electrolyzers participating in power system frequency regulation ancillary services is provided, specifically including:

[0021] Step S1: Collect AGC frequency regulation commands and power demand from the power grid, electrolyzer operating status data, and broadband impedance spectrum data to generate a synchronization dataset, including:

[0022] The system receives AGC command data packets from the power grid dispatch center via industrial communication protocols (such as IEC 61850 and IEC 104), parses these packets, and extracts the target power value PAGC(t) and the response time requirement t. _response In addition to information such as frequency modulation priority level, the accuracy and validity of the instructions are further ensured through CRC check and time stamp verification, and a standardized AGC instruction sequence is output.

[0023] The operating status data of the electrolytic cell is acquired at a sampling rate of not less than 1 kHz. This data constitutes an operating status data matrix, which includes at least the terminal voltage U of the electrolytic cell. cell (t), Operating current I cell (t), electrolyte temperature T _elec (t) and import / export pressure P _in / out(t). The acquired analog signal is digitized by an analog-to-digital converter (ADC) with at least 12-bit precision, and high-frequency noise can be eliminated using methods such as moving average filtering. Furthermore, a micro-perturbation AC signal is injected into the electrolytic cell. This signal has a wide frequency range (e.g., from 0.1 Hz to 100 kHz), and its amplitude is limited to a range that does not affect the normal operation of the electrolytic cell (e.g., not exceeding 5% of the rated current). The voltage and current response of the electrolytic cell are synchronously detected by devices such as a lock-in amplifier, and the complex impedance Z(ω) = U(ω) / I(ω) varying with frequency is calculated in real time and updated at a frequency of, for example, once every 10 milliseconds, forming broadband impedance spectrum data Z. spectrum (t, ω).

[0024] Using a high-precision clock source (such as a GPS clock) as a unified benchmark, the previously acquired AGC command sequence, operational status data matrix, and broadband impedance spectrum data are rigorously time-aligned. Missing data is imputed using interpolation algorithms; outliers are marked. This step aims to eliminate time asynchrony issues between multi-source data, ensuring the accuracy of subsequent calculations and outputting a synchronized dataset D containing all synchronization information. sync (t).

[0025] Step S2: Based on the synchronous dataset, extract millisecond-level electrochemical transient features and perform interface state estimation to obtain interface state estimates and transient feature vectors.

[0026] In this embodiment, the synchronous dataset acquired in the previous step, particularly the broadband impedance spectroscopy data, is used to quantify internal states that cannot be directly measured but are crucial to the performance and safety of the electrolyzer, such as the bubble coverage on the electrode surface and the ion concentration distribution in the electrolyte. An advanced state estimation algorithm (such as an extended Kalman filter) is constructed and applied to fuse various pieces of information to obtain an accurate depiction of these internal states. Optionally, the implementation details of this process will be described in detail in Embodiment Two.

[0027] Step S3: Combine the interface state estimate, transient feature vector, and AGC power requirement to perform transient sensing multi-objective optimization and determine the optimal power allocation scheme.

[0028] This step no longer simply responds to AGC power commands, but treats AGC requirements as an optimization objective, comprehensively considering the internal electrochemical state of the electrolyzer (characterized by the output of step S2). A perceptual multi-objective optimization algorithm (e.g., an improved adaptive particle swarm optimization algorithm) is run to meet grid frequency regulation requirements (such as frequency regulation accuracy and response speed) while also considering the electrolyzer's own operating efficiency and safety constraints. This algorithm can dynamically adjust its search strategy to adapt to the complex nonlinear behavior of the electrolyzer under different operating conditions. The implementation details of this process will be elaborated in Examples 3 and 4.

[0029] Step S4: Based on the optimal power allocation scheme and the interface state estimate, perform multi-time-scale coordinated control and output a power control signal for the electrolytic cell.

[0030] In this embodiment, this step receives the optimal power allocation scheme and, combined with real-time insight into transient disturbances within the electrolyzer, generates the final control signal through a carefully designed multi-timescale coordinated control architecture. This architecture can simultaneously perform coordinated control at multiple timescales, including milliseconds, hundreds of milliseconds, and seconds, enabling rapid compensation for internal electrochemical disturbances and precise tracking of the optimized power target, ensuring the electrolyzer provides fast and accurate frequency modulation services. The implementation details of this process will be elaborated in Embodiment 5.

[0031] Example 2 describes an optimal implementation method for extracting millisecond-level electrochemical transient features and performing interface state estimation based on a synchronous dataset to obtain the interface state estimate and transient feature vector. Specifically:

[0032] Step S21: Extracting millisecond-level electrochemical transient features and performing interface state estimation, including: analyzing broadband impedance spectral data in the synchronous dataset, decoupling and obtaining preliminary quantization values ​​characterizing bubble coverage and ion concentration gradient; constructing and applying an extended Kalman filter, receiving the preliminary quantization values ​​as observations; obtaining an accurate depiction of the electrochemical interface state of the electrolyzer through recursive estimation of the extended Kalman filter, and outputting interface state estimates; constructing and outputting transient feature vectors based on the interface state estimates; wherein, decoupling and obtaining preliminary quantization values ​​includes: dividing the broadband impedance spectral data into high-frequency and low-frequency bands in the frequency domain; for the high-frequency impedance data reflecting the rapid dynamics of the electrode interface, resolving parameters related to the effective electrode area, and calculating the preliminary quantization value of bubble coverage; for the low-frequency impedance data characterizing the slow ion diffusion process in the electrolyte, calculating parameters related to mass transfer limitation, and converting them into the preliminary quantization value of ion concentration gradient.

[0033] The acquired broadband impedance spectrum data Z _spectrum(t, ω) can be preprocessed, and a Savitzky-Golay filter can be applied to smooth the noise, dividing it into different frequency bands. For example, the high-frequency band can be set to 10kHz to 100kHz, and the low-frequency band can be set to 0.1Hz to 1kHz.

[0034] The specific implementation path for estimating bubble coverage includes: fitting high-frequency impedance data with a preset equivalent circuit model of the electrode interface to identify and extract real-time double-layer capacitance parameters; establishing a conversion relationship based on the physical model, comparing the real-time double-layer capacitance parameters with the pre-calibrated reference double-layer capacitance in a bubble-free state, and calculating the preliminary quantitative value of bubble coverage.

[0035] The equivalent circuit model can be a simplified Randle model, whose expression is:

[0036] Z eq =R s +1 / (jωC dl +1 / R ct 1) By fitting the high-frequency measurement data to the model using the complex nonlinear least squares method, the double-layer capacitance C can be solved in real time. dl (t). Since the size of the double-layer capacitance is proportional to the effective surface area of ​​the electrode, the attachment of bubbles will reduce the effective area. Therefore, the following relationship can be established to calculate the bubble coverage θ(t): θ(t) = 1 − (Cd l，0 C dl (t)) 1 / β Among them, C _dl，0 The reference double-layer capacitance value is measured in the bubble-free state after factory calibration or recent maintenance of the electrolytic cell and is used as a pre-stored parameter; β is the bubble shape factor, used to correct for errors caused by bubbles not being ideally covered by two-dimensional discs. Its value is usually set empirically, for example, 0.85 for spherical bubbles. Through this method, a preliminary quantitative value of bubble coverage can be obtained.

[0037] Based on this, the low-frequency impedance data is mainly affected by the diffusion process of ions in the electrolyte, which can be characterized by Warburg impedance. By fitting the data in the 0.1-1kHz frequency band to Warburg impedance, i.e., fitting Zw = σ / sqrt(w) × (1−j), the Warburg coefficient σ(t) can be extracted. This coefficient directly reflects the magnitude of the mass transfer resistance. According to Fick's diffusion law, the ion concentration gradient ∇c(t) can be calculated from parameters such as the diffusion coefficient and current density, and the diffusion coefficient is inversely proportional to the square of the Warburg coefficient σ(t). Therefore, by analyzing σ(t), a preliminary quantitative value of the ion concentration gradient ∇c(t) can be calculated.

[0038] This step utilizes dual-frequency impedance decoupling technology to simultaneously separate bubble coverage and ion concentration gradient information from a single broadband impedance spectrum. This method solves the problem of traditional methods being unable to obtain key electrochemical states inside the electrolyzer in real time and non-invasively. It leverages the principle that different physical processes (interfacial dynamics and diffusion dynamics) have different response speeds in the frequency domain, achieving multi-purpose use of a single source. Compared to relying on multiple expensive, complex, or delayed sensors, this invention provides a low-cost, fast-response (millisecond-level), and information-rich online diagnostic method. This not only provides unprecedentedly refined input for upper-level control, solving the problems of unclear or inaccurate readings, but also enables real-time monitoring of bubble coverage and ion concentration—two key indicators directly affecting electrolyzer efficiency and lifespan. This makes predictive maintenance and health management possible, indirectly extending the lifespan of expensive electrolyzer equipment while improving frequency tuning performance.

[0039] Step S22: Construct and apply an extended Kalman filter, and receive the preliminary quantization values ​​as observations;

[0040] A nonlinear state-space model is constructed to describe the interface dynamics of the electrolyzer, and its state vector x(k) contains at least four core state variables: x(k) = [θ(k), ∇c(k), dθ / dt(k), d∇c / dt(k)]. ^T , representing bubble coverage, ion concentration gradient, rate of change of bubble coverage, and rate of change of ion concentration gradient, respectively; the bubble coverage θ_ obtained by decoupling in step S21 meas(k) and ion concentration gradient ∇c_ meas(k) is then used as an observation of EKF, i.e., z(k) = [θ _meas (k),∇c _meas (k)] ^T .

[0041] The extended Kalman filter has a dynamically constructed process noise covariance matrix Q. A transient severity index characterizing the stability of the system (see Example 4) is introduced as an adjustment input. Based on the real-time value of this transient severity index, the amplitude of the basic process noise covariance matrix is ​​dynamically adjusted. This allows the confidence level of the filter in its internal state model to change inversely with the severity of the system transient, ensuring reliable estimation results under varying operating conditions.

[0042] Specifically, the process noise covariance matrix Q(k) is calculated as follows: Q(k) = Q base ×(1+K adapt ×S trans ); where Q _baseIt is a preset basic process noise covariance matrix, the value of which is set based on prior knowledge of the model uncertainty under normal system operating conditions, such as Q. _base =diag([0.001, 0.01, 0.0001, 0.001]). S _Trans The transient severity index, K, is calculated in subsequent steps (see Example 4 for details). This index quantifies the severity of the system's deviation from steady state; a larger value indicates a more unstable system. _adapt It is the adaptive gain coefficient, which can be set to 5 for example.

[0043] In this embodiment, when the system is running smoothly, S _Trans When the value is close to 0, Q(k) is approximately equal to Q. _base This indicates that the filter highly trusts its internal state transition model; when the system experiences severe transients (such as large power fluctuations causing drastic changes in bubble behavior), S _Trans Increasing the value of Q(k) increases the diagonal elements of Q(k). This implies that the current model has significant uncertainty, and when updating the state, the filter needs to increase the weights that depend on the new observations and decrease the dependence on the model's predictions. This adaptive adjustment mechanism enhances the filter's robustness and tracking accuracy under drastic changes in operating conditions.

[0044] Through the recursive prediction and update steps of EKF, an accurate description of the electrochemical interface state can be output, namely the interface state estimate x. est (t), and a 12-dimensional transient feature vector F that integrates the estimated value and other transient information. trans This provides a basis for decision-making in subsequent optimization control.

[0045] Step S23: Obtain an accurate description of the electrochemical interface state of the electrolyzer through recursive estimation using an extended Kalman filter, and output the interface state estimate.

[0046] In the embodiments of this application, this step achieves real-time monitoring and prediction of the internal state of the electrolyzer by fusing multi-source data. Specifically, a state-space model is constructed, with a state vector x = [θ, ∇c, dθ / dt, d∇c / dt]ᵀ, and the bubble coverage θ(t), ion concentration gradient ∇c(t), and transient event feature matrix E are read. _Trans The process noise covariance Q is recursively estimated using the nonlinear state transition equation x(k+1)=f(x(k), u(k))+w(k) and the observation equation z(k)=h(x(k))+v(k). The process noise covariance Q is adaptively adjusted according to the severity of the transient. The output interface state estimate is x. _est The equations P(t) and the estimation error covariance P(t) provide high-precision state information for subsequent model-based optimization control.

[0047] This step introduces a transient severity index to dynamically adjust the process noise covariance matrix Q of the EKF filter, solving the problem of traditional state estimation algorithms experiencing divergence or severe accuracy degradation due to model mismatch when facing sudden changes in system operating conditions. The Q matrix of a traditional EKF is fixed, representing its unchanging confidence in the system model. This invention makes this confidence intelligent: when the system is stable, it places greater trust in the model, ensuring the smoothness of the estimation; when the system undergoes drastic changes, it increases the Q value to reduce its confidence in the model, relying more on new measurement data, thus enhancing the state estimator's ability to track external disturbances and internal sudden changes and its robustness. This ensures that regardless of the electrolyzer's operating conditions, the subsequent optimization control module can obtain high-precision, high-reliability internal state data, which is the fundamental premise and guarantee for the successful implementation of the entire transient sensing control strategy.

[0048] Step S24: Based on the interface state estimation values, construct and output the transient feature vector;

[0049] Integrated interface state estimate x _est (t), transient event feature matrix E _Trans A 12-dimensional transient feature vector F is constructed from the temperature change rate dT / dt in the operating status data. _Trans =[θ,∇c,dθ / dt,d∇c / dt,ΔZ,τ _Trans The vectors dT / dt, ... are subjected to Min-Max normalization to ensure that each dimension is in the interval [0, 1], resulting in the normalized transient feature vector F. _norm This eliminates the influence of the dimensions of different features on subsequent analysis.

[0050] Example 3 describes the optimal implementation of the joint interface state estimate, transient feature vector, and AGC power requirement, and performs transient sensing multi-objective optimization to determine the optimal power allocation scheme. It also elaborates on the intelligent initialization and dynamic parameter adaptation steps in the transient sensing adaptive particle swarm optimization (PSO) method.

[0051] In this embodiment, performing transient sensing multi-objective optimization includes:

[0052] Step S31: Construct a directional initial particle swarm using transient feature vectors;

[0053] Constructing a guided initial particle swarm includes: parsing and transforming the system deviation information contained in the transient feature vector into transient offset predictions that characterize the drift direction and amplitude of the optimal region in the solution space; based on previously found historical optimal solutions and guided by the transient offset predictions, generating a core particle swarm within a target region with potential greater than a threshold; randomly distributing an exploration particle swarm within the algorithm's feasible region; and merging the core particle swarm and the exploration particle swarm to form the initial particle swarm. Specifically, this process can be decomposed into the following sub-steps:

[0054] The 12-dimensional transient eigenvector F norm Through the pre-trained mapping matrix M map (For example, a 12×3 matrix), mapped to a 3D search space corresponding to the optimization variables (such as power P, current I, voltage V), to obtain the feature mapping coordinates x. feature =M map ×F norm The mapping matrix M _map It can be obtained by dimensionality reduction methods such as principal component analysis (PCA) on historical data, aiming to retain, for example, more than 95% of the information in the original features.

[0055] Calculate the transient migration prediction Δx _Trans The prediction incorporates two types of information: one is the offset x, which is directly driven by the current system transient. feature Another type is the demand offset Δx generated to meet grid commands. demand Δx _demand It can be determined by the power deviation ΔP _demand =P _AGC -P _current The system's sensitivity matrix (representing the effect of changes in control variables on power) is calculated. The fusion formula can be: Δx trans =α trans ×x feature +(1−α trans Δx demand ; where the transient correction coefficient α _Trans With transient severity S _Trans Positive correlation, such as α _Trans =0.7×(1-exp(-2×S _Trans This means that the more severe the system transients, the greater the S... _Trans The larger α is _Trans The closer it is to 0.7, the more the direction of initialization depends on the direct guidance of transient characteristics.

[0056] Generate the initial particle swarm: Assume a total number of particles N=50, of which 80% (40) are core particles whose positions are at the previous optimal solution x. _best\_Prev Based on this, along Δx_Trans Generate by small-scale perturbation in the direction: x _core\_i =x _best\_Prev +Δx _Trans +ε _i , where ε _i The first 20% are normally distributed random numbers with a mean of 0 and a small standard deviation (e.g., 10% of the search step size). The remaining 20% ​​(10) are exploration particles, randomly generated within the entire feasible region: x explore_j =x min +rand()×(x max -x min ).

[0057] This step utilizes real-time transient eigenvectors to predict the offset direction of the optimal solution and constructs a directional initial particle swarm accordingly, overcoming the inherent defects of traditional PSO algorithms, such as slow convergence speed and susceptibility to local optima due to random initialization. This invention uses transient information before optimization begins to pre-identify high-potential regions where the optimal solution is most likely to exist, concentrating most computational resources there. This allows the optimization process to start from a point far superior to a random state, shortening the number of iterations and computation time required for the algorithm to converge to the global optimum. In power frequency regulation scenarios with extremely stringent response speed requirements, it can output high-quality power allocation schemes faster, improving the overall timeliness of frequency regulation response.

[0058] Step S32: During the iteration process, the control parameters for particle swarm optimization are dynamically adjusted based on the electrochemical environment reflected by the interface state estimation.

[0059] The control parameters for dynamically adjusting particle swarm optimization include: quantifying the impact of temperature on operating efficiency, the impact of bubble dynamics on control stability, and the limitation of ion concentration gradient on mass transfer capacity from electrolyzer operating status data and interface state estimates; constructing factors characterizing the multi-physics state based on these factors (integrating multi-physics interaction characteristics), including at least thermal deviation factor, bubble dynamics factor, and concentration gradient factor; establishing and utilizing a preset mapping relationship to couple this set of multi-physics factors into a real-time adjustment basis for particle swarm optimization control parameters, generating an adaptive PSO control parameter set.

[0060] The process of generating the adaptive PSO control parameter set specifically includes: multiplying the thermal deviation factor, the bubble dynamic factor, and the concentration gradient factor to form a comprehensive adjustment coefficient; and using the comprehensive adjustment coefficient to modulate the preset basic inertia weight to generate a dynamic inertia weight.

[0061] Specifically, the formula for calculating the dynamic inertia weight ω(t) is: ω(t) = ω base ×f thermal ×f bubble ×f conc; where ω _base This is the basic inertia weight, for example, 0.9. This multiplicative coupling mechanism means that any substantial deterioration in the physical field (corresponding to a smaller factor) will directly reduce the total inertia weight, slowing down the particle swarm's flight speed and making the search behavior more cautious. This effectively avoids system oscillations or deviations from the feasible region due to excessively large search steps under harsh conditions. Furthermore, the learning factors c1(t) and c2(t) can also be coupled with these factors, for example: c1(t) = c1 base ×(2−f bubble c2(t)=c2 base ×f conc This makes it possible when the bubble is dynamically stable (f _bubble When the concentration gradient is close to 1, particles are more likely to believe in their own historical best (c1 is larger); when the concentration gradient is limited (f _conc (smaller), particles rely more on the guidance of the global optimum (c2 is larger than c1).

[0062] This step addresses the fundamental challenge of adapting to the complex and variable electrochemical environment of the electrolyzer by mapping multiple physical field states within the electrolyzer, such as thermal deviation, bubble dynamics, and concentration gradient, to core control parameters like the inertia weight of the PSO algorithm in real time via a multiplicative coupling mechanism. For example, when bubble dynamics are volatile (bubble dynamics factor decreases), the inertia weight is adjusted accordingly, making the algorithm more cautious, reducing blind exploration, avoiding power oscillations, and enhancing control stability and robustness. Similarly, when temperature drift is severe (thermal deviation factor decreases), the algorithm adjusts its search strategy to prioritize ensuring the system returns to the optimal operating temperature range, effectively avoiding equipment damage risks caused by efficiency degradation or thermal runaway. This design, which deeply binds the algorithm's characteristics to physical reality, ensures that the optimization process always operates within the safe and efficient boundaries of the electrolyzer, achieving a dynamic balance between frequency modulation performance and long-term healthy operation of the equipment.

[0063] Step S33: For each particle in the particle swarm, in combination with the AGC power requirements, evaluate its (each particle in the particle swarm) comprehensive fitness in terms of frequency modulation accuracy, response speed and operating efficiency, and obtain the corresponding fitness value.

[0064] In a specific embodiment, this step quantifies the performance of particles under different control objectives through multidimensional performance indicators, and constructs an evaluation system that comprehensively reflects the frequency modulation capability of the electrolyzer by combining dynamic constraints; the fitness value, as the optimization objective of the PSO algorithm, directly affects the search direction and convergence speed of the particles.

[0065] Step S34: Based on the fitness value and dynamically adjusted control parameters, guide the particle swarm to explore and find a better solution space until the convergence condition is met, forming the global optimal solution, which is the optimal power allocation scheme.

[0066] Example 4 illustrates the multimodal search and intelligent particle relocation strategy of the transient sensing adaptive PSO method in Example 3 during the iterative process.

[0067] In a specific embodiment, guiding the particle swarm to explore and optimize towards a better solution space includes: extracting multi-dimensional transient information, including impedance change rate, temperature change rate, and power second derivative, from the electrolytic cell operating status data and transient feature vectors; generating a transient severity index characterizing the current stability of the system by weighted aggregation and quantification of the multi-dimensional transient information; and switching the particle swarm's optimization strategy between preset progressive search, extended search, and jump search modes based on the value of the transient severity index; wherein, when the optimization strategy is switched to the jump search mode, guiding the particle swarm to explore and optimize towards a better solution space includes intelligent relocation of specific particles in the particle swarm, including:

[0068] Choose one of the following strategies to generate a jump vector for relocation: Use the gradient of the objective function of the current optimization problem as the direction, and use the real-time impedance changes recorded in the transient event feature matrix to calibrate the jump amplitude and construct a transient relocation vector; perform pattern matching between the current operating condition features composed of interface state estimates and a pre-configured knowledge base; if similar historical operating conditions are found, refer to their corresponding adjustment strategies to construct an experience-guided jump vector; update the particle's position based on the jump vector.

[0069] Optionally, the intelligent relocation of specific particles in the particle swarm includes:

[0070] Step S41: Generate transient severity index;

[0071] In a specific embodiment, the process includes: normalizing information from different physical dimensions, such as impedance change rate, temperature change rate, and power second derivative, to eliminate dimensional differences and obtain a set of independent transient indices; performing a preset weighted summation on the set of independent transient indices to merge them into a single-dimensional transient severity index that can comprehensively reflect electrical, thermal, and power stability.

[0072] Specifically, the transient indices for each dimension are calculated:

[0073] Impedance transient index S Z :S Z =∣dZ / dt∣ / Z max Where |dZ / dt| is the real-time monitored rate of impedance change, Z _max The preset maximum allowable rate of change threshold, for example, 0.5 Ω / s; thermal transient index S T :S T =∣dT / dt∣ / T maxWhere |dT / dt| is the rate of change of electrolyte temperature, T _max A preset threshold, such as 0.5°C / s; power transient index S P :S P =∣d 2 P / dt 2 | / P max , where |d 2 P / dt 2 | is the second derivative of the actual power of the electrolyzer, reflecting the acceleration of the power change, P _max For a preset threshold, such as 50kW / s 2 ; Obtain the transient index vector S of each dimension _vec =[S _Z S _T S _P The comprehensive transient severity index S is obtained through weighted aggregation. trans :S trans =w Z ×S Z +w T ×S T +w P ×S P ;

[0074] Where, weight w = , w T w P Preset the parameters based on the importance of different transient effects on the system, for example, w Z =0.4, w T =0.3, w P =0.3, and w Z +w T +w P =1. S _Trans It is a dimensionless value between 0 and 1. The larger the value, the more severe the transient the system is currently experiencing and the worse its stability.

[0075] Step S41: Switch the search mode based on the transient severity index;

[0076] Based on the calculated S trans The algorithm can adaptively switch between multiple search modes:

[0077] When S trans When the value is less than 0.3, the system is stable. Switch to incremental search mode and set a smaller search step size. _size =0.05×search_range, perform a refined search around the current optimal solution; when 0.3≤S trans When the value is less than 0.7, the system exhibits a certain degree of disturbance. Switch to extended search mode and appropriately increase the set search step size._size =0.15×search_range, expanding the search range; when S trans When the value is ≥0.7, the system experiences a severe transient and switches to jump search mode. At this point, the algorithm considers the current optimal region to be invalid and needs to perform a large-scale jump search to find a new optimal region. The step size is set to step. _size =0.4×search_range; Calculate the directional weight w _dir =1+2×S _Trans Output the search mode identifier (mode) and the adaptive step size (step). _size and direction weight w _dir .

[0078] Step S43: When the optimization strategy switches to jump search mode, guide the particle swarm to explore and optimize towards a better solution space, including:

[0079] Intelligent relocation is performed on specific particles in the particle swarm; a strategy is selected from steps S44 and S45 to generate a jump vector for relocation; the particle's position is updated based on the jump vector; specifically, when the search mode identifier mode is jump search, the transient event feature matrix E is read. _Trans Z in _new and Z _old Calculate the impedance change ΔZ=Z _new -Z _old It directly reflects the physical strength of the disturbance; the gradient of the decision variables of the current objective function ∇J=[ΨJ / ΨP, ΨJ / ΨI, ΨJ / ΨV] is estimated by the finite difference method, which reflects the optimal improvement direction; the direction vector is calculated. _vector =∇J / |∇J|;Determine the adaptive gain K _Z =0.5×(1+S _Trans )×sign(ΔZ), whose value is positively correlated with the severity of the transient, is used to amplify the amplitude of the jump; generate the relocation vector Δx _Transient =K _Z ×|ΔZ|×direction _vector ×Direction weight w _dir Output transient relocation vector Δx _Transient .

[0080] The algorithm selects one of the strategies (e.g., prioritizing strategy two, and using strategy one if a match fails) to generate a jump vector to update the position of a specific particle: x new =x current +x jumpThis mechanism enables the PSO algorithm to perform intelligent relocation when faced with severe disturbances. Instead of simple random mutations, it adopts large jumps with clear objectives and physical or empirical justifications, enhancing the algorithm's global search capability and robustness in the face of unexpected events.

[0081] Step S44, Strategy 1: The current operating condition features, composed of interface state estimates, are matched with a pre-configured knowledge base. If similar historical operating conditions are found, the adjustment strategies corresponding to these (similar historical operating conditions) are used to construct an experience-guided jump vector. Specifically, the current operating condition (including the transient index vector S) is... _vec and interface state estimate x _est ) constitute the eigenvector F _current =[S _vec ,θ,∇c]; in the pre-stored historical knowledge base KB = (F _hist ΔP _opt Search within the range ... to calculate F. _current With each historical operating condition feature F in the database _ hist \ _i similarity score _i =exp(-||F _current -F _hist_i || 2 / σ 2 ), where σ=0.2; if the highest similarity max(score) _i If the value exceeds a threshold (e.g., 0.8), a match is successful, and the optimal power adjustment ΔP from the corresponding historical cases is extracted. _historical Generate Lévy flight stride: Lévy _step =step _size ×tan(π×(rand()-0.5))×direction _vector This is used to increase randomness to escape local optima; the empirically guided jump amount x is calculated. _jump_exp =x jump_exp =whist×ΔP historical +(1−whist)×Leˊvy _step =0.8×ΔP _historical +0.2×Lévy _step , where w _hist This represents the weight of historical experience (e.g., 0.8); the output is the experience jump vector x. _jump_exp Or match failure flag.

[0082] This step leverages historical experience to guide current decision-making. If the system has successfully handled similar crises in the past, directly referencing successful strategies from that time can lead to a faster and more accurate response.

[0083] Step S45, Strategy 2: Using the gradient of the objective function of the current optimization problem as the direction, the transition amplitude is calibrated using the real-time impedance change recorded in the transient event feature matrix, and a transient relocation vector is constructed; specifically, the particle position matrix X and the transient relocation vector Δx are read. _Transient and the experience jump vector x _jump_exp (If it exists), for each particle i that needs relocation, perform the following: if the match is successful, then x _new_i =x _current_i +x _jump_exp If the match fails, then x _new_i =x _current_i +Δx _Transient Check boundary constraints; if x _new_i Exceeding [x _min x _max ], execute intelligent reflection x _reflected_i =x _boundary -0.5×(x _new_i -x _boundary )×exp(-S _Trans This ensures that the particles are within the feasible region; it updates the particle positions and marks the relocated particle indices, outputting the relocated particle position matrix X. _relocated and relocation particle index set I _relocated .

[0084] Step S46: Search performance evaluation and knowledge base update;

[0085] Calculate the average fitness improvement rate η before and after relocation _improve =(f _after -f _before ) / f _before If η _improve If the similarity is greater than 0.1 and the current situation is a new one (similarity < 0.8), then the current feature vector F... _current and actual power adjustment ΔP _actual Add to the historical knowledge base; analyze the usage frequency and success rate of each search mode, update the mode switching threshold; output the search performance metric η. _improve and the updated knowledge base KB _updated .

[0086] This step constructs a transient severity index and adaptively switches the search mode accordingly (especially triggering jump search under severe disturbances). This solves the problem of response lag and tracking failure caused by traditional optimization algorithms (such as standard PSO) when facing sudden changes in power grid commands or drastic changes in the internal operating conditions of electrolyzers, due to excessively small search step size and conservative optimization direction. When the system is stable, a fine-grained asymptotic search is used to ensure optimization accuracy; once a drastic change is detected, the algorithm can immediately abandon local fine-tuning and instantly transfer computing resources (particles) to a more promising solution space region through intelligent relocation based on physical mutations (impedance changes) or historical experience. This shortens the system's confusion period under disturbances, achieves millisecond-level fast decision response and accurate and rapid tracking of high-power level commands, and improves the response speed and reliability of frequency regulation services. Furthermore, this ability to avoid wasting computing resources in ineffective regions also improves the overall computational efficiency of the control system, allowing high-quality optimization solutions to be completed within a shorter time window.

[0087] Example 5 describes an optimal implementation method for multi-timescale coordinated control based on the optimal power allocation scheme and interface state estimates, outputting power control signals for the electrolyzer. It also elucidates the implementation architecture and mechanism of multi-timescale coordinated control. Specifically:

[0088] Step S51: Fast control layer, based on the interface state estimate, generates feedforward compensation control quantity aimed at offsetting transient electrochemical disturbances;

[0089] This internal compensation layer actively counteracts rapid disturbances within the electrolyzer caused by electrochemical processes (such as sudden bubble aggregation or desorption, or local ion concentration depletion), preventing these internal disturbances from propagating to the outside and affecting the stability of power output. Its update cycle is extremely short, for example, 10 milliseconds.

[0090] Specifically, the fast control layer generates feedforward compensation control quantities, and the specific steps include:

[0091] Based on the bubble coverage change rate in the interface state estimate, the bubble compensation power is estimated to offset the power fluctuations caused by bubble aggregation or desorption; based on the second derivative of the ion concentration gradient in the interface state estimate, the concentration gradient compensation power is estimated to alleviate the local concentration polarization effect; the bubble compensation power and the concentration gradient compensation power are superimposed to form the final feedforward compensation control quantity.

[0092] Among them, the bubble compensation power ΔP _bubble The calculation aims to compensate for the increase or decrease in effective electrode area and the fluctuation in current density caused by changes in bubble coverage θ(t). For example, it can be calculated based on the following model: ΔP _bubble =-K _b1 • dθ / dt•V _cell •I_rated -K _b2 •Δi•V _cell In this formula, dθ / dt is derived from the interface state estimate x. _est The rate of change in bubble coverage obtained in (t); V _cell I is the terminal voltage of the electrolytic cell; _rated K is the rated current; Δi is the instantaneous current density increment caused by bubble coverage; K _b1 and K _b2 These are the fast response compensation coefficient and the steady-state compensation coefficient, respectively, which were calibrated experimentally.

[0093] Concentration gradient compensation power ΔP _conc The calculations aim to compensate for the increase in overpotential and decrease in efficiency caused by limited ion mass transfer (concentration polarization). For example, it can be calculated based on the following model: ΔP _conc =-K _c1 •∇ 2 c•A _electrode •i _lim In this formula, ∇ 2 c is the interface state estimate x _est The spatial second derivative of the ion concentration gradient obtained in (t) characterizes the non-uniformity of the concentration distribution; A _electrode i represents the electrode area; _lim K represents the local limiting mass transfer current density. _c1 To compensate for the gain coefficient, the two are superimposed and then filtered and clipped to obtain the output of the fast layer: u fast (t)=ΔP comp (t)=ΔP bubbl e+ΔP conc .

[0094] Step S52, Medium-speed control layer: Based on the deviation between the optimal power allocation scheme and the actual power of the electrolyzer, a feedback regulation control quantity is generated to achieve precise tracking of the target power.

[0095] This layer precisely executes the optimal power allocation scheme P opt (t); its update cycle is moderate, for example, 100 milliseconds; this layer typically uses a PI (proportional-integral) controller with feedforward, whose output u _medium (t) is calculated as follows: u medium (t)=K _P ⋅ e(t)+K _i ⋅∫e(t)dt+K _ff ⋅dP _AGC / dt; where the power deviation e(t) = P _opt (t)-P _actual (t), P_actual (t) is the actual output power of the electrolytic cell; K _P and K _i These are the proportional and integral gains of the PI controller, which can be optimized by PSO or pre-tuned; K _ ff It is the feedforward gain, dP AGC / dt is the rate of change of the power grid AGC command. This feedforward term can respond to changes in the command in advance, improving the tracking speed. To prevent integral saturation, upper and lower limits are set for the integral term.

[0096] Step S53: Slow control layer, based on the evolution trend of historical temperature data, generates preventive management control quantities aimed at maintaining the long-term thermal balance of the system;

[0097] This layer implements thermal management of the electrolyzer. Due to its high thermal inertia, its response and decision-making occur on timescales of seconds or longer; its update cycle is relatively long, for example, 1 second; this layer analyzes, for example, the temperature data sequence T from the past 60 seconds. history Predicting future temperature trends T pred If the predicted temperature will exceed the safe or efficient operating range (e.g., [T] _min T _max Then, a preventative power adjustment amount u is generated. _slow (t). For example, if the predicted temperature T _pred Will exceed T _max , then u slow (t)=−Kt⋅(T pred -T max ), where Kt is the thermal management adjustment coefficient, which suppresses further temperature rise by lowering the power reference.

[0098] Step S54: The three layers of control quantities are synergistically integrated to form a comprehensive control quantity, which serves as the basis for generating the electrolytic cell power control signal.

[0099] Specifically, the collaborative fusion to form a comprehensive control quantity includes: introducing a transient severity index characterizing the stability of the system as a decision-making basis; dynamically allocating weights to the control quantities output by the three control layers (fast, medium, and slow) based on the real-time value of the index; increasing the weight ratio of the fast control layer when the index value increases, and increasing the weight ratio of the medium control layer when the index value decreases; and performing algebraic summation on the dynamically weighted three-layer control quantities to obtain the comprehensive control quantity.

[0100] Furthermore, when the system is stable (S _Trans (approaching 0), weight allocation biased towards medium-speed layers (w) medium ≈0.6), the control is mainly based on accurately tracking the power target; when the system experiences severe transients (S _Trans (approaching 1), the weight allocation is biased towards the fast layer (w)_ fast ≈0.7), the control prioritizes counteracting internal disturbances and maintaining system stability; this dynamic fusion mechanism enables the control strategy of this invention to achieve an intelligent dynamic balance between accuracy and stability; the comprehensive control quantity u _Total (t) After safety verification, it is converted into specific control signals such as the PWM duty cycle of the electrolytic cell rectifier to drive the power change of the electrolytic cell (see Example 9 for details).

[0101] This step utilizes transient severity indicators to dynamically allocate the fusion weights of the three-layer control signals, establishing an intelligent arbitration mechanism: when the system is stable, the weight of the medium-speed layer is higher to pursue optimal power point tracking accuracy and economy; when the system suffers severe shocks, the weight of the fast feedforward compensation layer is immediately increased to prioritize system stability and rapid suppression of disturbances. This dynamic reconfiguration of the control strategy allows the system to automatically select the most suitable control mode under any operating condition, perfectly balancing the speed of frequency modulation response, the accuracy of tracking, and the stability of operation. It avoids the shortcomings of a single fixed strategy under complex operating conditions, improving the overall performance of the electrolyzer in frequency modulation.

[0102] Example 6 provides a specific numerical calculation case in a typical power grid frequency regulation event, specifically:

[0103] Step S61, Scene Setting: Rated Power P rated A 1000kW electrolytic cell is currently operating at P current The system is operating stably at a power output of 700kW. At time t0, it receives an AGC command from the power grid, requiring its power output to increase to P within 2 seconds. AGC =900kW; Due to the rapid power jump, a violent electrochemical transient was generated inside the electrolyzer.

[0104] Step S62, Transient severity assessment (corresponding to Example 4);

[0105] Within tens of milliseconds after t0, the system detected the following transient information:

[0106] The rate of change of impedance |dZ / dt| = 0.4Ω / s (preset Z) max =0.5Ω / s);

[0107] The rate of temperature change |dT / dt| = 0.25°C / s (preset T) max =0.5°C / s);

[0108] Second derivative of power |d 2 P / dt 2 |=40kW / s 2 (Preset P) max =50 kW / s2 );

[0109] Calculate transient indices for each dimension: S Z =0.4 / 0.5=0.8; S T =0.25 / 0.5=0.5; S P =40 / 50=0.8;

[0110] The comprehensive transient severity index S is calculated using a weighted formula. trans :S trans =0.4⋅S Z +0.3⋅S T +0.3⋅S P =0.4⋅0.8+0.3⋅0.5+0.3⋅0.8=0.32+0.15+0.24=0.71;

[0111] Step S63, PSO optimization strategy adjustment (corresponding to Example 4);

[0112] Due to S trans =0.71≥0.7, the PSO optimization algorithm determines that the system has entered the jump search mode; the algorithm will perform intelligent relocation for some particles; assuming that the relocation strategy based on real-time impedance mutation is selected at this time;

[0113] The system reads the impedance change of this event from the transient event characteristic matrix: ΔZ = 0.2 Ω; and calculates the adaptive gain: Kz = 0.5⋅(1+S trans =0.5⋅(1+0.71)=0.855; Assuming the optimal improvement direction given by the gradient ∇J of the current objective function is to increase power; the calculated transient relocation vector Δx _transient It will have a larger step size, pointing in the direction of increasing power, which is related to K. _Z It is proportional to the product of |ΔZ|;

[0114] Instead of blindly searching near the old optimal solution, the PSO algorithm decisively throws a portion of the search particles into a new, higher-power optimal region based on the strength and gradient direction of the physical perturbation, thus accelerating the search for the optimal operating point at 900kW.

[0115] Step S64, Multi-timescale control fusion (corresponding to Example 5);

[0116] The control layer utilizes S trans =0.71 to dynamically adjust the fusion weight:

[0117] w fast =0.3+0.4⋅0.71=0.3+0.284=0.584;

[0118] wmedium =0.6−0.3⋅0.71=0.6−0.213=0.387;

[0119] w slow =0.1−0.1⋅0.71=0.1−0.071=0.029;

[0120] Assume the outputs of the three-layer controller are as follows: Fast layer u fast Due to the power step causing a large number of bubbles to form, dθ / dt is positive. The fast layer outputs -30kW of compensation power to offset the overpotential caused by bubble adhesion and prevent voltage surges; the medium-speed layer... medium In order to track P opt (Approaching 900kW) target, with an output power regulation of +150kW; slow-speed layer u slow Temperature change is small, output is 0kW; calculate the overall control quantity: u total =0.584⋅(−30kW)+0.387⋅(+150kW)+0.029⋅(0kW)=−17.52kW+58.05kW=+40.53kW;

[0121] Although the goal for the mid-speed strata is to increase power (+150kW), the system is in a severe transient state (S... _Trans High), the weights of the fast layer (w) fast =0.584) was increased. Its -30kW compensation output effectively suppressed the aggressive adjustment of the mid-speed layer, making the final integrated control command (+40.53kW) smoother. This mechanism avoids overly drastic power regulation when the internal state of the system is unstable, effectively preventing the occurrence of malfunctions such as voltage overshoot and bubble layer rupture, and improving the stability and safety of the frequency modulation process while ensuring a fast response speed.

[0122] This numerical case clearly demonstrates how the present invention connects the high-level optimization algorithm search strategy with the low-level multi-scale control fusion mechanism through the core transient severity index, thereby achieving intelligent, coordinated, and robust control under complex working conditions.

[0123] Example 7 provides upstream technical support for concepts such as transient events and transient severity indicators, and explains how to accurately detect transient abrupt events from impedance time series before performing state estimation.

[0124] Step S71: Impedance time series differencing and smoothing;

[0125] To monitor the dynamic behavior of the system, we first used broadband impedance spectrum data Z... _spectrumSelect one or more characteristic frequency points (e.g., 1 kHz) in (t, ω) that are sensitive to changes in the system state, and extract the time series of the impedance magnitude |Z(t)| at these frequency points. To capture its dynamic changes, calculate the first-order difference (rate of change) and the second-order difference (acceleration) of the time series: dZ / dt=(|Z(t)|−|Z(t−Δt)|) / Δt; 2 Z / dt 2 =(dZ / dt(t)−dZ / dt(t−Δt)) / Δt; where Δt is the sampling time interval, for example, 10 milliseconds. Since differential operations amplify noise, a low-pass filter, such as a moving average filter with a window length of 50 milliseconds, needs to be applied to the calculated differential sequence to obtain a smoothed impedance change rate and acceleration sequence, providing a stable and reliable input for subsequent abrupt change detection.

[0126] Step S72: Apply the CUSUM adaptive threshold algorithm to perform mutation detection;

[0127] The CUSUM control chart is a sequential analysis technique for quickly detecting small but persistent process deviations. Compared to simple thresholding methods, it is more sensitive and has a lower false alarm rate.

[0128] In the specific implementation, two cumulative sum statistics S are initialized. + (0)=0 and S − (0) = 0; The CUSUM statistic is updated by recursively calculating the smoothed impedance change rate sequence dZ / dt: S + (t)=max(0, S+(t−1)+(dZ / dt−μ)−k⋅σ);

[0129] S − (t) = max(0, S−(t−1)−(dZ / dt−μ)−k⋅σ); where μ and σ are the mean and standard deviation of the dZ / dt sequence within the reference time window (e.g., the first 500 ms) under stable operating conditions, respectively, which constitute the background noise model of normal system fluctuations; k is a sensitivity parameter, usually taken as 0.5-1; the statistic S ⁺ (t) and S ⁻ (t) accumulates evidence of data deviations towards positive and negative directions, respectively. When either statistic exceeds the preset decision threshold h... _alarm For example, h _alarm When σ can be set to 5, a mutation event is determined to have occurred, and the current time is recorded as the mutation start time t. start .

[0130] When S ⁺ (t)>h _alarm or S ⁻(t)>h _alarm time (h) _alarm =5h), marking the mutation initiation time t _start Continue monitoring until S returns to normal, and record the end time t. _end Output mutation event label sequence t _start , t _end and CUSUM statistic series S⁺(t), S⁻(t);

[0131] In this embodiment, the use of an adaptive threshold enables the detection algorithm to be robust to baseline noise under different operating conditions, thus ensuring the accuracy of the detection.

[0132] Step S73: Feature quantification and feature matrix construction of mutation events;

[0133] Once a mutation event is detected, the system will capture the time window of that event (e.g., from t). start -50ms to the end of the mutation t end The smoothed impedance sequence Z within +50ms) _smooth (t). Analyze this data segment to quantify the characteristics of the mutation: mutation magnitude ΔZ: ΔZ = max(Z) - min(Z); rise time t rise :t _rise =argmax(Z)-t _start , where t _start This is the moment when the mutation begins and argmax(Z) reaches its peak; the duration τ trans τ _Trans =t _end -t _start， Among them, t _ end It is the end time;

[0134] Optionally, mutation classification can be performed based on features: if t _rise If the time is less than 20ms and ΔZ > 0.1Ω, it is marked as a bubble burst type;

[0135] If t _rise >100ms and τ _Trans If the time interval is >200ms, it is marked as concentration accumulation type; if it exhibits periodicity (FFT main frequency >1Hz), it is marked as oscillatory type; output abrupt change feature vector [ΔZ, t _rise , τ _Trans (type).

[0136] Step S74: Constructing the transient event feature matrix;

[0137] Summarize all mutation events within the last second, read the mutation feature vector of each event; construct the feature matrix E. _TransEach line corresponds to one event: [time, ΔZ, t] _rise , τ _Trans Type, Z _before Z _after ]; Calculate statistical characteristics: mutation frequency f _Trans =count / 1s, average amplitude ΔZ _avg σ ΔZ ; Mark the current dominant mutation type (the most frequently occurring type); Output the transient event feature matrix E _Trans Mutation statistical feature set f _Trans ΔZ _avg , σ ΔZ and dominant mutation type _ Type.

[0138] This matrix provides crucial, structured data input for calculating transient severity indices and for the PSO algorithm to select relocation strategies in jump search mode in subsequent embodiments.

[0139] Example 8 illustrates the objective function evaluation and iterative update method in the particle swarm optimization process, further including:

[0140] Step S81: Objective function evaluation and constraint handling;

[0141] In each iteration of PSO, the candidate solutions represented by each particle need to be evaluated to determine their quality, i.e., their fitness value needs to be calculated. This includes calculating multi-dimensional performance metrics and handling various operational constraints.

[0142] Specifically, the implementation process of step S81 includes:

[0143] Step S811: Calculate the multidimensional frequency modulation performance index;

[0144] For each particle i, the corresponding power scheme P i To evaluate its performance across multiple dimensions, a multidimensional performance index matrix is ​​formed [J1] _i J2 _i J3 _i J4 _i Frequency modulation accuracy index J1 i J1 i =∣P actual_i -P AGC | / P rated , where P _ actual_i It is based on model prediction, when the execution power P _i The actual response power at that time; the smaller this index, the better; response speed index J2 i J2i =max(0, 1−t) 90_i / t required ), where t _ 90_i This is the predicted time required to reach 90% of the target power, t_ "Required" refers to the response time required by the power grid; a higher value is better; the operational efficiency index J3 i J3 i =ηi / η rated , where η i The efficiency is obtained by looking up a table or simplifying the efficiency model, in power P i The higher the operating efficiency, the better.

[0145] Frequency modulation stability index J4 i J4 i =1 / (1+σ P_i / Prated ), where σ _ P_i This is the standard deviation of the predicted power fluctuation; the larger this indicator is, the better.

[0146] Step S812: Calculate dynamic constraint boundaries and perform softening treatment;

[0147] Read current operating status data and dynamically adjust constraint boundaries: power upper limit P _max_dyn =P _rated ×(1-0.2×θ) Considering the effect of bubbles, the lower power limit P _min_dyn =max(0.1×P _rated P _min _design);

[0148] Calculate the temperature-dependent power change rate limit: ramp _limit_T =ramp _base ×f _Thermal ;

[0149] Soften the hard constraints: Define the feasible region boundary layer thickness = 0.05 × (P) _max -P _min This means allowing for 5% boundary flexibility, balancing feasibility and search efficiency;

[0150] Applying a quadratic penalty function within the boundary layer: penalty _soft =(d / thickness) 2 , where d is the distance to the boundary;

[0151] Output dynamic constraint set C _dyn = P _min_dyn P _max_dyn , ramp _limit_T and softening parameter set.

[0152] Step S813: Constraint violation assessment and penalty function calculation;

[0153] Specifically, for each particle position x _i Read the corresponding power P _i rate of change dP / dt _i ;

[0154] Check power constraints: v _P ower _i =max(0, P) _min_dyn -P _i )+max(0,P _i -P _max_dyn );

[0155] Check the rate of change constraint: v _ramp_i =max(0,|dP / dt) _i |-ramp _limit_T );

[0156] Check temperature constraints: Predicted temperature T _Pred_i v _Temp_i =max(0,T) _min -T _Pred_i )+max(0,T _Pred_i -T _max );

[0157] Applying an adaptive penalty function: penalty _i =λ _P ×v _P ower _i 2 +λ_r×v _ramp_i 2 +λ _T ×v _Temp_i 2 The penalty factor λ increases with the number of iterations;

[0158] Output constraint violation degree vector V _c onstraint[N] and the penalty function value vector Penalty[N].

[0159] Step S814: Perform multi-objective aggregation and fitness calculation.

[0160] Read the multidimensional performance index matrix J and the penalty function value vector Penalty, and apply dynamic weight aggregation: w=[w1, w2, w3, w4]=[0.3, 0.3, 0.3, 0.1], which can be adjusted according to grid demand; calculate the weighted objective function: J _weighted_i =Σ(w _j ×J_j_i ); Integrating constraint penalties: F _raw_i =J _weighted_i +Penalty _i Perform a fitness transformation, converting it into a maximization problem: F _fitness_i =1 / (1+F _raw_i Apply fitness scaling to enhance selection pressure: F _scaled_i =F _fitness_i ^α Where α=1.5; output fitness value matrix F _fitness [N] and optimal fitness index i _best .

[0161] This step achieves multi-physics coupled modeling, integrating physical characteristics such as thermal deviation factor and bubble dynamics into the constraint boundary; dynamic weight adaptation adjusts the importance of each performance index in real time according to grid demand; and soft and hard constraint co-optimization balances search efficiency and feasibility through softening processing and adaptive penalty functions. Compared with traditional PSO, this invention improves the frequency regulation response speed of the electrolyzer, reduces energy consumption, and enhances the economy and reliability of AGC control.

[0162] Step S82: Particle iterative update and convergence judgment;

[0163] After calculating the fitness values ​​of all particles, the algorithm will update the individual optimal position p of each particle. _best_i and the global optimal position g of the entire population _best .

[0164] Specifically, based on the adaptive PSO parameter set⊙ _PSO (t), update particle velocity v _i (t+1)=ω×v _i (t)+c1×r1×(p _best_i -x _i )+c2×r2×(g _best -x _i ) and position x _i (t+1)=x _i (t)+v _i (t+1); where ω(t), c1(t), and c2(t) are the dynamic inertia weights and learning factors that adapt to the working conditions and are calculated in Example 3. r1 and r2 are two random numbers uniformly distributed in the interval [0, 1], introducing randomness into the search process; intelligent reflection is performed on particles that touch the boundary: x _reflected =x _boundary -0.5×(x _i -x _boundary )×exp(-S _Trans Output the optimal power allocation scheme P_opt (t) and predictive control parameter K _Pred =[K _P K _i ,K_d].

[0165] This iterative process is not indefinite, but terminates after a preset convergence condition is met; the possible termination conditions include at least one of the following:

[0166] Reaching the maximum number of iterations: For example, the number of iterations reaches 50; Fitness improvement stagnates: The fitness value of the global optimum improves by less than a very small threshold (e.g., 0.001) in several consecutive iterations (e.g., 3).

[0167] Computation timeout: To ensure real-time control, the total time consumed by the entire optimization calculation process exceeds the upper limit (e.g., 0.3 seconds); once any termination condition is met, the iteration stops, and the currently locked global optimal solution g is set. _best As the optimal power allocation scheme P _opt (t) is output.

[0168] The particle swarm converges efficiently to the global optimum under the guidance of dynamic parameters, while taking into account both real-time performance and robustness, providing a readily implementable optimal power allocation scheme for power systems.

[0169] Example 9: Description of obtaining the comprehensive control quantity u through dynamic weighted fusion. _Total Following (t), it further includes:

[0170] Step S91: Power control signal generation and security verification;

[0171] This step aims to transform the comprehensive control quantity u calculated at the algorithm level. total (t) is transformed into physical instructions that can be directly executed by power electronic hardware, and key security mechanisms are embedded in this process.

[0172] Specifically, the integrated control quantity is converted into a control signal (PWM, Pulse Width Modulation) for the electrolytic cell rectifier, and the PWM duty cycle is calculated: D = (P _base +u _Total ) / (U _dc ×I _max ); where P _base It is the current base power or the actual power at the previous moment; U _ dc It is the voltage of the rectifier's DC bus; I _max This is the maximum current that the rectifier is allowed to pass through. This formula linearly maps the power command to the duty cycle.

[0173] Before outputting the duty cycle D, a rigorous safety boundary check must be performed. The electrolyzer and its rectifier hardware have inherent, inviolable safe operating ranges. For example, to avoid control saturation and extremely low efficiency regions, the duty cycle D must be limited to a reasonable range, exemplarily [0.1, 0.95]; the system checks whether the calculated D violates this constraint.

[0174] Calculate the expected response time t _Pred and efficiency η _Pred If the safety constraints are violated (D<0.1 or D>0.95), the limiting protection mechanism will be activated; the electrolytic cell power control signal S will be output. _control (t); whereby the limiting protection mechanism will force the calculated duty cycle D to be clamped at the nearest boundary. For example, if D=0.98 is calculated, the final output duty cycle will be set to 0.95; if D=0.05 is calculated, the final output will be set to 0.1; this step is a key line of defense to ensure hardware safety and prevent equipment damage or hard protection shutdown caused by excessively large or small algorithm output instructions.

[0175] Step S91: Real-time evaluation and feedback of frequency modulation performance;

[0176] This step constructs a closed-loop learning and evaluation system, enabling the control method not only to perform frequency modulation tasks, but also to quantitatively evaluate its effectiveness after the task is completed, and to use the evaluation results for future decision optimization.

[0177] Specifically, after a frequency modulation event ends, the system reads the actual power response sequence P for that period. _actual (t) and the power demand sequence P of the grid AGC _AGC (t), and calculate a series of performance indicators: frequency modulation accuracy indicator D1: D1=1−∣P actual -P AGC | / P AGC This is used to measure how close the steady-state power at the end of the response is to the target power; the response speed index D2: D2 = min(1, t standard / t actual ), where t _actual It is the actual response time to achieve, for example, 90% of the target power, t_ "standard" refers to the standard response time required by the power grid specifications; tracking accuracy index D3: D3 = 1 − RMSE(P actual P AGC ) / P ratedRMSE is the root mean square error between the actual response trajectory and the command trajectory, which measures the tracking smoothness and accuracy of the entire dynamic process. Based on this, the comprehensive performance index D = w1⋅D1 + w2⋅D2 + w3⋅D3 can be calculated, for example, using weights w1 = 0.25, w2 = 0.25, and w3 = 0.5. The frequency modulation mileage M is calculated, defined as the integral of the absolute value of the actual power change rate over time: M = constant | dP actual / dt∣dt, which quantifies the total contribution of the electrolyzer in this service.

[0178] Most importantly, this step enables feedback to the historical knowledge base. The system combines the comprehensive performance index D of this event with the operating condition characteristics that triggered this event (such as the transient event feature matrix E in Example 7). _Trans and the transient severity index S in Example 4 _Trans The successful power adjustment strategies, along with the actual power adjustment strategies that ultimately prove successful, are packaged into a new success story. This story will be added to or used to update the historical knowledge base mentioned in Example 4. This feedback loop enables the entire control system to learn and self-optimize online, which is the core of achieving intelligent control.

[0179] Example 10 describes the method for determining key parameters and the maintenance mechanism of the historical knowledge base, including:

[0180] Explanation of the method for determining key parameters: The values ​​of several key parameters involved in this invention are not arbitrarily set, but are based on physical models, experimental data and control theory, and have clear engineering and scientific basis.

[0181] Regarding the shape factor β (e.g., 0.85): This parameter appears in the bubble coverage calculation formula of Example 2, and it corrects for model errors when the shielding effect of three-dimensional spherical or ellipsoidal bubbles on the electrodes is equated to two-dimensional planar coverage. Ideal two-dimensional planar coverage corresponds to β=1. In practical applications, the value of β needs to be calibrated through specialized experiments. A feasible method is to build an electrolytic cell test bench with a transparent observation window and simultaneously observe at different current densities using high-frequency impedance spectroscopy measurement technology and a high-speed camera. The true bubble coverage can be obtained through image processing algorithms. By fitting the true coverage with the coverage calculated by measuring the double-layer capacitance, the β value that best matches the actual situation of the electrolytic cell can be calibrated.

[0182] Regarding the gain / compensation coefficient (e.g., K) _b1 K _c1(e.g., ) These parameters appear in the fast control layer compensation model of Example 5. They are typical control system gain coefficients, and their values ​​are usually determined following the standard process of system identification and controller tuning. Specifically, firstly, system identification of the electrolytic cell is required, that is, applying a series of specific excitation signals (such as power steps, pseudo-random binary sequences, etc.) on the test platform and recording its dynamic response (such as voltage and impedance changes). Based on these input and output data, a mathematical model describing the dynamic behavior of the electrolytic cell can be established. Then, based on the model, using classical control theory (such as root locus method, frequency response method) or modern control theory, the controller parameters that enable the system to obtain the desired closed-loop performance (such as fast response, small overshoot, and strong robustness) are calculated. These parameters are K. _b1 K _c1 wait.

[0183] Initialization and maintenance of the historical knowledge base; the historical knowledge base is a key component in the PSO intelligent relocation strategy of Example 4, and its construction and maintenance mechanism is as follows:

[0184] The knowledge base's data structure: The knowledge base is a database that stores working condition-policy pairings; each record contains two parts: a working condition feature vector F. hist It consists of a series of quantitative data that can uniquely describe a specific working condition, such as transient severity indicators and transient event characteristics; and the corresponding optimal adjustment strategy ΔP, which has been verified as successful. opt For example, successful power adjustment or PSO jump vector.

[0185] Initialization (Cold Start): In the initial deployment of a new electrolyzer system, the knowledge base can be empty. At this time, the PSO's jump search will rely entirely on Strategy 1 based on real-time impedance mutations. Optionally, to accelerate the learning process, it can be pre-filled. The pre-filled data can come from extensive offline simulations of detailed electrochemical models of the electrolyzer, simulating various frequency modulation and fault conditions, and storing the optimal solutions found in the simulations into the knowledge base; alternatively, it can be from analyzing historical operating data from similar electrolyzers to extract typical successful response cases.

[0186] Online Updates and Maintenance: As described in Example 9, each time a frequency modulation service is successfully completed (judged by the high comprehensive performance index D), the system packages the service condition feature vector and success strategy into a new record. Before storing it in the knowledge base, a similarity comparison is performed with existing records in the base. Only when the service condition is considered novel (i.e., its similarity to all historical service conditions is below a certain threshold) will the record be added. This mechanism ensures that the knowledge base stores diverse and valuable experiences. Optionally, a forgetting mechanism can be introduced, such as assigning weights or freshness to each record. For old records that have not been successfully matched for a long time, or that have been updated or covered by better strategies, their weights can be gradually reduced or they can be removed to keep the knowledge base concise and efficient.

[0187] Example 11 describes a preferred implementation method for the process of generating feedforward compensation control quantities in a fast control layer, which more accurately reflects the intrinsic coupling relationship of the electrochemical process and improves the robustness of the control signal. Specifically, the optimization is as follows:

[0188] Step S111: Coupling correction of the compensation signal;

[0189] In Example 5, the bubble compensation power ΔP _bubble With concentration gradient compensation power ΔP _ conc These effects are calculated independently and then directly superimposed. However, in a real electrochemical environment, these two effects are not absolutely independent but rather mutually influential. For example, excessive bubble coverage on the electrode surface can hinder ion transport paths and exacerbate local concentration polarization. To reflect this physical coupling in the model, a coupling correction mechanism is introduced.

[0190] Specifically, after calculating the original bubble compensation power ΔP... _bubble_raw And the original concentration gradient compensation power ΔP _conc_raw Afterwards, instead of adding them directly, the coupling coefficient α is calculated first. _couple α couple =1−exp(−k couple ⋅θ⋅∣∇c∣); where θ and |∇c| are the real-time bubble coverage and ion concentration gradient, respectively; k _couple These are preset coupling strength parameters; it can be seen that when either θ or |∇c| is very small, α _couple A value close to 0 indicates that the coupling effect is not significant; when both are significant, α _couple A value close to 1 indicates a strong coupling effect; this coupling coefficient is used to correct the original compensation amount: ΔP bubble_coupled =ΔP bubble_raw ⋅(1+w b ⋅α couple );ΔP conc_coupled =ΔP conc_raw⋅(1+w c ⋅α couple ); where w b and w c It is to adjust the weights, for example, w can be taken. b =0.3, w c =0.2. This means that when the coupling effect is significant, both compensation powers will be moderately amplified to more accurately offset this superimposed, nonlinear disturbance; check the consistency of the compensation direction, if sign(ΔP) _bubble )≠ sign(ΔP _conc ), applying conflict resolution strategy: taking 80% of the larger absolute value; compensation power after coupling correction. ΔP _bubble_coupled ΔP _conc_coupled Read the compensation power after coupling correction and calculate the total compensation ΔP. _comp_raw =ΔP _bubble_coupled +ΔP _conc_coupled Processed by a second-order Butterworth low-pass filter: H(s) = 1 / (1 + sqrt(2) × s / ω) _c +s 2 / ω _c2 Application output limiting: ΔP _comp_limited =sign(ΔP _comp_filtered )×min(|ΔP _comp_filtered |,0.3×P _rated Add dead zone handling, |ΔP _comp Output 0 when |<1kW; Output transient compensation power ΔP _comp (t) and the compensation activation flag comp_active.

[0191] Step S111: Variable bandwidth adaptive filtering of the compensation signal;

[0192] To avoid the feedforward compensation signal being too sharp and causing system oscillation, it needs to be low-pass filtered. However, using a fixed-bandwidth filter presents two major challenges: if the bandwidth is too wide, the noise suppression capability is poor, and it may lead to overcompensation during severe transients; if the bandwidth is too narrow, the response speed of the compensation signal will decrease, making it unable to keep up with rapid electrochemical disturbances.

[0193] Therefore, this embodiment adopts a variable bandwidth self-adaptive filtering strategy. The core idea of ​​this strategy is to make the filter cutoff frequency f c It is dynamically related to the system's stability. Specifically, the filter's cutoff frequency f c Compared with the transient severity index S calculated in Example 4 _Trans Negative correlation: f c (t)=f c，base⋅(1−k filter ⋅S trans (t)); where f _c，base This is the fundamental cutoff frequency, such as 100Hz, which represents the highest expected compensation speed of the system in a steady state; k _ filter This involves adjusting the sensitivity, for example, to 0.5; and reducing the bandwidth to avoid overcompensation during severe transients.

[0194] The filter used can be a second-order Butterworth low-pass filter, whose transfer function is H(s) = 1 / (1 + sqrt(2) (s / ω) c )+(s / ω c ) 2 ), where ω c (t)=2πf c (t); Total compensation power ΔP after coupling correction comp_raw =ΔP bubble_coupled +ΔP conc_coupled The transient compensation power ΔP will be processed through this variable bandwidth filter to obtain the final smooth and robust transient compensation power ΔP. comp (t), which is the output u of the fast layer. fast (t); Subsequent processing includes: applying output limiting: ΔP _comp_limited =sign(ΔP _comp_filtered )×min(|ΔP _comp_filtered |,0.3×P _rated Add dead zone handling, |ΔP _comp | Output 0 when < 1kW; Output transient compensation power ΔP _comp (t) and the compensation activation flag comp_active.

[0195] This method realizes the intelligence of the feedforward compensation channel: when the system is running smoothly (S _Trans Small, f c High-frequency noise (HF) and rapid compensation channel response effectively filter out high-frequency noise and quickly compensate for minor disturbances; when the system experiences severe transients (S... _Trans Big, f c (Low bandwidth) By reducing bandwidth, the compensation signal is smoothed, avoiding the exacerbation of system oscillations due to excessive feedforward, thus improving the stability of the control system.

[0196] 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 control method for an electrolytic cell participating in frequency regulation ancillary services of a power system, characterized in that, include: Collect AGC frequency regulation commands and power demand from the power grid, electrolyzer operating status data, and broadband impedance spectrum data to generate a synchronization dataset; Based on the synchronous dataset, millisecond-level electrochemical transient features are extracted and interface state estimation is performed to obtain interface state estimates and transient feature vectors. By combining the interface state estimate, transient feature vector, and AGC power requirement, transient perception multi-objective optimization is performed to determine the optimal power allocation scheme. Based on the optimal power allocation scheme and the interface state estimate, multi-time-scale coordinated control is performed, and a power control signal for the electrolyzer is output. The transient sensing multi-objective optimization includes: constructing a guided initial particle swarm using transient feature vectors; dynamically adjusting the control parameters of the particle swarm optimization based on the electrochemical environment reflected by the interface state estimates during the iteration process; evaluating the comprehensive fitness of each particle in the particle swarm in terms of frequency modulation accuracy, response speed, and operating efficiency in conjunction with the AGC power requirements, and obtaining the corresponding fitness value; and guiding the particle swarm to explore and find a better solution space based on the fitness value and the dynamically adjusted control parameters until the convergence condition is met, forming the global optimal solution, which is the optimal power allocation scheme. The process of extracting millisecond-level electrochemical transient features and performing interface state estimation includes: analyzing broadband impedance spectral data in the synchronous dataset, decoupling and obtaining preliminary quantization values ​​characterizing bubble coverage and ion concentration gradient; constructing and applying an extended Kalman filter to receive the preliminary quantization values ​​as observations; obtaining an accurate depiction of the electrochemical interface state of the electrolyzer through recursive estimation using the extended Kalman filter, and outputting interface state estimates; and constructing and outputting transient feature vectors based on the interface state estimates. The multi-timescale coordinated control includes: a fast control layer, which generates feedforward compensation control quantities to offset transient electrochemical disturbances based on interface state estimates; a medium-speed control layer, which generates feedback regulation control quantities to achieve target power tracking based on the deviation between the optimal power allocation scheme and the actual power of the electrolyzer; and a slow control layer, which generates preventive management control quantities to maintain the long-term thermal balance of the system based on the evolution trend of historical temperature data. The three control quantities are then synergistically integrated to form a comprehensive control quantity.

2. The method according to claim 1, characterized in that, Guiding the particle swarm to explore and find better solutions, including: Multi-dimensional transient information, including impedance change rate, temperature change rate, and power second derivative, is extracted from the electrolytic cell operating status data and transient feature vectors. By performing weighted aggregation and quantitative evaluation on multi-dimensional transient information, a transient severity index characterizing the current stability of the system is generated; Based on the value of the transient severity index, the particle swarm optimization strategy is switched between preset progressive search, extended search and jump search modes.

3. The method according to claim 2, characterized in that, When the optimization strategy switches to jump search mode, the particle swarm is guided to explore and find better solutions, and intelligent relocation is performed on specific particles in the swarm, including: Choose one of the following strategies to generate the jump vector for relocation: Using the gradient of the objective function of the current optimization problem as the direction, the transition amplitude is calibrated by using the real-time impedance change recorded in the transient event feature matrix, and a transient relocation vector is constructed. The current working condition features, which are composed of interface state estimates, are matched with a pre-configured knowledge base. If similar historical working conditions are found, the corresponding adjustment strategies are used to construct an experience-guided jump vector. Update the particle's position based on the jump vector.

4. The method according to claim 1, characterized in that, Dynamically adjusting the control parameters for particle swarm optimization includes: From the electrolyzer operating status data and interface state estimates, the effects of temperature on operating efficiency, the impact of bubble dynamics on control stability, and the limitations of ion concentration gradient on mass transfer capacity are quantified respectively. Based on this, factors characterizing the multi-physics state are constructed, including at least the thermal deviation factor, the bubble dynamics factor, and the concentration gradient factor. By establishing and utilizing a preset mapping relationship, multiple physics field factors are coupled into a real-time adjustment basis for particle swarm optimization control parameters, generating an adaptive PSO control parameter set.

5. The method according to claim 4, characterized in that, Generate an adaptive PSO control parameter set, including: The thermal deviation factor, bubble dynamic factor, and concentration gradient factor are multiplied and fused into a comprehensive adjustment coefficient. The preset basic inertia weight is modulated using a comprehensive adjustment coefficient to generate a dynamic inertia weight.

6. The method according to claim 1, characterized in that, Constructing a guided initial particle swarm includes: The system deviation information contained in the transient feature vector is analyzed and transformed into a transient offset prediction quantity that characterizes the drift direction and amplitude of the optimal region in the solution space. Based on the previously found historical best solution, and guided by the transient offset prediction, the core particle swarm is generated in the target area with potential greater than the threshold. Randomly disperse the exploration particle swarm within the algorithm's feasible region; The core particle swarm and the exploration particle swarm are merged to form the initial particle swarm.

7. The method according to claim 1, characterized in that, Decoupling and obtaining preliminary quantization values ​​includes: Broadband impedance spectrum data is divided into high-frequency and low-frequency bands in the frequency domain; For high-frequency impedance data that reflects the rapid dynamics of the electrode interface, parameters related to the effective electrode area are extracted and a preliminary quantitative value of the bubble coverage is calculated. Based on the low-frequency impedance data characterizing the slow ion diffusion process in the electrolyte, parameters related to mass transfer limitation are calculated, and a preliminary quantitative value of the ion concentration gradient is derived.

8. The method according to claim 7, characterized in that, The estimation of bubble coverage includes: The high-frequency impedance data is fitted to the preset equivalent circuit model of the electrode interface to identify and extract the real-time double-layer capacitance parameters. A conversion relationship based on a physical model is established, and the real-time double-layer capacitance parameters are compared with the pre-calibrated reference double-layer capacitance under bubble-free conditions to calculate the preliminary quantitative value of bubble coverage.

Images