Method for participation of modular large-scale hydrogen production power in primary frequency regulation control of power system

Abstract

The application provides a method for participating in primary frequency modulation control of a power system based on a modular large-scale hydrogen production power supply, comprising: obtaining original data streams from a high-precision synchronous phasor measurement unit, an electrolytic cell and a sub-module, performing multi-scale decomposition and noise reduction processing after time synchronization and abnormality detection repair; using unscented Kalman filtering for state estimation, combining empirical mode decomposition and Hilbert transform to extract frequency disturbance characteristics and estimate grid inertia; using an established electrolytic cell dynamic model to identify and optimize parameters through recursive least squares; designing a control strategy based on a multi-time scale control framework, optimizing virtual synchronous generator parameters through online Gaussian process regression, combining scenario tree stochastic model predictive control and minimax regret learning to realize frequency modulation control; finally, coordinating control through a distributed consistency algorithm and performing multi-objective performance evaluation optimization. The system's frequency modulation performance and operating efficiency are improved.

Description

Technical Field

[0001] This invention relates to power electronics technology, and in particular to a method for primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source. Background Technology

[0002] With the large-scale grid integration of renewable energy, the frequency fluctuation problem in power systems has become increasingly prominent. Traditional primary frequency regulation mainly relies on the rotational inertia of conventional units and primary frequency regulation characteristics. However, with the integration of a high proportion of renewable energy, the system inertia is significantly reduced, making traditional frequency regulation methods insufficient to meet system demands. Large-scale hydrogen production, as an important adjustable load, features fast response and a wide adjustment range, making its introduction into power system primary frequency regulation significant. In particular, modular hydrogen production systems, through flexible power regulation capabilities, can effectively participate in grid frequency regulation while also achieving efficient energy storage and utilization, providing crucial technical support for building new power systems.

[0003] Currently, research on hydrogen production systems' participation in power frequency regulation mainly focuses on several aspects: first, the study of hydrogen production load characteristics based on simplified models, describing the dynamic characteristics of the electrolyzer through first- or second-order inertial elements; second, the adoption of single-time-scale control strategies, primarily achieving power point tracking through PI control or model predictive control; and third, centralized control methods for small-scale hydrogen production systems, achieving system scheduling through a unified control center. In data processing, traditional filtering and state estimation methods, such as Kalman filtering and extended Kalman filtering, are mainly used; in model identification, linear system identification methods, such as least squares and subspace identification, are primarily employed; and in control strategies, optimization methods based on single performance indicators are mainly used.

[0004] However, existing technologies face several challenges when applied to large-scale modular hydrogen production systems: First, in terms of data processing, traditional time synchronization and anomaly detection methods struggle to handle high sampling rates and heterogeneous data streams from multiple sources, particularly when dealing with sudden anomalies, resulting in poor data repair quality. Second, existing electrolyzer models fail to accurately describe the system's dynamic characteristics under different operating conditions, especially the nonlinear response characteristics during rapid power changes. Third, single-time-scale control strategies cannot simultaneously balance rapid system response and economic efficiency, particularly when coordinating the control of multiple subsystems, making it difficult to balance real-time performance and optimization. Furthermore, traditional centralized control methods suffer from heavy communication burdens and high computational complexity when dealing with numerous distributed subsystems, limiting system scalability. Finally, existing performance evaluation methods often focus on a single indicator, failing to comprehensively reflect the system's dynamic and economic performance across multiple time scales. These technical challenges severely restrict the effectiveness of large-scale hydrogen production systems in grid primary frequency regulation. Summary of the Invention

[0005] The purpose of this invention is to provide a method for primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source, in order to solve one of the aforementioned problems in the prior art.

[0006] The technical solution provides a method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source, comprising the following steps:

[0007] Step S1: Obtain the original PMU data stream P(t) from the high-precision synchronous phasor measurement unit, the original electrolytic cell data stream E(t) from the electrolytic cell, and the original submodule current data stream I(t) from the submodule; input P(t), E(t), and I(t) into the precise time protocol module for time synchronization to obtain the time-aligned data stream D(t), and input it into the robust principal component analysis module for anomaly detection and repair to obtain the repaired data stream D'(t); apply the biorthogonal wavelet transform to the repaired data stream D'(t) for multi-scale decomposition and noise reduction, and output the noise-reduced wavelet coefficients W' and eigenvector F;

[0008] Step S2: Input the denoised wavelet coefficients W' and eigenvector F into an unscented Kalman filter to obtain the system state vector x_est; extract the frequency component f_est from the state vector x_est and perform empirical mode decomposition and Hilbert transform to obtain the frequency disturbance feature matrix F_dist; based on the frequency disturbance feature matrix F_dist and the system state vector x_est, calculate the grid inertia estimate M_est using the improved area method.

[0009] Step S3: Input the electrolyzer data stream E(t) into the improved Hammerstein-Wiener model to obtain the electrolyzer dynamic model structure HW; apply the recursive least squares method to the model structure HW to identify parameters and obtain the model parameter vector θ_est; input the model parameter vector θ_est into the model validation module for parameter optimization and obtain the final model parameter vector θ_final.

[0010] Step S4: Input the system state vector x_est, frequency disturbance feature matrix F_dist, grid inertia estimate M_est, and model parameter vector θ_final into the multi-time-scale control framework to obtain the control framework structure MSC; apply online Gaussian process regression to the control framework MSC to obtain the virtual synchronous generator parameters J_vsg and D_vsg; construct and solve the stochastic model predictive control problem based on the scenario tree to obtain the optimal control sequence u*; apply the online gradient descent method to solve the minimax regret learning problem to obtain the long-term frequency regulation strategy L_opt.

[0011] Beneficial results: Through multi-objective optimization and performance evaluation, a balance between system operating efficiency and frequency regulation performance was achieved, providing reliable control support for large-scale hydrogen production systems to participate in grid frequency regulation. Attached Figure Description

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

[0013] Figure 2 This is a topology diagram of a hydrogen production power source.

[0014] Figure 3 It is the equivalent circuit of an electrolytic cell.

[0015] Figure 4 This is a block diagram of the hydrogen production power source participating in frequency modulation control. Detailed Implementation

[0016] like Figure 1 As shown, according to one aspect of this application, a method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source is provided, comprising the following steps:

[0017] Step S1: Obtain the original PMU data stream P(t) from the high-precision synchronous phasor measurement unit, the original electrolytic cell data stream E(t) from the electrolytic cell, and the original submodule current data stream I(t) from the submodule; input P(t), E(t), and I(t) into the precise time protocol module for time synchronization to obtain the time-aligned data stream D(t), and input it into the robust principal component analysis module for anomaly detection and repair to obtain the repaired data stream D'(t); apply the biorthogonal wavelet transform to the repaired data stream D'(t) for multi-scale decomposition and noise reduction, and output the noise-reduced wavelet coefficients W' and eigenvector F;

[0018] Step S2: Input the denoised wavelet coefficients W' and eigenvector F into an unscented Kalman filter to obtain the system state vector x_est; extract the frequency component f_est from the state vector x_est and perform empirical mode decomposition and Hilbert transform to obtain the frequency disturbance feature matrix F_dist; based on the frequency disturbance feature matrix F_dist and the system state vector x_est, calculate the grid inertia estimate M_est using the improved area method.

[0019] Step S3: Input the electrolyzer data stream E(t) into the improved Hammerstein-Wiener model to obtain the electrolyzer dynamic model structure HW; apply the recursive least squares method to the model structure HW to identify parameters and obtain the model parameter vector θ_est; input the model parameter vector θ_est into the model validation module for parameter optimization and obtain the final model parameter vector θ_final.

[0020] Step S4: Input the system state vector x_est, frequency disturbance feature matrix F_dist, grid inertia estimate M_est, and model parameter vector θ_final into the multi-time-scale control framework to obtain the control framework structure MSC; apply online Gaussian process regression to the control framework MSC to obtain the virtual synchronous generator parameters J_vsg and D_vsg; construct and solve the stochastic model predictive control problem based on the scenario tree to obtain the optimal control sequence u*; apply the online gradient descent method to solve the minimax regret learning problem to obtain the long-term frequency regulation strategy L_opt.

[0021] According to one aspect of this application, step S5 is also included:

[0022] Step S5: Input the control framework MSC, virtual synchronous generator parameters J_vsg and D_vsg, optimal control sequence u*, and long-term frequency regulation strategy L_opt into the distributed consensus algorithm module to obtain the global coordinated control instruction set C; input the coordinated control instruction set C into the real-time control execution module and collect feedback data to obtain the control execution result R(t); input the control execution result R(t) into the multi-objective performance evaluation module for optimization, and output the optimized control parameter set P_opt and performance evaluation report E_report.

[0023] In this embodiment, firstly, unified processing and high-quality feature extraction of multi-source heterogeneous data are achieved, providing a reliable data foundation for system control. Secondly, through accurate state estimation and model identification, the dynamic characteristics of the system and the grid state are accurately grasped. Thirdly, a multi-timescale control framework is adopted to achieve rapid response and economic optimization of the hydrogen production system to the grid frequency. Finally, through distributed coordinated control and multi-objective optimization, the efficient operation of the large-scale system is ensured. This systematic solution not only significantly improves the performance of the hydrogen production system in grid frequency regulation but also achieves economic efficiency and reliability in system operation, providing important technical support for the application of large-scale hydrogen production systems in power systems. Especially when dealing with large-scale distributed systems, it can effectively coordinate the operation of multiple subsystems, optimize overall performance, and ensure the scalability and adaptability of the system.

[0024] According to one aspect of this application, step S1 specifically comprises:

[0025] Step S11: Collect basic data including voltage, current, and phase angle from the synchronous phasor measurement unit, preprocess it through the adaptive sampling rate adjustment module to obtain the original PMU data stream P(t) with a sampling frequency of a predetermined frequency; collect parameters such as voltage, current, temperature, pressure, and electrolyte concentration from the electrolytic cell, perform preliminary noise reduction through a multivariable composite Kalman pre-filter to obtain the original electrolytic cell data stream E(t) with a sampling frequency of 10Hz; collect inductor current data from N sub-modules, compress the data through the distributed data compression sampling module to obtain the original sub-module current data stream I(t) with a sampling frequency of 10kHz; input P(t), E(t), and I(t) into the data quality evaluation module, calculate the data quality index by combining it with pre-stored historical data, and obtain the data stream set D_q(t) with quality score; N is a predetermined natural number;

[0026] Step S12: Input the data stream set D_q(t) into the precise time protocol module, and perform timestamp alignment using the pre-stored time synchronization parameter set T to obtain the time-aligned data stream D(t); input D(t) into the tensor decomposition anomaly detection module, and obtain the low-rank tensor L and sparse tensor S through multi-dimensional tensor decomposition, and identify the abnormal data point set A from S; input the abnormal data point set A into the deep autoencoder network for data repair to obtain the repaired dataset R; merge R with the normal data in D(t), and perform local manifold alignment processing to obtain the repaired data stream D'(t); apply adaptive multi-scale filtering to D'(t) for preliminary noise reduction to obtain the filtered data stream D_f(t);

[0027] Step S13: Input the filtered data stream D_f(t) into the variational mode decomposition module to obtain multiple intrinsic mode functions (IMFs)_k; apply bioorthogonal wavelet transform to each IMF_k for multi-scale decomposition to obtain approximation coefficients a_J and detail coefficients d_j; input all coefficients into the adaptive threshold calculation module to generate the optimal threshold λ_j based on the noise standard deviation of each layer, and perform soft thresholding on the coefficients to obtain the denoised wavelet coefficients W'; input W' into the multi-dimensional feature extraction module to calculate statistical features, entropy features, and time-frequency features to obtain the feature vector F; where statistical features include mean, variance, skewness, and kurtosis; entropy features include sample entropy, fuzzy entropy, and permutation entropy; and time-frequency features include Hilbert spectrum and wavelet packet energy.

[0028] By employing a high-precision synchronous phasor measurement unit to acquire PMU data streams, electrolytic cell data streams, and submodule current data streams, and performing time synchronization and data repair, unified processing of multi-source heterogeneous data is achieved. Specific innovations include: First, an adaptive sampling rate adjustment module is used to preprocess PMU data to ensure consistent sampling frequencies; second, a multivariate composite Kalman pre-filter is used to perform preliminary noise reduction on electrolytic cell data, improving data quality; third, a distributed data compression sampling module is used to compress submodule current data, reducing data transmission burden. Time synchronization is achieved through a precise time protocol module, achieving microsecond-level data alignment accuracy. Tensor decomposition anomaly detection and a deep autoencoder network are used for data repair, significantly improving data reliability. Finally, variational mode decomposition and bioorthogonal wavelet transform are used for multi-scale decomposition and noise reduction, extracting rich statistical, entropy, and time-frequency features. This multi-level data processing scheme not only ensures data time consistency and reliability but also achieves effective noise reduction and feature extraction, providing a high-quality data foundation for subsequent control decisions. Especially when processing high-frequency sampling data in large-scale modular hydrogen production systems, it can effectively reduce the burden of data transmission and processing while maintaining the integrity of key feature information.

[0029] According to one aspect of this application, step S2 specifically comprises:

[0030] Step S21: Input the denoised wavelet coefficients W' and eigenvector F into the adaptive state expansion module to construct the expanded state vector x, which includes node voltage magnitude, phase angle and system frequency; apply adaptive weighted iterative unscented transformation to the expanded state vector x to generate 2n+1 sigma point sets X_i and corresponding weights W_i; input X_i into the nonlinear state prediction equation, and calculate the predicted state x(k+1|k) by combining the pre-stored process noise statistical characteristics; input x(k+1|k) into the nonlinear measurement prediction equation, and calculate the predicted observation y(k+1|k) by combining the observation noise covariance matrix R; update the state estimate through sigma point inverse transformation and weighted summation to obtain the system state vector x_est;

[0031] Step S22: Extract the frequency component f_est from the system state vector x_est; input f_est into the variational mode decomposition module, and obtain a series of intrinsic mode functions IMF_k by iteratively optimizing the center frequency ω_k and bandwidth constraint τ_k; apply adaptive multi-scale permutation entropy analysis to each intrinsic mode function IMF_k, calculate the permutation entropy PE_k(τ) at different scales, and extract the entropy value features; input IMF_k into the time-varying fuzzy information entropy module, construct the time-varying fuzzy membership function μ_k(t), and calculate the fuzzy information entropy FIE_k(t); apply Hilbert transform to each intrinsic mode function IMF_k, and calculate the instantaneous frequency IF_k(t) and instantaneous amplitude IA_k(t); integrate all features using a multi-layer feature fusion network to obtain the frequency perturbation feature matrix F_dist;

[0032] Step S23: Input the frequency disturbance feature matrix F_dist and the system state vector x_est into the power balance calculation module. Combine the pre-stored power generation P_gen and load power P_load data to calculate the power imbalance ΔP. Extract the frequency change rate data from the frequency disturbance feature matrix F_dist and apply the adaptive time window selection algorithm to determine the optimal analysis window [t1,t2]. Within the selected window, use the inertia estimator to calculate the system inertia and obtain the preliminary inertia estimate M_raw. Input the preliminary inertia estimate M_raw into the dynamic correction module. Combine the historical inertia estimation results and the system operating status, use the recursive Bayes estimation method to correct it and obtain the grid inertia estimate M_est.

[0033] By combining unscented Kalman filtering, variational mode decomposition, and inertia estimation, accurate estimation of the power grid state and extraction of frequency disturbance features were achieved. Specifically, firstly, adaptive state dimension expansion and unscented transformation were used to generate the optimal sigma point set and weights, improving the accuracy of state estimation. Secondly, variational mode decomposition was employed to achieve accurate decomposition of frequency components through iterative optimization of the center frequency and bandwidth constraints. Thirdly, multi-dimensional features of frequency disturbances were extracted through multi-scale permutation entropy analysis and time-varying fuzzy information entropy. Finally, an inertia estimator and dynamic correction module were used to accurately estimate the power grid inertia. This multi-level state estimation and feature extraction scheme significantly improves the accuracy and real-time performance of system state perception. Especially when dealing with large-scale hydrogen production systems participating in power grid frequency regulation, it can accurately capture the dynamic characteristics and disturbance features of the system, providing reliable state information for subsequent control strategy optimization. This scheme not only improves the accuracy of state estimation but also enhances the system's ability to perceive disturbances, enabling the hydrogen production system to better respond to the power grid's frequency regulation requirements.

[0034] According to one aspect of this application, step S3 specifically comprises:

[0035] Step S31: Input the electrolytic cell data stream E(t) into the data segmentation module, and use the dynamic programming algorithm to perform optimal segmentation to obtain the segmented dataset E_seg; apply a polynomial kernel function mapping to each segment of the segmented dataset E_seg to construct a nonlinear static input mapping f(u(t)); input the mapped data into the linear dynamic system identification module, and obtain the linear dynamic block transfer function G(s) through the subspace identification algorithm; apply the generalized Hammerstein structure mapping to the system output to construct a nonlinear output mapping h(x(t)); combine f(u(t)), G(s) and h(x(t)) in series to obtain the electrolytic cell dynamic model structure HW;

[0036] Step S32: Input the model structure HW and the electrolytic cell data stream E(t) into the dual recursive least squares module; construct the augmented state vector x_aug, which contains the system state and the parameters to be identified; generate multiple sets of Sigma points {X_i} and calculate the corresponding weights W_i; apply the state prediction equation f(X_i) and the observation prediction equation h(X_i) to each Sigma point; calculate the predicted state mean x_pred and the predicted observation mean y_pred; calculate the Kalman gain K using the new observation data and update the state estimate x_update; extract the parameter estimates from x_update to obtain the model parameter vector θ_est;

[0037] Step S33: Input the parameter vector θ_est into the adaptive moment estimation module to construct a sliding time window W(k); calculate the higher-order statistical moments M_1 to M_4 of the data within the window; adaptively adjust the forgetting factor λ(k) according to the trend of moment changes; update the moment estimate M_update and calculate the parameter correction Δθ; input the corrected parameters into the multi-model adaptive estimator to design N candidate models {M_1, M_2, ..., M_N}; calculate the likelihood probability p(y(k)|M_i) of each model and update the model weights w_i(k); optimize the parameters using particle filtering to obtain the final model parameter vector θ_final.

[0038] Step S34: Input the final parameter vector θ_final and the electrolytic cell validation dataset E_val(t) into the model evaluation module; construct the parameter constraint set Θ and design a constraint optimization problem based on projection; use a quadratic programming solver to obtain the parameter estimates θ_constrained that satisfy the constraints; input θ_constrained into the online model validator and calculate the predicted output y_pred(t); compare y_pred(t) with the actual output y(t) to evaluate the model performance index J_perf; based on the performance evaluation results, use an iterative optimization algorithm to fine-tune the parameters and output the final optimized model parameter vector θ_opt.

[0039] Accurate modeling of the dynamic characteristics of electrolyzers was achieved using the Hammerstein-Wiener model and multiple recursive least squares method. Specific innovations include: First, optimal data segmentation was performed using dynamic programming, and nonlinear static mapping was conducted using a composite kernel function, improving the model's nonlinear expressive power. Second, generalized matrix fraction decomposition and subspace identification algorithms were employed to accurately identify linear dynamic characteristics. Third, real-time updates of model parameters were achieved through double recursive least squares method and adaptive moment estimation. Finally, robustness of parameter estimation was improved by utilizing multi-model adaptive estimation and particle filter optimization. This comprehensive modeling scheme not only accurately describes the static and dynamic characteristics of electrolyzers but also enables real-time parameter updates and optimization. Particularly in large-scale hydrogen production systems, it can accurately capture dynamic characteristic changes under different operating conditions, providing precise model support for the system's participation in grid frequency regulation. By considering multi-model fusion and parameter constraints, the adaptability and reliability of the model are significantly improved, enabling the hydrogen production system to better adapt to the dynamic demands of grid frequency regulation.

[0040] According to one aspect of this application, step S4 specifically comprises:

[0041] Step S41: Input the system state vector x_est, frequency disturbance feature matrix F_dist, and grid inertia estimate M_est into the hierarchical control structure design module; construct a fast response layer controller and output fast control command u_f; design a medium speed regulation layer controller and output medium speed control command u_m; configure a slow speed optimization layer controller and output slow speed control command u_s; input the three-layer control commands into the coordinated control module, apply the time-sharing hierarchical optimization algorithm, and obtain the multi-time-scale control framework structure MSC; wherein the response time of the fast response layer controller is 0.1-1s; the response time of the medium speed regulation layer controller is 1-10s; and the response time of the slow speed optimization layer controller is 10-60s.

[0042] Step S42: Read the historical VSG parameter dataset D_hist from the pre-stored database; construct the input feature vector x using the grid inertia estimate M_est and the frequency disturbance feature matrix F_dist; apply a deep kernel network to the feature vector x to perform feature transformation, obtaining the transformed feature x_trans; input x_trans into the optimal kernel function selection module, and select the optimal kernel function k_opt from the candidate kernel function set K using the marginal likelihood maximization criterion; calculate the covariance matrix K_opt and the eigenvalue decomposition result {U, Λ}; perform parameter prediction using sparse Gaussian process regression, select the induced point set Z and calculate the cross-covariance matrix K_xz; update the VSG parameters using the adaptive particle swarm optimization algorithm to obtain the virtual rotational inertia J_vsg and the virtual damping coefficient D_vsg.

[0043] Step S43: Input the control frame MSC, VSG parameters [J_vsg, D_vsg], and model parameter vector θ_final into the dynamic scene tree generation module; construct the frequency perturbation probability distribution p(w) using the kernel density estimation method; generate the scene tree S_tree using a hierarchical sampling strategy; perform a scene reduction algorithm on S_tree to obtain the simplified scene set S_reduced and the probability set Π_reduced; construct a robust stochastic model predictive control problem, design the objective function J and constraint condition g; solve the optimization problem using the distributed alternating direction multiplier method to obtain the optimal control sequence u*.

[0044] Step S44: Input the optimal control sequence u into the adaptive learning module; construct the value function approximator V(x, θ_V) and the policy function approximator π(x, θ_π); update the value function and policy function parameters using the deep deterministic policy gradient algorithm; establish an experience replay buffer D to store transfer samples; sample mini-batch from D for policy update; fuse the updated policy function π(x, θ_π) with u to design an adaptive weight function w(x); calculate the fused control input u_fusion; apply input constraint processing to u_fusion to obtain the final control input u_final;

[0045] Step S45: Read historical control effect data H_ctrl from the pre-stored database; input u_final, the output of step S44, into the long-term optimization strategy generation module; construct the decision set D and the disturbance set Ω; design the minimax regret learning objective function and update the decision variables using the online gradient descent method; adopt the adaptive learning rate adjustment strategy η(t); input the optimization result into the strategy evaluation module and calculate the long-term performance index J_long; dynamically adjust the strategy based on J_long and output the final long-term frequency modulation strategy L_opt.

[0046] A multi-timescale control framework and online Gaussian process regression were employed to achieve coordinated control of the hydrogen production system's participation in grid frequency regulation. The specific innovations are as follows: First, a hierarchical control structure comprising a fast response layer, a medium-speed regulation layer, and a slow optimization layer was constructed to achieve coordinated control across different timescales. Second, adaptive optimization of virtual synchronous generator parameters was achieved through deep kernel networks and online Gaussian process regression. Third, stochastic model predictive control based on scenario trees improved the robustness of the control strategy. Finally, optimization of the long-term frequency regulation strategy was achieved through minimum-maximum regret learning. This multi-timescale coordinated control scheme not only enables the hydrogen production system to respond quickly to grid frequencies but also ensures the system's economic efficiency and reliability. Especially when dealing with large-scale hydrogen production systems, it effectively coordinates the operation of multiple subsystems, optimizing the overall system performance. By considering the control requirements of multiple timescales, the system's frequency regulation performance and operating efficiency are significantly improved.

[0047] According to one aspect of this application, step S5 specifically comprises:

[0048] Step S51: Input the control strategies MSC, J_vsg, D_vsg, u*, L_opt, and the submodule current data I(t) into the distributed coordination control module; read the communication topology information G and the initial weight matrix W from the pre-stored database; construct the local state vector x_i, which includes the inductor current I_i, the output voltage V_i, and the local frequency f_i; input the weight matrix W into the adaptive weight optimizer, apply the gradient projection method to update the communication weights, and obtain the optimized weight matrix W_opt; calculate the Laplacian matrix L, input L into the distributed H∞ controller, and obtain the local control law u_i.

[0049] Step S52: Input the local control law u_i into the event triggering mechanism design module; construct the trigger function f_i(e_i, x_i), calculate the relationship between the measurement error e_i and the state x_i; dynamically adjust the trigger threshold σ_i(t) according to the designed adaptive trigger threshold update rule; construct the Zeno behavior exclusion mechanism to ensure the minimum trigger interval τ_min; input the trigger result into the robustness analysis module, use the polynomial chaotic expansion method to handle parameter uncertainties, and obtain the extended stability region R_extended; based on the trigger result and stability analysis, generate the distributed model predictive control problem P_i, and apply the alternating direction multiplier method to solve for the local optimal control sequence u_i*;

[0050] Step S53: Input the local optimal control sequence {u_i*} into the adaptive multi-scale evaluation module; apply wavelet packet transform to each performance index J_i for multi-scale decomposition to obtain performance components {J_i,k} in different frequency bands; calculate the weights w_i,k of each scale to construct the comprehensive performance index J_comp_i; input J_comp_i into the dynamic Pareto front generation module; initialize the population P(0) using the NSGA-III algorithm, and design a selection operation based on reference points; apply simulated binary crossover and polynomial mutation operators to generate offspring Q(t); select the next generation population through the reference point association mechanism to obtain the Pareto optimal solution set P*;

[0051] Step S54: Input the Pareto optimal solution set P* into the adaptive weighted Kriging surrogate model; generate an initial sample point set X using Latin hypercube sampling; evaluate the objective function value Y of the sample points; construct the Kriging model to predict the function value and prediction variance of the unknown points; design the expected improvement criterion EI(x); dynamically adjust the EI weight w_EI using an adaptive weight strategy; solve the weighted EI maximization problem to obtain new sample points x_new; update the sample set and retrain the model to obtain the final Kriging surrogate model M_Kriging;

[0052] Step S55: Input the surrogate model M_Kriging into the multi-objective Bayesian optimization module; design the acquisition function α(x) and use the Thompson sampling strategy to sample from the posterior distribution; solve the sub-optimal problem to obtain the optimal point x_opt; evaluate the actual objective function value y_opt; update the Gaussian process model; apply the Sobol method to calculate the global sensitivity index {S_i, ST_i}; construct a response surface model to analyze the impact of decision variables; calculate the robustness index RI and confidence interval CI; generate a comprehensive performance evaluation report E_report, which includes the multi-objective evaluation results, Pareto optimal solution, sensitivity analysis and optimization suggestions, and outputs the optimized control parameter set P_opt.

[0053] This study achieves coordinated control and performance optimization of a large-scale hydrogen production system through distributed consensus algorithms and multi-objective performance evaluation. Specific innovations include: first, improving the convergence and robustness of distributed control by employing graph analysis and weighted optimization; second, reducing communication burden through event-triggered mechanisms and Zeno behavior exclusion; third, realizing multi-objective coordinated optimization using multi-scale performance evaluation and Pareto optimization; and finally, improving optimization efficiency through an adaptive weighted Kriging surrogate model and multi-objective Bayesian optimization. This distributed coordinated control scheme not only achieves efficient coordinated operation of the large-scale hydrogen production system but also ensures continuous optimization of system performance. Especially when handling a large number of distributed subsystems, it effectively reduces communication burden and improves system scalability. Through multi-objective optimization and performance evaluation, a balance between system operating efficiency and frequency regulation performance is achieved, providing reliable control support for large-scale hydrogen production systems to participate in grid frequency regulation.

[0054] According to one aspect of this application, step S11 specifically comprises:

[0055] Step S111: Obtain the original voltage, current, and phase angle data V_raw(t), I_raw(t), and θ_raw(t) from the high-precision synchronous phasor measurement unit; input these data into the information entropy evaluator to calculate the information entropy H(t) for different time windows; calculate the signal fluctuation index σ(t) using sliding variance analysis; construct an adaptive sampling rate adjustment function f_s(t) based on H(t) and σ(t); apply a variable sampling rate reconstruction algorithm to the original data to obtain the resampled data V_s(t), I_s(t), and θ_s(t);

[0056] Step S112: Construct the observation vector y(t) from the output of step S111 using V_s(t), I_s(t), and θ_s(t); read the initial state estimate x(0) and covariance matrix P(0) from the pre-stored parameter library; construct the noise covariance matrix Q(t) and R(t) using the adaptive attenuation factor α(t); apply the unscented transformation to generate the sigma point set χ_i and weight W_i; perform state prediction and observation prediction for each sigma point; update the gain K(t) by calculating the cross covariance; obtain the pre-filtered PMU data stream P(t).

[0057] Step S113: Collect raw operating parameters from the electrolytic cell, including voltage V_e(t), current I_e(t), temperature T_e(t), pressure P_e(t), and electrolyte concentration C_e(t); apply wavelet dictionary learning algorithm to construct sparse representation basis D_i for each parameter; calculate compression coefficient α_i using orthogonal matching pursuit algorithm; design an adaptive compression ratio controller to dynamically adjust the compression ratio r_i(t) according to the signal importance index; perform distributed compression sampling to obtain the compressed electrolytic cell data stream E(t);

[0058] Step S114: Input the pre-filtered PMU data stream P(t) and the compressed electrolytic cell data stream E(t) into the multi-dimensional data quality assessment module; collect inductor current I_L(t) from 20 sub-modules to construct a current data matrix I(t); calculate the integrity index C_i, accuracy index A_i, and consistency index U_i for P(t), E(t), and I(t) respectively; apply the fuzzy comprehensive evaluation method and combine it with the pre-stored evaluation rule set R to calculate the comprehensive quality score Q_i; perform weighted processing on the data according to Q_i and output the data stream set D_q(t) with quality score.

[0059] First, an adaptive sampling rate adjustment function is constructed using information entropy assessment and sliding variance analysis to dynamically optimize the PMU data sampling rate. This ensures higher sampling accuracy during rapid system state changes and reduces the sampling rate to decrease data burden when the system is stable. Second, a multivariate composite Kalman pre-filter is used for initial noise reduction of the electrolyzer data. This filter dynamically adjusts the noise covariance matrix through an adaptive attenuation factor, significantly improving the signal-to-noise ratio of the electrolyzer operating parameters. Finally, a distributed data compression sampling module is used to compress the submodule current data. Efficient data compression is achieved through wavelet dictionary learning and orthogonal matching pursuit algorithms, significantly reducing data transmission burden while preserving key information. This multi-level data acquisition and preprocessing scheme not only improves data quality but also optimizes data transmission efficiency, making it particularly suitable for the multi-source heterogeneous data processing needs of large-scale hydrogen production systems. In particular, the multi-dimensional data quality assessment module calculates integrity, accuracy, and consistency indicators, enabling quantitative assessment of data quality and providing a reliable data foundation for subsequent processing.

[0060] According to one aspect of this application, step S12 specifically comprises:

[0061] Step S121: Decompose the data stream set D_q(t) output in step S11 into PMU data stream P_q(t), electrolytic cell data stream E_q(t), and submodule current data stream I_q(t); read the reference clock source T_ref and delay compensation parameter set Δt from the pre-stored time synchronization parameter library; apply a bidirectional delay measurement algorithm to each data stream to calculate the transmission delay τ_i and clock deviation δ_i; use an adaptive Kalman filter to estimate and correct τ_i and δ_i in real time; based on the corrected time parameters, perform time alignment interpolation on the three types of data streams to obtain the preliminary aligned data stream D_sync(t);

[0062] Step S122: Reconstruct the D_sync(t) output in step S121 into a three-dimensional tensor X∈R^{I×J×K}; apply Tucker decomposition to X to obtain the core tensor G and factor matrices U, V, W; calculate the tensor reconstruction error E = X - G×_1U×_2V×_3W; calculate the anomaly detection threshold λ(t) using an adaptive threshold algorithm; reconstruct the tensor decomposition problem through sparse constraint optimization to obtain the low-rank tensor L and sparse tensor S; extract the set of abnormal data points A from S, and output the normal data tensor L and the abnormal dataset A.

[0063] Step S123: Input the abnormal dataset A output in step S122 into a multilayer autoencoder network; load the pre-trained encoder parameters θ_e and decoder parameters θ_d from the pre-stored model library; perform feature encoding on the data points in A to obtain the hidden layer representation h = f_e(A;θ_e); select similar patterns as references from the normal data tensor L using an attention mechanism; combine the reference patterns and the hidden layer representations to generate repair data R = f_d(h;θ_d) through the decoder; apply residual learning to optimize the repair results to obtain the final repair dataset R_opt;

[0064] Step S124: Take the normal data tensor L output in step S122 and the repaired dataset R_opt output in step S123 as input; construct a k-nearest neighbor graph G_k and calculate the local geometric structure of each data point; design an adaptive weight matrix W_l to quantify the local relationships between data points; construct a global alignment objective function J = tr(X^TLX) using the manifold regularization algorithm, where L is the graph Laplacian matrix; solve the alignment problem by alternating optimization to obtain the aligned data stream D'(t); apply adaptive multi-scale filtering to D'(t) and output the filtered data stream D_f(t).

[0065] First, a bidirectional delay measurement algorithm is used to accurately calculate transmission delay and clock skew, and an adaptive Kalman filter is used for real-time correction, ensuring microsecond-level time synchronization accuracy. Second, the Tucker decomposition method is used to reconstruct and decompose the synchronized data into tensors, and an adaptive threshold algorithm is used to accurately identify abnormal data, significantly improving the accuracy of anomaly detection. Third, a multi-layer autoencoder network combined with an attention mechanism is used to repair abnormal data. This method selects similar patterns from normal data as references, achieving high-quality repair of abnormal data. Finally, by constructing a k-nearest neighbor graph and designing an adaptive weight matrix, local manifold alignment of the data is achieved, ensuring the spatiotemporal consistency of the repaired data. This comprehensive data synchronization and repair scheme not only ensures the time alignment accuracy of multi-source data in large-scale hydrogen production systems but also improves data reliability and continuity. Especially when dealing with high-sampling-rate current data and low-sampling-rate electrolyzer data, it can effectively coordinate data synchronization issues at different time scales, providing high-quality data support for the coordinated control of the system.

[0066] According to one aspect of this application, step S21 specifically comprises:

[0067] Step S211: Input the denoised wavelet coefficients W' and eigenvector F into the adaptive state expansion module; construct the expanded state vector x, which includes node voltage amplitude, phase angle and system frequency; read the state expansion parameter set P_ext from the pre-stored parameter library; calculate the state correlation matrix C using the Cramer decomposition method; dynamically adjust the expansion strategy based on the correlation matrix C to obtain the optimized expanded state vector x_opt; input the state vector x_opt into the covariance estimator to calculate the initial state covariance matrix P(0);

[0068] Step S212: Input the extended state vector x_opt into the spherical-radial sampling module; calculate the optimal sampling scaling factor κ, and use the Householder transformation to generate 2n+1 uniformly distributed sigma point sets X_i; design an adaptive weight allocation strategy based on state sensitivity, and calculate the weight coefficient W_i for each sampling point; read the process noise statistical characteristics Q from the pre-stored parameter library, and construct the process noise covariance matrix Q_new that considers state correlation; input the sigma point set {X_i} and the weights {W_i} into the adaptive sampling point optimizer to obtain the optimized sampling point set {X_i_opt};

[0069] Step S213: Input the optimized sampling point set {X_i_opt} into the nonlinear state prediction module; read the state transition parameter θ_t from the pre-stored model library; apply an odorless transformation to each sampling point and calculate the one-step predicted state X_pred_i; design an anomaly prediction detector based on Mahalanobis distance to identify and correct anomaly prediction points; use the second-order moment matching method to calculate the predicted state mean x(k+1|k) and prediction error covariance P(k+1|k);

[0070] Step S214: Input the predicted state x(k+1|k) into the nonlinear measurement prediction module; read the observation noise covariance matrix R from the pre-stored parameter library; construct the observation prediction equation based on ellipsoidal constraints and calculate the predicted observation y(k+1|k); design an adaptive measurement noise estimator based on maximum likelihood estimation and dynamically update the R matrix; calculate the cross-covariance matrix P_xy between the predicted state and the observation using the singular value decomposition method.

[0071] Step S215: Input the predicted observation y(k+1|k) and the actual observation y(k+1) into the gain calculation module; calculate the inverse matrix of the innovation covariance based on the Cholesky decomposition method; design a gain matrix calculation method with fading memory characteristics to obtain the optimal filter gain K(k+1); update the state components one by one using sequential processing technology to obtain the state estimate x_est(k+1); perform consistency verification between the state estimate result and the prediction result to determine the reliability of the estimate; output the system state vector x_est.

[0072] First, the state correlation matrix is ​​calculated using the Cramer decomposition method. Based on this correlation matrix, the dimensionality expansion strategy is dynamically adjusted, improving the dimensionality adaptability of the state estimation. Second, a uniformly distributed sigma point set is generated using a spherical-radial sampling module, and the sampling effect is optimized through an adaptive weight allocation strategy based on state sensitivity. Third, a tasteless transformation is used for state prediction, and the reliability of the prediction is improved by combining it with an anomaly prediction detector based on Mahalanobis distance. Finally, a gain matrix calculation method with fading memory characteristics is used to achieve stable updates of the state estimation. This innovative state estimation scheme not only improves estimation accuracy but also enhances the robustness of the algorithm. Especially when dealing with multivariate state estimation in large-scale hydrogen production systems, it effectively handles the coupling relationships between state variables, providing accurate state information for system control. Sequential processing techniques and consistency checks further ensure the reliability of the state estimation results, enabling the hydrogen production system to respond more accurately to grid frequency regulation requirements.

[0073] According to one aspect of this application, step S22 specifically comprises:

[0074] Step S221: Input the frequency component f_est in the system state vector x_est into the adaptive bandwidth optimization module; construct the variational mode decomposition objective function, which includes the center frequency ω_k and the bandwidth constraint τ_k; use the alternating direction multiplier method to iteratively solve the optimization problem, and update the Lagrange multiplier μ and the penalty factor α for each iteration; design an adaptive mode number selection algorithm to dynamically adjust the number of decomposition layers K according to the signal energy distribution; perform variational mode decomposition on f_est to obtain a series of mode functions {u_1,u_2,...,u_K} and the corresponding center frequencies {ω_1,ω_2,...,ω_K}.

[0075] Step S222: Input the modal function set {u_k} output in step S221 into the multi-scale analysis module; construct coarse-grained time series {x_τ} at different scales for each modality u_k; design an adaptive embedding dimension selector to determine the optimal embedding dimension m_k and time delay τ_k based on the mutual information criterion; calculate the permutation pattern probability distribution p_i at each scale; calculate the permutation entropy PE_k(τ) using the information entropy formula; extract the permutation entropy curve features, including the entropy slope s_k and curvature c_k; output the permutation entropy feature vector PE_features.

[0076] Step S223: Input the modal function {u_k} from step S221 and the PE_features from step S222 into the time-frequency analysis module; apply the generalized S-transform to each modality to calculate the time-frequency distribution TFD_k; construct an adaptive time-frequency resolution controller to dynamically adjust the time-frequency window parameters; extract time-frequency features including instantaneous frequency IF_k(t), instantaneous amplitude IA_k(t), and energy distribution ED_k(t); design a multi-feature fusion network, including a feature mapping layer, an attention layer, and a fully connected layer; obtain the fused feature vector F_fused through network forward propagation.

[0077] Step S224: Input the fused feature vector F_fused output from step S223 into the deep learning module; construct a deep neural network containing multiple residual blocks; load pre-trained parameters θ_pre from the pre-stored model library; process the feature data using local response normalization; extract multi-level features through residual connections and skip connections; design a feature selection gating unit to dynamically filter important features; apply dropout regularization to prevent overfitting; perform principal component analysis to reduce the dimensionality of the network output and retain the main components; organize the final frequency perturbation feature matrix F_dist.

[0078] First, the center frequency and bandwidth constraints are dynamically adjusted by an adaptive bandwidth optimization module, improving the accuracy of mode decomposition. Second, the complexity features of the frequency signal are extracted using multi-scale permutation entropy analysis, and the accuracy of feature extraction is ensured by an adaptive embedding dimension selector. Third, time-frequency analysis is performed using a generalized S-transform, combined with an adaptive time-frequency resolution controller to achieve accurate extraction of time-frequency features. Finally, effective fusion of multi-level features is achieved through deep learning networks and feature selection gating units. This comprehensive frequency feature extraction scheme can not only accurately capture the frequency response characteristics of the hydrogen production system at different time scales, but also identify the dynamic behavior patterns of the system during frequency regulation. Especially when dealing with large-scale hydrogen production systems participating in grid frequency regulation, it can accurately analyze the system's frequency response capability and dynamic characteristics, providing rich feature information for subsequent control optimization. Through multi-feature fusion and dimensionality reduction, the effectiveness of feature representation and computational efficiency are significantly improved.

[0079] According to one aspect of this application, step S31 specifically comprises:

[0080] Step S311: Input the electrolytic cell data stream E(t) into the data segmentation module; read the segmentation parameter set S_param from the pre-stored parameter library; perform multi-scale analysis on the data using wavelet transform to identify data feature points; construct a segmentation boundary detector based on the Kolmogorov-Smirnov test; design a dynamic programming algorithm for optimal segmentation, and use the BIC criterion to determine the optimal number of segments; apply local linear regression to each data segment for preprocessing to obtain the segmented dataset E_seg;

[0081] Step S312: Input the segmented dataset E_seg into the nonlinear static mapping module; read the kernel function parameter set K_param from the pre-stored model library; construct a composite kernel function based on Mercer's theorem, including radial basis kernel and polynomial kernel; perform feature extraction using kernel principal component analysis to obtain the kernel feature mapping φ(u); design a sparse kernel mapping method based on L1 regularization to construct the nonlinear static input mapping f(u(t));

[0082] Step S313: Input the mapped data f(u(t)) into the linear dynamic system identification module; construct the initial structure of the system transfer function using the generalized matrix fraction decomposition method; estimate the system state space matrix using the subspace identification algorithm; design a model order selection criterion based on the Hankel matrix; and perform numerical calculations using the QR decomposition algorithm to obtain the linear dynamic block transfer function G(s).

[0083] Step S314: Compare the linear dynamic response output with the original data; construct a nonlinear output mapping based on the generalized Hammerstein structure; approximate the nonlinear characteristics using a piecewise linear approximator; design smoothing constraints based on B-spline functions to construct the nonlinear output mapping h(x(t)); and combine f(u(t)), G(s), and h(x(t)) in series to obtain the electrolytic cell dynamic model structure HW.

[0084] First, a piecewise boundary detector based on the Kolmogorov-Smirnov test and a dynamic programming algorithm were employed to achieve optimal data segmentation, improving the local accuracy of the modeling. Second, by constructing a composite kernel function based on Mercer's theorem and kernel principal component analysis, accurate mapping of nonlinear static characteristics was achieved. Third, generalized matrix fraction decomposition and subspace identification algorithms were used to accurately identify the dynamic characteristics of the system. Finally, a complete nonlinear dynamic model was constructed using generalized Hammerstein structure and B-spline function smoothing constraints. This comprehensive modeling scheme not only accurately describes the static and dynamic characteristics of the electrolyzer but also ensures the model's adaptability under different operating conditions. Especially when dealing with large-scale hydrogen production systems, it can effectively capture the dynamic response characteristics of the system under different operating states, providing reliable model support for the system's participation in grid frequency regulation. The combination of piecewise modeling and nonlinear mapping significantly improves the model's local accuracy and global adaptability.

[0085] According to one aspect of this application, step S32 specifically comprises:

[0086] Step S321: Input the model structure HW and the electrolytic cell data stream E(t) into the augmented state construction module; design the state vector x containing physical quantities such as voltage and current, and construct the parameter vector θ containing nonlinear mapping coefficients and dynamic parameters; combine x and θ to form the augmented state vector x_aug; generate 2n+1 Sigma points {X_i}, and calculate the corresponding weight coefficients {W_m_i,W_c_i}; apply the nonlinear state transition equation f(X_i) to each Sigma point to predict the state at the next time step; calculate the predicted observation value based on the observation equation h(X_i); obtain the state prediction x_pred and the observation prediction y_pred by weighted summation;

[0087] Step S322: Take the predicted state x_pred output from step S321 as input; construct a sliding time window W(t), with the window length L adaptively adjusted according to the signal change rate; calculate the first moment M_1 to the fourth moment M_4 of the data within the window; design a forgetting factor update rule based on the moment change rate to obtain the adaptive forgetting factor λ(t); apply the exponential weighted recursive algorithm to update the moment estimate, M_new = λ(t)M_old + (1-λ(t))M_current; calculate the parameter correction amount Δθ based on the higher-order moment characteristics to obtain the corrected parameter estimate θ_corr;

[0088] Step S323: Input the θ_corr output in step S322 into the multi-model library; construct N candidate models {M_1,M_2,...,M_N} based on historical data and expert knowledge; calculate the likelihood function p(y(k)|M_i) for each model; design a weight update mechanism based on minimizing prediction error; calculate the model weights w_i(k) using the softmax function; construct an interactive multi-model estimation framework to handle the complementarity between models; fuse the parameter estimates of each model to obtain the weighted parameter vector θ_weighted;

[0089] Step S324: Use θ_weighted, output from step S323, as the initial parameter; initialize Np particles {θ_p} and their weights {w_p}; apply the state-space model to predict each particle; calculate the particle likelihood based on the latest observation data; design an adaptive resampling strategy based on the effective number of particles; construct a moving resampling algorithm to avoid particle degradation; update and normalize the particle weights; calculate the mean and variance of the posterior distribution; and obtain the final parameter vector θ_final through iterative optimization.

[0090] First, an augmented state vector containing physical states and parameters to be identified is constructed, and the accuracy of parameter estimation is improved by generating an optimal Sigma point set. Second, nonlinear state prediction and observation prediction equations, combined with adaptive window length and forgetting factor, enable real-time parameter updates. Third, higher-order moment estimation and adaptive weighting strategies improve the robustness of parameter estimation. Finally, likelihood calculation and model weight updates enable parameter optimization under multi-model fusion. This innovative parameter identification scheme not only improves the accuracy of parameter estimation but also enhances the algorithm's adaptability. Especially when dealing with parameter identification in large-scale hydrogen production systems, it effectively addresses the time-varying characteristics and uncertainties of system parameters, providing reliable parameter support for precise system control. Multi-model fusion and adaptive moment estimation significantly improve the robustness and real-time performance of parameter identification.

[0091] According to one aspect of this application, step S33 specifically comprises:

[0092] Step S331: Input the parameter vector θ_est into the adaptive moment estimation module; read the initial window parameter W_init from the pre-stored parameter library; construct a dynamic window size determination criterion using generalized extremum theory; calculate the higher-order statistical moments M_1 to M_4 of the data within the window based on local polynomial fitting; design a window update mechanism based on the minimum description length criterion to dynamically adjust the sliding window W(k); calculate the covariance matrix of the moment estimation using Cholesky decomposition.

[0093] Step S332: Input the statistical moment sequence {M_1,M_2,M_3,M_4} into the forgetting factor optimization module; read the initial value λ_0 of the forgetting factor from the pre-stored parameter library; construct a moment change detector using Kullback-Leibler divergence; calculate the changing trend of the moment sequence based on the exponentially weighted recursive least squares method; design an adaptive threshold based on the generalized likelihood ratio to dynamically adjust the forgetting factor λ(k); update the moment estimate M_update using the Kalman filtering algorithm.

[0094] Step S333: Input the updated moment estimate M_update into the parameter correction module; read the parameter constraint Θ_con from the pre-stored parameter library; calculate the parameter gradient using the quasi-Newton method to obtain the correction amount Δθ; construct a parameter constraint checker based on Lyapunov stability; design a parameter update law with dead zone characteristics to avoid parameter jitter; apply sequential quadratic programming to solve the constrained parameter optimization problem.

[0095] Step S334: Input the corrected parameters into the multi-model adaptive estimator; read N candidate models {M_1,M_2,...,M_N} from the pre-stored model library; calculate the posterior probability p(M_i|y(k)) of each model using the variational Bayesian method; evaluate the model performance based on the generalized predictable error criterion; design a model weight update mechanism based on the entropy weight method and calculate the model weight w_i(k); update the model parameters using maximum likelihood estimation.

[0096] Step S335: Input the multi-model estimation results into the parameter optimization module; construct a parameter sampler based on Markov chain Monte Carlo; generate a parameter sample set using the Metropolis-Hastings algorithm; design a sampling strategy based on adaptive acceptance probability; construct a kernel density estimator for the posterior distribution of the parameters; determine the confidence interval of the parameters using the quantile regression method; output the final model parameter vector θ_final.

[0097] First, higher-order statistical moments are calculated using a sliding time window, and the forgetting factor is dynamically adjusted based on the moment's changing trend, improving the real-time performance of parameter estimation. Second, multiple candidate models are constructed using a multi-model adaptive estimator, and the model weights are dynamically updated through likelihood probability calculation. Third, particle filtering is used for parameter optimization, and a moving resampling algorithm is employed to avoid particle degradation. Finally, parameter constraints and online model validation ensure the physical feasibility of parameter estimation. This comprehensive parameter optimization scheme not only guarantees the accuracy of parameter estimation but also achieves real-time parameter updates and constraint handling. Especially when dealing with parameter optimization in large-scale hydrogen production systems, it effectively handles parameter uncertainties and time-varying characteristics, providing reliable parameter support for stable system operation. Through multi-model fusion and constraint optimization, the reliability and adaptability of parameter estimation are significantly improved.

[0098] According to one aspect of this application, step S41 specifically comprises:

[0099] Step S411: Input the system state vector x_est, frequency disturbance feature matrix F_dist, and grid inertia estimate M_est into the hierarchical control structure design module; read the control hierarchy parameters set L_param from the pre-stored parameter library; design the inter-hierarchical decoupling strategy using the optimal separation theorem; evaluate the system controllability based on the singular value decomposition method; design a hierarchical scheduling mechanism based on time scale separation; and construct the fast response layer controller structure C_f.

[0100] Step S412: Input the fast response layer controller C_f into the medium speed regulation layer design module; read the medium speed control parameters M_param from the pre-stored parameter library; design a robust controller using H∞ control theory; construct the controller parameter optimization problem based on linear matrix inequalities; design a feedforward compensator based on generalized predictive control; solve the controller parameters using sequential quadratic programming; construct the medium speed regulation layer controller structure C_m.

[0101] Step S413: Input the medium-speed regulation layer controller C_m into the slow-speed optimization layer design module; read the slow-speed control parameter S_param from the pre-stored parameter library; construct the optimal control problem using Pontryagin's maximum principle; design an economic scheduling strategy based on dynamic programming algorithm; design a constraint handling mechanism based on model prediction; solve the optimization problem using the interior point method; construct the slow-speed optimization layer controller structure C_s.

[0102] Step S414: Input the three-layer controller {C_f,C_m,C_s} into the coordination control module; read the coordination control parameter C_param from the pre-stored parameter library; design a hierarchical coordination mechanism using Lyapunov stability theory; construct the optimal coordination strategy based on the Hamilton-Jacobi equation; design a control quantity allocation algorithm based on adaptive weights; solve the coordination control problem using the gradient projection method; output the multi-timescale control framework structure MSC.

[0103] First, a hierarchical control structure comprising a fast response layer, a medium-speed regulation layer, and a slow-speed optimization layer was constructed, achieving control coordination across different time scales. Second, a controller for the medium-speed regulation layer was designed using H∞ control theory and generalized predictive control, improving the system's robustness. Third, the economic dispatch strategy for the slow-speed layer was optimized using Pontryagin's maximum principle and dynamic programming algorithms. Finally, time-sharing hierarchical optimization achieved coordination of the multi-layered control. This multi-time-scale control scheme not only meets the rapid response requirements of the hydrogen production system participating in grid frequency regulation but also ensures the economic efficiency of system operation. Especially when dealing with large-scale hydrogen production systems, it effectively coordinates control objectives across different time scales, achieving overall system performance optimization. Through hierarchical control and coordinated optimization, the system's control performance and operational efficiency are significantly improved.

[0104] According to one aspect of this application, step S42 specifically comprises:

[0105] Step S421: Read the historical VSG parameter dataset D_hist from the pre-stored database; combine the power grid inertia estimate M_est and the frequency disturbance feature matrix F_dist to construct the input feature vector x; design a multi-layer deep kernel network, including a kernel feature extraction layer and a feature transformation layer; input x into the first-layer kernel function to generate the initial kernel feature mapping φ_1(x); obtain the deep kernel feature φ_d(x) through multi-layer nonlinear transformation and residual connection; apply the attention mechanism to weight the features to obtain the representation x_trans that highlights the key features; optimize the network parameters using the backpropagation algorithm and output the transformed feature vector x_opt.

[0106] Step S422: Input the feature vector x_opt output in step S421 into the kernel function selection module; construct a candidate kernel function set K={k_SE,k_Matern,k_RQ,k_Linear}; evaluate the performance of each kernel function using the marginal likelihood maximization criterion; select the optimal kernel function k_opt to construct the basic kernel matrix; apply the Nyström method to approximate the large-scale kernel matrix; select m induced points Z={z_1,...,z_m}; calculate the cross-covariance matrices K_xz and K_zz; construct the sparse approximate prediction equation to obtain the VSG parameter prediction values ​​[J_pred,D_pred] and prediction variance [v_J,v_D].

[0107] Step S423: Use the parameter prediction values ​​[J_pred, D_pred] output in step S422 as the initial solution; initialize the particle swarm P={p_1,...,p_N}; design an adaptive inertia weight update strategy w(t); construct a comprehensive performance evaluation function f(J,D), considering frequency deviation, rate of change, and kinetic energy loss; evaluate the fitness value f(p_i) for each particle; update the individual optimal position p_best_i and the global optimal position g_best; dynamically adjust the learning factors c1 and c2 according to the performance gradient; apply the velocity update formula to avoid local optima; iteratively optimize to obtain the parameter optimization values ​​[J_opt, D_opt].

[0108] Step S424: Input the optimized parameters [J_opt, D_opt] output in step S423 into the robustness analysis module; construct the state-space model of the VSG control system; define the parameter disturbance ranges ΔJ and ΔD; apply structural singular value analysis to evaluate the system stability margin; calculate the H∞ norm and stability index of the system; design a Monte Carlo simulation scheme to verify parameter sensitivity; construct a parameter adaptive compensator and generate compensation quantities [ΔJ_comp, ΔD_comp]; verify the stability of the compensated parameters; output the final VSG parameters [J_vsg, D_vsg] and the stability evaluation report R_stable.

[0109] First, feature transformation is performed through a multi-layer deep kernel network, combined with an attention mechanism to improve the effectiveness of feature representation. Second, the optimal kernel function is selected using the marginal likelihood maximization criterion, and parameter prediction is achieved through sparse Gaussian process regression. Third, an adaptive particle swarm optimization algorithm is used to update the VSG parameters, improving the efficiency of parameter optimization. Finally, robustness analysis ensures the stability of the parameters. This innovative parameter optimization scheme not only improves the adaptability of VSG parameters but also enhances the dynamic performance of the system. Especially when dealing with large-scale hydrogen production systems, it can accurately capture changes in the dynamic characteristics of the system and achieve real-time optimization of VSG parameters. By combining feature learning and parameter optimization, the frequency modulation performance of the system is significantly improved.

[0110] Step S51 specifically involves:

[0111] Step S511: Read the initial communication topology information G and weight matrix W from the pre-stored database; construct the node connectivity matrix D and adjacency matrix A; calculate the network algebraic connectivity λ_2 as an evaluation index; design the topology optimization objective function J_topo based on graph analysis; search for the optimal topology using the discrete particle swarm optimization algorithm; construct communication quality constraints considering communication delay τ_ij and packet loss rate p_loss; apply the branch and bound method to solve the discrete optimization problem; obtain the optimized topology G_opt and the initial weight matrix W_init.

[0112] Step S512: Take the initial weight matrix W_init output from step S511 as input; construct a local state vector x_i, including inductor current I_i, output voltage V_i, and local frequency f_i; design a weight update gradient based on consistency error; construct a Lyapunov function V(x,W) to evaluate convergence; apply the projection gradient method to update the communication weights, ensuring that W_ij≥0 and the row sum is 1; design an adaptive learning rate η(t) to dynamically adjust the update step size; iterate and optimize until the convergence criterion is met; output the optimized weight matrix W_opt and Laplacian matrix L.

[0113] Step S513: Input the Laplacian matrix L output in step S512 into the controller design module; construct a generalized state-space model, including state equations and output equations; design a distributed state observer to estimate the unknown state ẑ_i; construct the H∞ performance index γ to form a linear matrix inequality (LMI) constraint; solve the LMI optimization problem using the interior-point method to obtain the feedback gain matrix K_i; design the event-triggered observer gain L_i; integrate the state feedback and observer design to obtain the complete controller structure C_i.

[0114] Step S514: Input the controller structure C_i output in step S513 into the distributed optimization module; construct the local cost function J_i and coupling constraints g_i; design a distributed solution framework based on ADMM; introduce augmented Lagrangian quantities to handle coupling constraints; construct a prediction model for each subsystem to predict the state trajectory x_i(k+j|k); iteratively and alternately optimize and update local variables, global variables, and Lagrange multipliers; check convergence and perform error compensation; integrate the optimization results of each subsystem to generate a distributed control instruction set u_i*; output the final local control law u_i.

[0115] First, the communication topology is optimized through graph analysis, improving the convergence of distributed control. Second, an adaptive weight optimizer is used to dynamically update communication weights, enhancing control robustness. Third, local control optimization is achieved through the design of a distributed H∞ controller. Finally, an event-triggered mechanism reduces the communication burden. This distributed coordinated control scheme not only achieves efficient control of large-scale systems but also ensures control reliability. Especially when handling a large number of distributed subsystems, it effectively reduces communication requirements and improves system scalability. Through topology optimization and distributed control, the overall system performance is significantly improved.

[0116] According to one aspect of this application, step S52 specifically comprises:

[0117] Step S521: Input the local control law u_i into the event trigger design module; read the trigger parameter set T_param from the pre-stored parameter library; construct the state-dependent trigger function f_i(e_i,x_i) using the Lyapunov-Krasovskii functional; evaluate the relationship between the measurement error e_i and the state x_i based on the generalized Wirtinger inequality; design the dynamic dead-zone characteristic based on exponential decay; construct a trigger compensation mechanism that considers communication delay.

[0118] Step S522: Input the trigger function f_i(e_i,x_i) into the threshold optimization module; read the initial threshold value σ_0 from the pre-stored parameter library; construct an adaptive threshold update law using a random approximation method; design a dynamic threshold adjustment strategy based on local minima analysis; design a normalized trigger condition based on relative error; apply an extreme value optimization algorithm to solve for the optimal threshold σ_i(t); construct a threshold smoothing mechanism based on exponential moving average.

[0119] Step S523: Input the optimized trigger threshold σ_i(t) into the Zeno behavior analysis module; read the minimum trigger interval parameter τ_min from the pre-stored parameter library; construct the trigger interval constraint using the Dwell-time method; analyze the periodicity of the trigger sequence based on the Poincaré mapping; design a trigger limiter based on the energy function; solve the minimum trigger interval τ_min_new using the convex optimization method; and construct a protection mechanism to prevent frequent triggering.

[0120] Step S524: Input the trigger result into the robustness analysis module; read the uncertainty parameter set U_param from the pre-stored parameter library; construct the uncertainty propagation model using the polynomial chaotic expansion method; analyze the influence of parameter perturbation based on Wiener-Askey expansion; design a stability criterion based on the probability ellipsoid; evaluate the system robustness using the Monte Carlo method; obtain the extended stability region R_extended.

[0121] Step S525: Input the extended stability region R_extended into the distributed optimization module; read the optimization parameter set O_param from the pre-stored parameter library; construct a distributed solution framework using the alternating direction multiplier method; design a parallel optimization algorithm based on dual decomposition theory; design a distributed coordination mechanism based on consistency constraints; apply the proximal point method to solve the subproblem; output the local optimal control sequence u_i*.

[0122] First, the triggering timing of the control signal is optimized through state-dependent trigger function design and adaptive threshold updating. Second, the Zeno behavior exclusion mechanism ensures the minimum triggering interval, avoiding frequent triggering. Third, parameter uncertainty is analyzed using a polynomial chaotic expansion method, expanding the system's stability region. Finally, local optimal control is achieved using distributed model predictive control. This not only reduces the system's communication burden but also ensures the real-time performance and reliability of the control. Dynamic adjustment of the adaptive triggering threshold achieves efficient utilization of communication resources, while robustness analysis ensures stable operation of the system under parameter uncertainty. Especially in handling distributed control of large-scale hydrogen production systems, it effectively balances control performance and communication overhead, providing reliable technical support for large-scale system deployment. The combination of event-triggered mechanisms and distributed optimization significantly improves the system's control efficiency and scalability while ensuring the robustness of the control system in the face of uncertainty. The distributed solution strategy based on the alternating direction multiplier method enables each subsystem to achieve optimal control under limited communication conditions, further enhancing the overall system performance.

[0123] According to one aspect of this application, step S53 specifically comprises:

[0124] Step S531: Input the local optimal control sequence {u_i*} into the performance decomposition module; read the wavelet packet basis function set W_base from the pre-stored parameter library; construct a multi-scale decomposition framework using the wavelet packet transform algorithm; select the optimal decomposition layer based on the Shannon entropy criterion; design a feature extraction strategy based on energy distribution; and calculate the performance components {J_i,k} of different frequency bands using the orthogonal projection method.

[0125] Step S532: Input the performance components {J_i,k} into the weight optimization module; read the initial weight value W_init from the pre-stored parameter library; construct an objective weight calculation model using the entropy weight method; design a subjective weight evaluation criterion based on the analytic hierarchy process; design a dynamic weight adjustment mechanism based on the coefficient of variation; apply the quadratic programming method to solve the weight optimization problem; obtain the weights w_i,k at each scale.

[0126] Step S533: Input the weights w_i,k of each scale into the index fusion module; read the fusion parameter set F_param from the pre-stored parameter library; construct a nonlinear index fusion model using Choquet integral; design an interaction characteristic evaluation method based on fuzzy measure theory; design a dynamic fusion strategy based on synergy; apply a convex optimization algorithm to solve for the fusion weights; obtain the comprehensive performance index J_comp_i.

[0127] Step S534: Input the comprehensive performance index J_comp_i into the Pareto optimization module; read the population parameter set P_param from the pre-stored parameter library; initialize the reference point set using the improved NSGA-III algorithm; design an individual selection strategy based on crowding order; design a genetic operation based on adaptive crossover rate; update the population using a non-dominated sorting method; output the Pareto optimal solution set P*.

[0128] First, a wavelet packet transform algorithm is used for multi-scale decomposition, and the optimal decomposition level is selected using the Shannon entropy criterion, enabling refined analysis of performance indicators. Second, a comprehensive evaluation system of objective and subjective weights is constructed using the entropy weight method and the analytic hierarchy process (AHP), ensuring the scientific nature of weight allocation. Third, the Choquet integral is used to effectively fuse nonlinear indicators, considering the interactions between indicators. Finally, the NSGA-III algorithm is used to solve the multi-objective optimization problem, obtaining the Pareto optimal solution set. This innovative performance evaluation scheme not only achieves comprehensive quantification of system performance but also provides multi-dimensional optimization directions. Especially when evaluating the frequency regulation performance of large-scale hydrogen production systems, it can accurately capture performance characteristics at different time and spatial scales, providing a reliable decision-making basis for system optimization. Through multi-dimensional feature fusion and multi-objective optimization, the accuracy and guidance of performance evaluation are significantly improved.

[0129] According to one aspect of this application, step S54 specifically comprises:

[0130] Step S541: Input the Pareto optimal solution set P* into the sampling design module; read the sampling parameter set S_param from the pre-stored parameter library; construct the initial sampling scheme using Latin hypercube design; optimize the sampling point distribution based on the minimum-maximum distance criterion; design a sampling point supplementation strategy based on space filling; apply cluster analysis to evaluate the representativeness of the samples; obtain the initial sample point set X.

[0131] Step S542: Input the sample point set X into the Kriging modeling module; read the kernel function parameter set K_param from the pre-stored parameter library; identify the Kriging model parameters using the maximum likelihood estimation method; estimate the trend function coefficients based on the generalized least squares method; design the correlation function structure based on anisotropy; solve the correlation matrix using the diagonalization algorithm; construct the Kriging prediction equation.

[0132] Step S543: Input the Kriging model into the adaptive optimization module; read the weight parameter set W_param from the pre-stored parameter library; construct the sampling point evaluation function using the expected improvement criterion; design an exploration-utilization balancing strategy based on the confidence limit ratio; design a local search mechanism based on gradient information; apply the sequential quadratic programming method to solve for the optimal point; obtain the new sampling point x_new;

[0133] Step S544: Input the new sampling point x_new into the model update module; read the observation data y_new from the pre-stored database; update the model parameters using the recursive least squares method; evaluate the prediction accuracy based on the cross-validation method; design an anomaly detection strategy based on residual analysis; apply the incremental learning algorithm to update the surrogate model; output the final Kriging surrogate model M_Kriging.

[0134] By constructing a response surface model and calculating robustness indices, the reliability of the optimization results was assessed, ensuring the effectiveness of the optimization scheme in practical applications.

[0135] In this embodiment, the variational mode decomposition method introduces adaptive bandwidth optimization and dynamic center frequency adjustment mechanisms, enabling the mode decomposition process to automatically adjust decomposition parameters according to signal characteristics. Specifically, an adaptive mode number selection algorithm based on signal energy distribution is designed, and more accurate signal decomposition is achieved by iteratively optimizing the center frequency ω_k and bandwidth constraint τ_k. This is more adaptable to the dynamic changes of frequency signals than the traditional variational mode decomposition method that uses fixed parameters.

[0136] The Hammerstein-Wiener model, based on the traditional HW model, adds dynamic programming for optimal piecewise segmentation and a composite kernel function mapping mechanism. Specifically, it first uses a dynamic programming algorithm to optimally segment the data, and then constructs a composite kernel function by combining radial basis function kernels and polynomial kernels, improving the model's ability to describe the nonlinear characteristics of the electrolyzer. Simultaneously, it employs B-spline functions to implement smoothing constraints, avoiding abrupt changes in the model output.

[0137] A dual recursive least squares method, combined with adaptive moment estimation and a sliding time window, enables real-time parameter updates. Specifically, it constructs an augmented state vector containing the system state and the parameters to be identified, dynamically adjusts the forgetting factor by calculating higher-order statistical moments (first to fourth order moments), and uses the Metropolis-Hastings algorithm to generate parameter samples, thereby improving the accuracy of parameter estimation.

[0138] The multi-timescale control framework divides the control system into a fast response layer (0.1-1s), a medium-speed regulation layer (1-10s), and a slow optimization layer (10-60s), achieving coordinated control across different timescales. Specifically, the fast layer uses H∞ control for rapid response, the medium-speed layer uses generalized predictive control for regulation, and the slow layer performs economic optimization based on Pontryagin's maximum principle, thus achieving global performance optimization.

[0139] An adaptive weighted Kriging surrogate model is proposed, which introduces dynamic weight adjustment and expected improvement criteria to improve the prediction accuracy of the surrogate model. Specifically, Latin hypercube sampling is used to generate initial samples, Kriging model parameters are identified through maximum likelihood estimation, and an adaptive weight strategy is designed to dynamically adjust the expected improvement weights, thereby achieving optimal selection of sampling points.

[0140] Multi-layer deep kernel networks, based on traditional kernel networks, add multi-layer feature extraction and attention mechanisms. Specifically, they design multi-layer nonlinear transformations and residual connection structures, and use attention mechanisms to highlight key features, achieving more effective feature learning and parameter prediction.

[0141] An event-triggered distributed control mechanism, combined with adaptive trigger thresholds and Zeno behavior exclusion mechanisms, optimizes the triggering strategy for control signals. Specifically, it designs state-dependent trigger functions, dynamically adjusts trigger conditions through adaptive thresholds, and ensures a minimum trigger interval, thus avoiding frequent communication issues.

[0142] A distributed multi-objective performance evaluation method, integrating wavelet packet transform and Choquet integral, achieves comprehensive evaluation of multi-dimensional performance indicators. Specifically, it selects the optimal decomposition level using the Shannon entropy criterion, determines the weights using the entropy weight method and the analytic hierarchy process, and considers the interactions between indicators using the Choquet integral, thus achieving a comprehensive performance evaluation.

[0143] The NSGA-III algorithm adds a reference point-based selection operation and an adaptive crossover rate design. Specifically, it selects the next generation population through a reference point association mechanism and dynamically adjusts the crossover and mutation probabilities, improving the convergence and diversity of multi-objective optimization.

[0144] Multi-dimensional tensor decomposition anomaly detection, combined with Tucker decomposition and an adaptive threshold algorithm, improves the accuracy of anomaly detection. Specifically, it constructs a three-dimensional tensor from multi-source data, obtains the core tensor and factor matrix through tensor decomposition, and uses an adaptive threshold algorithm to achieve accurate identification and repair of anomalies.

[0145] In another embodiment of this application, the hardware component is structured as follows:

[0146] The hydrogen production power supply topology mainly includes: a 24-pulse rectifier circuit, a DC-side filter capacitor, 20 sub-modules, an output filter capacitor, and a corresponding control system. The 20 sub-modules are connected in parallel sequentially, then in parallel with the output filter capacitor, and finally connected to the electrolyzer. The control signal generation unit includes: an outer loop subtractor 10, an outer loop proportional-integral controller 11, a current-sharing loop subtractor 12, a current-sharing loop proportional-integral controller 13, an adder 14, a limiter 15, a carrier phase-shift modulator 16, and a comparator 17. One input terminal of the outer loop subtractor 10 receives an inductor current reference signal I. Lref One input terminal of the current sharing ring subtractor 12 receives the average inductor current, and the other input terminal is connected to the input terminal of the outer loop proportional-integral controller 11. The other input terminal of the current sharing ring subtractor 12 receives the average inductor current, and the other input terminal receives the collected inductor current values ​​of each phase. The output terminal is connected to the input terminal of the current sharing ring proportional-integral controller 13. One input terminal of the adder 14 is connected to the output terminal of the outer loop proportional-integral controller 11, and the other input terminal is connected to the output terminal of the current sharing ring proportional-integral controller 13. The output terminal is connected to the input terminal of the limiter 15. One input terminal of the comparator 17 is connected to the output terminal of the adder 14, and the other input terminal is connected to the output terminal of the carrier phase-shift modulator 16. The output is a PWM wave.

[0147] The 24-pulse rectifier circuit is composed of four 6-pulse rectifier circuits connected in parallel. Each 6-pulse rectifier circuit includes a first diode, a second diode, a third diode, a fourth diode, a fifth diode, and a sixth diode. The anode of the first diode is connected to a three-phase voltage source A, the anode of the third diode is connected to a three-phase voltage source B, and the anode of the fifth diode is connected to a three-phase voltage source C. The cathodes of the second diode, the fourth diode, and the sixth diode are connected to a three-phase voltage source B and C, respectively. The first, third, and fifth diodes, when connected together, serve as the positive output terminal of the 6-pulse rectifier circuit, and the anodes of the second, fourth, and sixth diodes, when connected together, serve as the negative output terminal of the 6-pulse rectifier circuit.

[0148] The submodule adopts a Buck topology, which includes an IGBT, a diode, and an inductor L. The emitter of the IGBT is connected to the cathode of the diode and the anode of the inductor L. The collector of the IGBT serves as the positive input terminal of the Buck circuit, the anode of the diode serves as the negative input terminal of the Buck circuit, the output of the inductor L serves as the positive output terminal of the Buck circuit, and the anode of the diode serves as the negative output terminal of the Buck circuit.

[0149] According to one aspect of this application, a transfer function from the effective value of the hydrogen production power supply input voltage to the rectified average voltage is obtained, expressed by the following formula:

[0150] Where V d V is the rectified voltage. in U1 is the effective value of the 24-pulse input voltage, and U2 is the effective value of the phase voltage on the secondary side of the transformer.

[0151] According to one aspect of this application, the transfer function from the average inductor current of the interleaved parallel Buck to the duty cycle reference value can be obtained, expressed by the following formula:

[0152] Among them, i Lavg d is the average inductor current, d is the duty cycle reference value, C is the output capacitor, L is the inductance of each phase, R is the load resistance, and N is the number of phases in parallel.

[0153] The transfer function from the output voltage to the average inductor current can be obtained, expressed by the following formula:

[0154] The transfer function from the difference between the k-th phase inductor current and the average inductor current to the compensation duty cycle can be obtained, expressed by the following formula:

[0155] Among them, i Lk Let p be the inductor current of a certain phase, and p be the compensation duty cycle.

[0156] The expression for the terminal voltage of the electrolytic cell can be obtained, represented by the following formula:

[0157] Among them, V el V is the terminal voltage of the electrolytic cell. rev For a reversible potential, V ohm For overvoltage, V act Activation voltage, V diff This represents the concentration polarization potential.

[0158] According to one aspect of this application, expressions for reversible potential, overvoltage, and activation voltage can be obtained, expressed by the following formulas:

[0159]

[0160] Where V0 rev, T is the reversible voltage at standard pressure (1 bar) as a function of temperature, H is the universal gas constant (8.314 J / mol / K), and P v，KOH It is the vapor pressure of the electrolyte, α H α is the water activity in the electrolyte, r is a parameter related to the temperature of the electrolyzer, A is the cross-sectional area of ​​the region through which the current flows in the electrolyzer, I0 is the equilibrium current in which the oxidation reaction rate equals the reduction reaction rate, and α is the charge transfer coefficient.

[0161] According to one aspect of this application, an expression for the electrolytic cell current can be obtained, represented by the following formula:

[0162]

[0163] According to one aspect of this application, a complex domain expression for the current flowing through the electrolytic cell and the terminal voltage of the electrolytic cell can be obtained, expressed by the following formula:

[0164]

[0165] like Figure 3 As shown, the generation of the control signal includes the following steps:

[0166] Step 1: Acquire the system frequency value to obtain the system frequency deviation, and the inductor current I of the 20 sub-modules. L1 I L2 I L2 …I L20 ;

[0167] Step 2: Calculate the average value i of the submodule inductor current based on the submodule inductor current collected in Step 1. Lavg ;

[0168] Step 3: Obtain the inductor current command value based on the system frequency deviation and compare it with the average inductor current of the submodule obtained in Step 2. Lavg The difference is calculated by the outer loop subtractor 10, and the error signal is input to the outer loop proportional-integral controller 11;

[0169] Step 4: The average value of the submodule inductor current obtained in Step 2 is compared with the submodule inductor current obtained in Step 1. The difference is then input into the current sharing loop subtractor 12, and the error signal is input into the current sharing loop proportional-integral controller 13.

[0170] Step 5: Use the output of the outer loop proportional-integral controller 11 obtained in step 3 and the output of the current sharing loop proportional-integral controller 13 obtained in step 4 as the inputs of the adder 14, respectively.

[0171] Step 6: Use the output of adder 14 obtained in step 5 as the input of limiter 15;

[0172] Step 7: Use the output of the limiter 15 obtained in step 6 and the output of the carrier phase-shift modulator 16 as the input of the comparator 17, and the output is a PWM wave.

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

Claims

1. A method for primary frequency regulation control of a power system based on modular large-scale hydrogen production power sources, characterized in that: Includes the following steps: Step S1: Obtain the original PMU data stream P(t), electrolytic cell data stream E(t), and original submodule current data stream I(t); perform time synchronization, anomaly detection and repair, and noise reduction processing, and output the noise-reduced wavelet coefficients W' and feature vector F; Step S2: Input the denoised wavelet coefficients W' and eigenvector F into the unscented Kalman filter, and perform empirical mode decomposition and Hilbert transform to obtain the frequency disturbance feature matrix F_dist; calculate the power grid inertia estimate M_est based on this. Step S3: Input the electrolyzer data stream E(t) into the improved Hammerstein-Wiener model to obtain the electrolyzer dynamic model structure HW; and perform parameter identification and optimization to obtain the final model parameter vector θ_final. Step S4: Input the system state vector x_est, frequency disturbance feature matrix F_dist, grid inertia estimate M_est, and model parameter vector θ_final into the multi-time-scale control framework to obtain the control framework structure MSC; construct and solve the stochastic model predictive control problem based on the scenario tree to obtain the optimal control sequence u*.

2. The method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source as described in claim 1, characterized in that, It also includes step S5: Step S5: Input the control framework MSC, virtual synchronous generator parameters J_vsg and D_vsg, optimal control sequence u*, and long-term frequency regulation strategy L_opt into the distributed consensus algorithm module to obtain the global coordinated control instruction set C; and optimize it to output the optimized control parameter set P_opt.

3. The method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source as described in claim 2, characterized in that, Step S11: Collect basic data including voltage, current and phase angle from the synchronous phasor measurement unit, and obtain the raw PMU data stream P(t) after preprocessing; Voltage, current, temperature, pressure, and electrolyte concentration are collected from the electrolytic cell. After preliminary noise reduction, the original electrolytic cell data stream E(t) is obtained. Inductor current data is collected from N sub-modules, and the data is compressed to obtain the original sub-module current data stream I(t); Based on P(t), E(t), and I(t), the data quality index is calculated to obtain the data stream set D_q(t) with quality score; Step S12: Perform timestamp alignment, anomaly detection, data repair, and preliminary noise reduction on the data stream set D_q(t) to obtain the filtered data stream D_f(t); Step S13: Perform multi-scale decomposition, denoising, and feature extraction on the filtered data stream D_f(t) to obtain the feature vector F and the denoised wavelet coefficients.

4. The method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source as described in claim 2, characterized in that, Step S2 specifically involves: Step S21: Expand the dimension and perform nonlinear state prediction on the denoised wavelet coefficients W' and eigenvector F to obtain the system state vector x_est; Step S22: Extract the frequency component f_est from the system state vector x_est; perform mode decomposition, extraction and fusion to obtain the frequency perturbation feature matrix F_dist; Step S23: Based on the frequency disturbance feature matrix F_dist and the system state vector x_est, calculate the power loss ΔP and the system inertia and correct them to obtain the grid inertia estimate M_est.

5. The method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source as described in claim 2, characterized in that, Step S31: The electrolytic cell data stream E(t) is segmented, mapped by a polynomial kernel function, and identified by a system. Then, it is combined in series to obtain the dynamic model structure HW of the electrolytic cell. Step S32: Perform state estimation on the model structure HW and the electrolytic cell data stream E(t) to obtain the model parameter vector θ_est; Step S33: Optimize the parameter vector θ_est using particle filtering to obtain the final model parameter vector θ_final; Step S34: Evaluate the model using the final parameter vector θ_final and the electrolytic cell validation dataset E_val(t), and output the final optimized model parameter vector θ_opt.

6. The method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source as described in claim 2, characterized in that, Step S41: Input the system state vector x_est, frequency disturbance feature matrix F_dist, and grid inertia estimate M_est into the hierarchical control structure design module to obtain the multi-time-scale control framework structure MSC; including: the response time of the fast response layer controller is 0.1-1s; the response time of the medium-speed regulation layer controller is 1-10s; and the response time of the slow optimization layer controller is 10-60s. Step S42: Read the historical VSG parameter dataset D_hist from the pre-stored database; construct the grid inertia estimate M_est and the frequency disturbance feature matrix F_dist as the input feature vector x to obtain the virtual rotational inertia J_vsg and the virtual damping coefficient D_vsg; Step S43: Input the control frame MSC, VSG parameters [J_vsg, D_vsg] and model parameter vector θ_final into the dynamic scene tree generation module to obtain the optimal control sequence u*; Step S44: Input the optimal control sequence u into the adaptive learning module to obtain the final control input u_final; Step S45: Read historical control effect data H_ctrl from the pre-stored database; input u_final into the long-term optimization strategy generation module to obtain the final long-term frequency modulation strategy L_opt.

7. The method for participating in primary frequency regulation control of a power system based on a modular large-scale hydrogen production power source as described in claim 2, characterized in that, Step S5 specifically involves: Step S51: Input the control strategies MSC, J_vsg, D_vsg, u*, L_opt, and the submodule current data I(t) into the distributed coordination control module; obtain the local control law u_i; Step S52: Input the local control law u_i into the event triggering mechanism design module to generate the distributed model predictive control problem P_i, and apply the alternating direction multiplier method to solve for the local optimal control sequence u_i*; Step S53: Input the local optimal control sequence {u_i*} into the adaptive multi-scale evaluation module to obtain the Pareto optimal solution set P*; Step S54: Input the Pareto optimal solution set P* into the adaptive weight Kriging surrogate model to obtain the final Kriging surrogate model M_Kriging; Step S55: Input the surrogate model M_Kriging into the multi-objective Bayesian optimization module; Output the optimized set of control parameters P_opt.

Images