Primary frequency modulation control method based on participation of modular large-scale hydrogen production power supply in power system

Abstract

The invention provides a method for participating in primary frequency modulation control of a power system based on a modular large-scale hydrogen production power supply, and the method comprises the steps: obtaining an original data flow from a high-precision synchronous phasor measurement unit, an electrolytic cell and a submodule, and carrying out the multi-scale decomposition and noise reduction processing after time synchronization and anomaly detection restoration; performing state estimation by using unscented Kalman filtering, extracting frequency disturbance characteristics in combination with empirical mode decomposition and Hilbert transform, and estimating the inertia of the power grid; an electrolytic cell dynamic model is established, and parameter identification and optimization are carried out through a recursive least square method; a control strategy is designed based on a multi-time scale control framework, virtual synchronous generator parameters are optimized through online Gaussian process regression, and frequency modulation control is realized in combination with scene tree random model predictive control and minimum-maximum regret learning; and finally, carrying out coordination control through a distributed consistency algorithm and carrying out multi-target performance evaluation optimization. And the frequency modulation performance and the operation efficiency of the system are improved.

Description

Technical Field

[0001] The invention relates to power electronics technology, in particular to a method for controlling primary frequency regulation of a power system based on modular large-scale hydrogen production power supply. Background Art

[0002] With the large-scale grid connection of renewable energy, the frequency fluctuation problem of power system has become increasingly prominent. Traditional primary frequency regulation mainly relies on the rotational inertia and primary frequency regulation characteristics of conventional units. However, when a high proportion of renewable energy is connected, the system inertia is significantly reduced, and traditional frequency regulation methods are difficult to meet system requirements. As an important adjustable load, large-scale hydrogen production has the characteristics of fast response speed and large adjustment range. It is of great significance to introduce it into the primary frequency regulation of the power system. In particular, the modular hydrogen production system can effectively participate in the frequency regulation of the power grid through its flexible power regulation capability, while also realizing efficient storage and utilization of energy, providing important technical support for the construction of a new power system.

[0003] At present, the research on hydrogen production system participating in power frequency regulation mainly focuses on several aspects: first, the research on hydrogen production load characteristics based on simplified models, describing the dynamic characteristics of the electrolyzer through first-order or second-order inertia links; second, the use of a single time scale control strategy, mainly through PI control or model predictive control to achieve power tracking; third, a centralized control method for small-scale hydrogen production systems, realizing system scheduling through a unified control center. In terms of data processing, traditional filtering and state estimation methods such as Kalman filtering and extended Kalman filtering are mainly used; in terms of model identification, linear system identification methods such as least squares and subspace identification are mainly used; in terms of control strategy, optimization methods based on a single performance indicator are mainly used.

[0004] However, the existing technologies have the following problems when applied to large-scale modular hydrogen production systems: First, in terms of data processing, traditional time synchronization and anomaly detection methods are difficult to handle high sampling rate, multi-source heterogeneous data streams, especially when facing sudden abnormal data, the data repair quality is not high; second, the existing electrolyzer model is difficult to accurately describe the dynamic characteristics of the system under different working conditions, especially the nonlinear response characteristics when the power changes rapidly; third, the control strategy of a single time scale cannot take into account the rapid response and economy of the system at the same time, especially when dealing with the coordinated control of multiple subsystems, it is difficult to balance real-time and optimization; in addition, the traditional centralized control method has the problems of heavy communication burden and high computational complexity when facing a large number of distributed subsystems, and the scalability of the system is limited; finally, the existing performance evaluation methods often only focus on a single indicator, which is difficult to fully reflect the dynamic performance and economic performance of the system at multiple time scales. These technical problems seriously restrict the effect of large-scale hydrogen production systems participating in the primary frequency regulation of the power grid. Summary of the invention

[0005] The purpose of the invention is to provide a method for controlling the primary frequency regulation of a power system based on a modular large-scale hydrogen production power source, in the hope of solving one of the above-mentioned problems existing in the prior art.

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

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

[0008] Step S2, input the denoised wavelet coefficient W' and the eigenvector F into the 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 processing to obtain the frequency disturbance characteristic matrix F_dist; according to the frequency disturbance characteristic matrix F_dist and the system state vector x_est, apply the improved area method to calculate the power grid inertia estimation value M_est;

[0009] Step S3, inputting the electrolytic cell data stream E(t) into the improved Hammerstein-Wiener model to obtain the electrolytic cell dynamic model structure HW; applying the recursive least squares method to the model structure HW for parameter identification to obtain the model parameter vector θ_est; inputting the model parameter vector θ_est into the model verification module for parameter optimization to obtain the final model parameter vector θ_final;

[0010] Step S4, input the system state vector x_est, the frequency disturbance characteristic matrix F_dist, the grid inertia estimation value M_est, and the 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 modulation strategy L_opt.

[0011] Beneficial effect: Through multi-objective optimization and performance evaluation, a balance between system operation efficiency and frequency regulation performance is achieved, providing reliable control support for large-scale hydrogen production systems to participate in grid frequency regulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 It is a flow chart of the present invention.

[0013] Figure 2 This is the hydrogen production power supply topology diagram.

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

[0015] Figure 4 This is the block diagram of the frequency regulation control of the hydrogen production power supply. DETAILED DESCRIPTION

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

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

[0018] Step S2, input the denoised wavelet coefficient W' and the eigenvector F into the 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 processing to obtain the frequency disturbance characteristic matrix F_dist; according to the frequency disturbance characteristic matrix F_dist and the system state vector x_est, apply the improved area method to calculate the power grid inertia estimation value M_est;

[0019] Step S3, inputting the electrolytic cell data stream E(t) into the improved Hammerstein-Wiener model to obtain the electrolytic cell dynamic model structure HW; applying the recursive least squares method to the model structure HW for parameter identification to obtain the model parameter vector θ_est; inputting the model parameter vector θ_est into the model verification module for parameter optimization to obtain the final model parameter vector θ_final;

[0020] Step S4, input the system state vector x_est, the frequency disturbance characteristic matrix F_dist, the grid inertia estimation value M_est, and the 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 modulation strategy L_opt.

[0021] According to one aspect of the present application, the method further comprises step S5:

[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 consistency 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 realized, providing a reliable data basis for system control; secondly, through precise state estimation and model identification, the dynamic characteristics of the system and the state of the power grid are accurately grasped; thirdly, a multi-time scale control framework is adopted to realize the rapid response and economic optimization of the hydrogen production system to the power grid frequency; finally, through distributed coordinated control and multi-objective optimization, the efficient operation of large-scale systems is guaranteed. This systematic solution not only significantly improves the performance of the hydrogen production system in participating in the power grid frequency regulation, but also realizes the economy and reliability of the system operation, and provides 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 to achieve overall performance optimization, while ensuring the scalability and adaptability of the system.

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

[0025] Step S11, collecting basic data including voltage, current and phase angle from the synchronous phasor measurement unit, preprocessing through the adaptive sampling rate adjustment module, and obtaining the original PMU data stream P(t) with a sampling frequency of a predetermined frequency; collecting parameters such as voltage, current, temperature, pressure, electrolyte concentration, etc. from the electrolyzer, and performing preliminary noise reduction through a multivariable composite Kalman prefilter to obtain the original electrolyzer data stream E(t) with a sampling frequency of 10Hz; collecting inductor current data from N submodules, compressing the data through a distributed data compression sampling module, and obtaining the original submodule current data stream I(t) with a sampling frequency of 10kHz; inputting P(t), E(t) and I(t) into the data quality assessment module, calculating the data quality index in combination with the pre-stored historical data, and obtaining the data stream set D_q(t) with quality scores; N is a predetermined natural number;

[0026] Step S12, input the data stream set D_q(t) into the precise time protocol module, use the pre-stored time synchronization parameter set T to perform timestamp alignment, and obtain the time-aligned data stream D(t); input D(t) into the tensor decomposition anomaly detection module, obtain the low-rank tensor L and the 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 data set R; merge R with the normal data in D(t), and obtain the repaired data stream D'(t) after local manifold alignment; apply adaptive multi-scale filtering to D'(t) for preliminary denoising 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 IMF_k; apply biorthogonal wavelet transform to each intrinsic mode function IMF_k for multi-scale decomposition to obtain approximate coefficients a_J and detail coefficients d_j; input all coefficients into the adaptive threshold calculation module, generate the optimal threshold λ_j according to the standard deviation of each layer of noise, and perform soft threshold processing on the coefficients to obtain the denoised wavelet coefficients W'; input W' into the multi-dimensional feature extraction module, calculate the statistical features, entropy features and time-frequency features, and obtain the feature vector F; wherein the statistical features include mean, variance, skewness and kurtosis; the entropy features include sample entropy, fuzzy entropy and permutation entropy; the time-frequency features include Hilbert spectrum and wavelet packet energy.

[0028] By adopting high-precision synchronous phasor measurement units to obtain PMU data streams, electrolyzer data streams and sub-module current data streams, and performing time synchronization and data repair, unified processing of multi-source heterogeneous data is achieved. The specific innovations are: first, the adaptive sampling rate adjustment module is used to pre-process the PMU data to ensure the consistency of the sampling frequency; second, the multivariable composite Kalman pre-filter is used to perform preliminary noise reduction on the electrolyzer data to improve the data quality; third, the distributed data compression sampling module is used to compress the sub-module current data to reduce the data transmission burden. Time synchronization is performed through the precise time protocol module to achieve microsecond data alignment accuracy. Tensor decomposition anomaly detection and deep autoencoder network are used for data repair, which significantly improves data reliability. Finally, multi-scale decomposition and noise reduction are performed through variational mode decomposition and biorthogonal wavelet transform, and rich statistical features, entropy features and time-frequency features are extracted. This multi-level data processing scheme not only ensures the time consistency and reliability of the data, but also realizes effective noise reduction and feature extraction of the data, 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 data transmission and processing burden while maintaining the integrity of key characteristic information.

[0029] According to one aspect of the present application, step S2 is specifically:

[0030] Step S21, input the denoised wavelet coefficients W' and the eigenvector F into the adaptive state dimension expansion module, construct the dimension expansion state vector x, including the node voltage amplitude, phase angle and system frequency; apply the adaptive weight iterative traceless transformation to the dimension expansion state vector x, generate 2n+1 sigma point sets Xi and corresponding weights Wi; input Xi into the nonlinear state prediction equation, and calculate the predicted state x(k+1|k) in combination with 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) in combination with the observation noise covariance matrix R; update the state estimate through the sigma point inverse transformation and weighted summation to obtain the system state vector x_est;

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

[0032] Step S23, input the frequency disturbance characteristic matrix F_dist and the system state vector x_est into the power balance calculation module, and calculate the power imbalance ΔP in combination with the pre-stored power generation power P_gen and load power P_load data; extract the frequency change rate data from the frequency disturbance characteristic matrix F_dist, apply the adaptive time window selection algorithm, and determine the optimal analysis window [t1, t2]; within the selected window, use the inertia estimator to calculate the system inertia to 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 operation status, and use the recursive Bayesian estimation method to correct it to obtain the power grid inertia estimate M_est.

[0033] By combining unscented Kalman filtering, variational mode decomposition and inertia estimation, accurate estimation of power grid state and extraction of frequency disturbance features are achieved. Specifically, firstly, the optimal sigma point set and weight are generated through adaptive state dimension expansion and unscented transformation, which improves the accuracy of state estimation; secondly, the variational mode decomposition method is adopted to achieve accurate decomposition of frequency components by iteratively optimizing the center frequency and bandwidth constraints; thirdly, the multi-dimensional features of frequency disturbance are extracted through multi-scale permutation entropy analysis and time-varying fuzzy information entropy; finally, the inertia estimator and dynamic correction module are used to achieve accurate estimation of 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 characteristics of the system, and provide 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, so that the hydrogen production system can better respond to the frequency regulation needs of the power grid.

[0034] According to one aspect of the present application, step S3 is specifically:

[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 a segmented data set E_seg; apply polynomial kernel function mapping to each segmented data in the segmented data set 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 an augmented state vector x_aug, including the system state and the parameters to be identified; generate multiple groups 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 estimate from x_update to obtain the model parameter vector θ_est;

[0037] Step S33, input the parameter vector θ_est into the adaptive moment estimation module, construct a sliding time window W(k); calculate the high-order statistical moments M_1 to M_4 of the data in the window; adaptively adjust the forgetting factor λ(k) according to the changing trend of the moment; update the moment estimate M_update, and calculate the parameter correction Δθ; input the corrected parameters into the multi-model adaptive estimator, and 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 weight w_i(k); use particle filtering to optimize the parameters and obtain the final model parameter vector θ_final;

[0038] Step S34, input the final parameter vector θ_final and the electrolytic cell verification data set E_val(t) into the model evaluation module; construct a parameter constraint set θ, and design a projection-based constrained optimization problem; use a quadratic programming solver to obtain a parameter estimate θ_constrained that satisfies the constraints; input θ_constrained into the online model verifier to calculate the predicted output y_pred(t); compare y_pred(t) with the actual output y(t), and evaluate the model performance indicator J_perf; according to the performance evaluation results, use an iterative optimization algorithm to fine-tune the parameters and output the final optimized model parameter vector θ_opt.

[0039] The Hammerstein-Wiener model and multiple recursive least squares method are used to achieve accurate modeling of the dynamic characteristics of the electrolyzer. Specific innovations include: first, the optimal data segmentation is performed through the dynamic programming algorithm, and the composite kernel function is used for nonlinear static mapping, which improves the nonlinear expression ability of the model; second, the generalized matrix fractional decomposition and subspace identification algorithm are used to achieve accurate identification of linear dynamic characteristics; third, the real-time update of model parameters is achieved through the double recursive least squares method and adaptive moment estimation; finally, the robustness of parameter estimation is improved by using multi-model adaptive estimation and particle filter optimization. This comprehensive modeling scheme not only accurately describes the static and dynamic characteristics of the electrolyzer, but also realizes real-time update and optimization of parameters. Especially in large-scale hydrogen production systems, it can accurately capture the changes in dynamic characteristics under different operating conditions, providing accurate model support for the system to participate in grid frequency regulation. By considering multi-model fusion and parameter constraints, the adaptability and reliability of the model are significantly improved, so that the hydrogen production system can better adapt to the dynamic needs of grid frequency regulation.

[0040] According to one aspect of the present application, step S4 is specifically:

[0041] Step S41, input the system state vector x_est, the frequency disturbance characteristic matrix F_dist and the power grid inertia estimation value M_est into the hierarchical control structure design module; construct a fast response layer controller, output a fast control instruction u_f; design a medium-speed regulation layer controller, output a medium-speed control instruction u_m; configure a slow optimization layer controller, output a slow control instruction u_s; input the three-layer control instructions into the coordinated control module, apply the time-sharing and hierarchical optimization algorithm, and obtain a 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 optimization layer controller is 10-60s;

[0042] Step S42, read the historical VSG parameter data set D_hist from the pre-stored database; construct the grid inertia estimation value M_est and the frequency disturbance characteristic matrix F_dist as the input feature vector x; apply the deep kernel network to the feature vector x for feature transformation to obtain the transformation 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 by the marginal likelihood maximization criterion; calculate the covariance matrix K_opt and the characteristic decomposition result {U, Λ}; use sparse Gaussian process regression to perform parameter prediction, select the induced point set Z and calculate the cross-covariance matrix K_xz; update the VSG parameters through the adaptive particle swarm optimization algorithm to obtain the virtual moment of inertia J_vsg and the virtual damping coefficient D_vsg.

[0043] Step S43, input the control framework MSC, VSG parameters [J_vsg, D_vsg] and model parameter vector θ_final into the dynamic scene tree generation module; use the kernel density estimation method to construct the frequency disturbance probability distribution p(w); apply the stratified sampling strategy to generate the scene tree S_tree; execute the 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 conditions g; apply the distributed alternating direction multiplier method to solve the optimization problem and obtain the optimal control sequence u*.

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

[0045] Step S45, read the historical control effect data H_ctrl from the pre-storage database; input u_final output from 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 use the online gradient descent method to update the decision variables; adopt the adaptive learning rate adjustment strategy η(t); input the optimization result into the strategy evaluation module, and calculate the long-term performance indicator J_long; dynamically adjust the strategy based on J_long, and output the final long-term frequency modulation strategy L_opt.

[0046] The coordinated control of the hydrogen production system participating in the grid frequency regulation is realized by adopting a multi-time scale control framework and online Gaussian process regression. The specific innovations are: first, a hierarchical control structure including a fast response layer, a medium-speed regulation layer and a slow optimization layer is constructed to realize coordinated control at different time scales; second, the adaptive optimization of the parameters of the virtual synchronous generator is realized through a deep kernel network and online Gaussian process regression; third, the robustness of the control strategy is improved by the random model predictive control based on the scenario tree; finally, the optimization of the long-term frequency regulation strategy is realized through the minimax regret learning. This multi-time scale coordinated control scheme not only realizes the rapid response of the hydrogen production system to the grid frequency, but also ensures the economy and reliability of the system operation. Especially when dealing with large-scale hydrogen production systems, it can effectively coordinate the operation of multiple subsystems to optimize the overall performance of the system. By considering the control requirements of multiple time scales, the frequency regulation performance and operation efficiency of the system are significantly improved.

[0047] According to one aspect of the present application, step S5 is specifically:

[0048] Step S51, input the control strategy MSC, J_vsg, D_vsg, u*, L_opt, and sub-module current data I(t) into the distributed coordinated 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, including 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 weight, 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 trigger 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 chaos expansion method to process the parameter uncertainty, 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 use the alternating direction multiplier method to solve and obtain 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 Ji 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 frontier generation module; use the NSGA-III algorithm to initialize the population P(0), and design the reference point-based selection operation; 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 weight Kriging proxy model; use Latin hypercube sampling to generate the initial sample point set X; evaluate the objective function value Y of the sample point; construct a Kriging model to predict the function value and prediction variance of the unknown point; design the expected improvement criterion EI(x); use the adaptive weight strategy to dynamically adjust the EI weight w_EI; solve the weighted EI maximization problem to obtain a new sampling point x_new; update the sample set and retrain the model to obtain the final Kriging proxy model M_Kriging;

[0052] Step S55, input the proxy model M_Kriging into the multi-objective Bayesian optimization module; design the acquisition function α(x), and adopt the Thompson sampling strategy to sample from the posterior distribution; solve the sub-optimization 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 multi-objective evaluation results, Pareto optimal solution, sensitivity analysis and optimization suggestions, and output the optimized control parameter set P_opt.

[0053] Through distributed consensus algorithms and multi-objective performance evaluation, coordinated control and performance optimization of large-scale hydrogen production systems are achieved. Specific innovations include: first, the convergence and robustness of distributed control are improved by using graph analysis and weight optimization; second, the communication burden is reduced by event triggering mechanism and Zeno behavior exclusion; third, multi-objective coordinated optimization is achieved by using multi-scale performance evaluation and Pareto optimization; finally, the optimization efficiency is improved by using adaptive weighted Kriging agent model and multi-objective Bayesian optimization. This distributed coordinated control scheme not only realizes the efficient coordinated operation of large-scale hydrogen production systems, but also ensures the continuous optimization of system performance. Especially when dealing with a large number of distributed subsystems, it can effectively reduce the communication burden and improve the scalability of the system. Through multi-objective optimization and performance evaluation, a balance between system operation 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 the present application, step S11 is specifically:

[0055] Step S111, obtaining 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; inputting these data into the information entropy evaluator to calculate the information entropy H(t) of different time windows; using sliding variance analysis to calculate the signal fluctuation index σ(t); constructing an adaptive sampling rate adjustment function f_s(t) based on H(t) and σ(t); applying 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 V_s(t), I_s(t), θ_s(t) output by step S111 into the observation vector y(t); read the initial state estimate x(0) and the 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 untraceable transformation to generate the sigma point set χ_i and the weight W_i; perform state prediction and observation prediction for each sigma point; update the gain K(t) by cross-covariance calculation; obtain the pre-filtered PMU data stream P(t);

[0057] Step S113, collecting original 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); applying the wavelet dictionary learning algorithm to each parameter to construct a sparse representation basis D_i; using the orthogonal matching pursuit algorithm to calculate the compression coefficient α_i; designing an adaptive compression rate controller to dynamically adjust the compression rate r_i(t) according to the signal importance index; performing distributed compression sampling to obtain a compressed electrolytic cell data stream E(t);

[0058] Step S114, input the pre-filtered PMU data stream P(t) and the compressed electrolyzer data stream E(t) into the multidimensional data quality assessment module; collect the inductor current I_L(t) from 20 sub-modules to construct the 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, combined 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 a data stream set D_q(t) with quality score.

[0059] Firstly, an adaptive sampling rate adjustment function is constructed by using information entropy evaluation and sliding variance analysis to achieve dynamic optimization of the PMU data sampling rate, ensuring that higher sampling accuracy can be obtained when the system state changes rapidly, and the sampling rate is reduced to reduce the data burden when the system is stable; secondly, the electrolyzer data is preliminarily denoised by a multivariate composite Kalman prefilter, which 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, and efficient data compression is achieved through the wavelet dictionary learning algorithm and the orthogonal matching pursuit algorithm, which greatly reduces the data transmission burden while maintaining key information. This multi-level data acquisition and preprocessing scheme not only improves data quality, but also optimizes data transmission efficiency, which is particularly suitable for the multi-source heterogeneous data processing requirements in large-scale hydrogen production systems. In particular, the integrity, accuracy and consistency indicators are calculated through the multidimensional data quality assessment module, which realizes the quantitative evaluation of data quality and provides a reliable data basis for subsequent processing.

[0060] According to one aspect of the present application, step S12 is specifically:

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

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

[0063] Step S123, input the abnormal data set A output from step S122 into the multi-layer 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); use the attention mechanism to select similar patterns from the normal data tensor L as a reference; combine the reference pattern and the hidden layer representation, and generate repair data R = f_d(h;θ_d) through the decoder; apply residual learning to optimize the repair result to obtain the final repair data set R_opt;

[0064] Step S124, taking the normal data tensor L output by step S122 and the repaired data set R_opt output by step S123 as input; constructing a k-nearest neighbor graph G_k, and calculating the local geometric structure of each data point; designing an adaptive weight matrix W_l, and quantifying the local relationship between data points; using the manifold regularization algorithm to construct a global alignment objective function J = tr(X^TLX), where L is the graph Laplacian matrix; solving the alignment problem by alternating optimization to obtain the aligned data stream D'(t); applying adaptive multi-scale filtering to D'(t), and outputting the filtered data stream D_f(t).

[0065] Firstly, the transmission delay and clock deviation are accurately calculated by the two-way delay measurement algorithm, and the adaptive Kalman filter is used to realize real-time correction, ensuring the microsecond time synchronization accuracy; secondly, the Tucker decomposition method is used to reconstruct and decompose the synchronization data tensor, and the adaptive threshold algorithm is used to realize the accurate identification of abnormal data, which significantly improves the accuracy of anomaly detection; thirdly, the multi-layer autoencoder network combined with the attention mechanism is used to repair the abnormal data. This method achieves high-quality repair of abnormal data by selecting similar patterns from normal data as references; finally, by constructing a k-nearest neighbor graph and designing an adaptive weight matrix, the local manifold alignment of the data is realized, 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 the reliability and continuity of the data. Especially when processing high-sampling current data and low-sampling electrolyzer data, it can effectively coordinate the data synchronization problems of different time scales, providing high-quality data support for the coordinated control of the system.

[0066] According to one aspect of the present application, step S21 is specifically:

[0067] Step S211, input the denoised wavelet coefficients W' and the eigenvector F into the adaptive state dimension expansion module; construct the dimension expansion state vector x, including the node voltage amplitude, phase angle and system frequency; read the state dimension 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 dimension expansion strategy based on the correlation matrix C to obtain the optimized dimension expansion state vector x_opt; input the state vector x_opt into the covariance estimator to calculate the state initial covariance matrix P(0);

[0068] Step S212, input the expanded 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 of each sampling point; read the process noise statistical characteristics Q from the pre-stored parameter library, and construct a process noise covariance matrix Q_new considering state correlation; input the sigma point set {X_i} and the weight {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 transfer parameter θ_t from the pre-stored model library; apply the unscented transformation to each sampling point to calculate the one-step prediction state X_pred_i; design an abnormal prediction detector based on Mahalanobis distance to identify and correct abnormal prediction points; use the second-order moment matching method to calculate the predicted state mean x(k+1|k) and the 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 the ellipsoid constraint, 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; use the singular value decomposition method to calculate the cross-covariance matrix P_xy between the predicted state and the observation;

[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 a gradual memory characteristic to obtain the optimal filter gain K(k+1); use the sequential processing technology to update the state components one by one to obtain the state estimate x_est(k+1); perform a consistency check on the state estimation result and the prediction result to determine the estimation reliability; output the system state vector x_est.

[0072] Firstly, the state correlation matrix is ​​calculated by the Cramer decomposition method, and the dimension expansion strategy is dynamically adjusted based on the correlation matrix to improve the dimensional adaptability of the state estimation. Secondly, the spherical-radial sampling module is used to generate a uniformly distributed sigma point set, and the sampling effect is optimized by an adaptive weight allocation strategy based on state sensitivity. Thirdly, the unscented transformation is used for state prediction, and the prediction reliability is improved by combining the abnormal prediction detector based on the Mahalanobis distance. Finally, the gain matrix calculation method with the fading memory characteristic is used to achieve stable update of the state estimation. This innovative state estimation scheme not only improves the estimation accuracy, but also enhances the robustness of the algorithm. Especially when dealing with multivariable state estimation of large-scale hydrogen production systems, it can effectively handle the coupling relationship between state variables and provide accurate state information for system control. Through sequential processing technology and consistency verification, the reliability of the state estimation results is further guaranteed, so that the hydrogen production system can respond more accurately to the grid frequency regulation needs.

[0073] According to one aspect of the present application, step S22 is specifically:

[0074] Step S221, input the frequency component f_est in the system state vector x_est into the adaptive bandwidth optimization module; construct a variational mode decomposition objective function, including 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 decomposition level K according to the signal energy distribution; perform variational mode decomposition on f_est to obtain a series of modal functions {u_1,u_2,...,u_K} and corresponding center frequencies {ω_1,ω_2,...,ω_K}.

[0075] Step S222, input the modal function set {u_k} output by step S221 into the multi-scale analysis module; construct a coarse-grained time series {x_τ} of different scales for each mode u_k; design an adaptive embedding dimension selector, and 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 the curvature c_k; output the permutation entropy feature vector PE_features.

[0076] Step S223, input the modal function {u_k} of step S221 and the PE_features of step S222 into the time-frequency analysis module; apply the generalized S transform to each mode 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 by 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; use local response normalization to process feature data; extract multi-level features through residual connections and jump connections; design feature selection gating units to dynamically screen important features; apply dropout regularization to prevent overfitting; perform principal component analysis on the network output to reduce dimensionality 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 through the adaptive bandwidth optimization module to improve the accuracy of modal decomposition; secondly, the complexity characteristics of the frequency signal are extracted by multi-scale permutation entropy analysis, and the accuracy of feature extraction is ensured by the adaptive embedding dimension selector; thirdly, the generalized S transform is used for time-frequency analysis, and the accurate extraction of time-frequency features is achieved by combining the adaptive time-frequency resolution controller; finally, the effective fusion of multi-level features is achieved through the deep learning network and feature selection gating unit. 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 mode of the system during the frequency modulation process. Especially when dealing with large-scale hydrogen production systems participating in power grid frequency modulation, it can accurately analyze the frequency response capability and dynamic characteristics of the system, providing rich feature information for subsequent control optimization. Through multi-feature fusion and dimensionality reduction processing, the effectiveness and computational efficiency of feature expression are significantly improved.

[0079] According to one aspect of the present application, step S31 is specifically:

[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; use wavelet transform to perform multi-scale analysis on the data and identify data feature points; construct a segmentation boundary detector based on 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 a segmented data set E_seg;

[0081] Step S312, input the segmented data set 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; use kernel principal component analysis to perform feature extraction to obtain kernel feature mapping φ(u); design a sparse kernel mapping method based on L1 regularization to construct a 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; 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; use a piecewise linear approximator to approximate the nonlinear characteristics; design a smooth constraint based on the B-spline function 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.

[0084] Firstly, the segmented boundary detector and dynamic programming algorithm based on Kolmogorov-Smirnov test are used to achieve the optimal segmentation of data and improve the local accuracy of modeling. Secondly, the accurate mapping of nonlinear static characteristics is achieved by constructing a composite kernel function based on Mercer's theorem and kernel principal component analysis. Thirdly, the dynamic characteristics of the system are accurately identified by using generalized matrix fraction decomposition and subspace identification algorithm. Finally, a complete nonlinear dynamic model is constructed by 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 adaptability of the model under different working conditions. Especially when dealing with large-scale hydrogen production systems, it can effectively capture the dynamic response characteristics of the system under different operating conditions, providing a reliable model support for the system to participate in grid frequency regulation. The combination of segmented modeling and nonlinear mapping significantly improves the local accuracy and global adaptability of the model.

[0085] According to one aspect of the present application, step S32 is specifically:

[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 to include physical quantities such as voltage and current, and construct the parameter vector θ to include nonlinear mapping coefficients and dynamic parameters; combine x and θ to form an 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 transfer equation f(X_i) to each Sigma point to predict the state at the next moment; 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, taking the predicted state x_pred outputted by step S321 as input; constructing a sliding time window W(t), wherein the window length L is adaptively adjusted according to the signal change rate; calculating the first-order moment M_1 to the fourth-order moment M_4 of the data in the window; designing a forgetting factor update rule based on the moment change rate, and obtaining an adaptive forgetting factor λ(t); applying an exponentially weighted recursive algorithm to update the moment estimate, M_new = λ(t)M_old + (1-λ(t))M_current; calculating the parameter correction Δθ according to the high-order moment characteristics, and obtaining the corrected parameter estimate θ_corr;

[0088] Step S323, inputting the θ_corr outputted from step S322 into the multi-model library; constructing N candidate models {M_1, M_2, ..., M_N} according to historical data and expert knowledge; calculating the likelihood function p(y(k)|M_i) for each model; designing a weight update mechanism based on minimizing the prediction error; using the softmax function to calculate the model weight w_i(k); constructing an interactive multi-model estimation framework to handle the complementarity between models; fusing 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 each particle for prediction; calculate the particle likelihood based on the latest observation data; design an adaptive resampling strategy based on the number of effective particles; construct a moving resampling algorithm to avoid particle degradation; update the particle weights and normalize them; calculate the mean and variance of the posterior distribution; obtain the final parameter vector θ_final through iterative optimization.

[0090] First, an augmented state vector containing the physical state and the parameters to be identified is constructed, and the accuracy of parameter estimation is improved by generating the optimal Sigma point set; secondly, the nonlinear state prediction and observation prediction equations are used, combined with adaptive window length and forgetting factor, to achieve real-time update of parameters; thirdly, the robustness of parameter estimation is improved through high-order moment estimation and adaptive weight strategy; finally, parameter optimization under multi-model fusion is achieved by likelihood calculation and model weight update. 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 of large-scale hydrogen production systems, it can effectively deal with the time-varying characteristics and uncertainty of system parameters, and provide reliable parameter support for the precise control of the system. Through multi-model fusion and adaptive moment estimation, the robustness and real-time performance of parameter identification are significantly improved.

[0091] According to one aspect of the present application, step S33 is specifically:

[0092] Step S331, input parameter vector θ_est into adaptive moment estimation module; read window initial parameter W_init from pre-stored parameter library; construct dynamic window size determination criterion using generalized extreme value theory; calculate high-order statistical moments M_1 to M_4 of data in window based on local polynomial fitting; design window update mechanism based on minimum description length criterion, dynamically adjust sliding window W(k); calculate covariance matrix of 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 of the forgetting factor λ_0 from the pre-stored parameter library; construct a moment change detector using Kullback-Leibler divergence; calculate the change 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); use the Kalman filtering algorithm to update the moment estimate M_update;

[0094] Step S333, input the updated moment estimate M_update into the parameter correction module; read the parameter constraint condition Θ_con from the pre-stored parameter library; calculate the parameter gradient using the quasi-Newton method to obtain the correction value Δθ; 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; use the variational Bayes method to calculate the posterior probability p(M_i|y(k)) of each model; 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 parameters; determine the confidence interval of the parameters using the quantile regression method; and output the final model parameter vector θ_final.

[0097] First, the high-order statistical moments are calculated through a sliding time window, and the forgetting factor is dynamically adjusted based on the changing trend of the moment, which improves the real-time performance of parameter estimation. Secondly, multiple candidate models are constructed using a multi-model adaptive estimator, and the dynamic update of model weights is achieved through likelihood probability calculation. Thirdly, particle filtering is used for parameter optimization, and the particle degradation problem is avoided through a moving resampling algorithm. Finally, parameter constraints and online model verification are used to ensure the physical feasibility of parameter estimation. This comprehensive parameter optimization scheme not only ensures the accuracy of parameter estimation, but also realizes real-time update and constraint processing of parameters. Especially when dealing with parameter optimization of large-scale hydrogen production systems, it can effectively handle the uncertainty and time-varying characteristics of parameters, providing reliable parameter support for the stable operation of the system. Through multi-model fusion and constraint optimization, the reliability and adaptability of parameter estimation are significantly improved.

[0098] According to one aspect of the present application, step S41 is specifically:

[0099] Step S411, input the system state vector x_est, the frequency disturbance characteristic matrix F_dist and the grid inertia estimation value M_est into the hierarchical control structure design module; read the control level parameter set L_param from the pre-stored parameter library; design the inter-level 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 a 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 parameter M_param from the pre-stored parameter library; design a robust controller using H∞ control theory; construct a controller parameter optimization problem based on linear matrix inequality; design a feedforward compensator based on generalized predictive control; solve the controller parameters using sequential quadratic programming; and construct a medium speed regulation layer controller structure C_m;

[0101] Step S413, input the medium-speed regulation layer controller C_m into the slow optimization layer design module; read the slow control parameter S_param from the pre-stored parameter library; construct the optimal control problem using the Pontryagin maximum principle; design the economic dispatch strategy based on the dynamic programming algorithm; design the constraint processing mechanism based on model prediction; apply the interior point method to solve the optimization problem; construct the slow optimization layer controller structure C_s;

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

[0103] Firstly, a hierarchical control structure including a fast response layer, a medium-speed regulation layer and a slow optimization layer was constructed to achieve control coordination at different time scales; secondly, a medium-speed regulation layer controller was designed through H∞ control theory and generalized predictive control to improve the robustness of the system; thirdly, the economic dispatch strategy of the slow layer was optimized using the Pontryagin maximum principle and dynamic programming algorithm; finally, the coordination of multi-layer control was achieved through time-sharing and hierarchical optimization. This multi-time scale control scheme not only meets the rapid response requirements of the hydrogen production system participating in the grid frequency regulation, but also ensures the economy of the system operation. Especially when dealing with large-scale hydrogen production systems, it can effectively coordinate the control objectives of different time scales and achieve the overall optimization of system performance. Through hierarchical control and coordinated optimization, the control performance and operation efficiency of the system are significantly improved.

[0104] According to one aspect of the present application, step S42 is specifically:

[0105] Step S421, read the historical VSG parameter data set 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 an initial kernel feature map φ_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 a representation x_trans that highlights the key features; use the back propagation algorithm to optimize the network parameters 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 induction points Z={z_1,...,z_m}; calculate the cross-covariance matrix K_xz and K_zz; construct a sparse approximate prediction equation to obtain the VSG parameter prediction value [J_pred,D_pred] and prediction variance [v_J,v_D].

[0107] Step S423, take the parameter prediction value [J_pred, D_pred] output by 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, c2 according to the performance gradient; apply the speed update formula to avoid local optimality; iteratively optimize to obtain the parameter optimization value [J_opt, D_opt].

[0108] Step S424, input the optimization parameters [J_opt, D_opt] output from 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 the structural singular value analysis method 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 to generate compensation [Δ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, and the effectiveness of feature expression is improved by combining the attention mechanism; secondly, the optimal kernel function is selected by using the marginal likelihood maximization criterion, and parameter prediction is achieved through sparse Gaussian process regression; thirdly, the adaptive particle swarm optimization algorithm is used to update the VSG parameters, which improves the efficiency of parameter optimization; finally, the stability of the parameters is ensured through robustness analysis. 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 the changes in the dynamic characteristics of the system and realize real-time optimization of VSG parameters. Through the combination of feature learning and parameter optimization, the frequency modulation performance of the system is significantly improved.

[0110] The step S51 is specifically as follows:

[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 the adjacency matrix A; calculate the network algebraic connectivity λ_2 as an evaluation index; design a topology optimization objective function J_topo based on graph analysis; use a discrete particle swarm algorithm to search for the optimal topology structure; consider the communication delay τ_ij and the packet loss rate p_loss to construct the communication quality constraint; apply the branch and bound method to solve the discrete optimization problem; obtain the optimized topology structure G_opt and the initial weight matrix W_init;

[0112] Step S512, take the initial weight matrix W_init outputted in step S511 as input; construct a local state vector x_i, including the inductor current I_i, the output voltage V_i and the 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 projected 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 from 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; use the interior point method to solve the LMI optimization problem to obtain the feedback gain matrix K_i; design the observer gain L_i based on event triggering; integrate state feedback and observer design to obtain a complete controller structure C_i. ;

[0114] Step S514, input the controller structure C_i output from step S513 into the distributed optimization module; construct the local cost function J_i and the coupling constraint g_i; design a distributed solution framework based on ADMM; introduce the augmented Lagrangian to process the coupling constraint; construct a prediction model for each subsystem to predict the state trajectory x_i(k+j|k); update local variables, global variables and Lagrange multipliers through iterative alternating optimization; 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 to improve the convergence of distributed control. Second, the communication weight is dynamically updated using an adaptive weight optimizer to enhance the robustness of control. Third, local control optimization is achieved through the design of a distributed H∞ controller. Finally, the communication burden is reduced based on an event trigger mechanism. This distributed coordinated control scheme not only achieves efficient control of large-scale systems, but also ensures the reliability of control. Especially when dealing with a large number of distributed subsystems, it can effectively reduce communication requirements and improve the scalability of the system. Through topology optimization and distributed control, the overall performance of the system is significantly improved.

[0116] According to one aspect of the present application, step S52 is specifically:

[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 a dynamic dead zone characteristic based on exponential decay; construct a trigger compensation mechanism considering 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 threshold dynamic adjustment strategy based on local minimum analysis; design a normalized trigger condition based on relative error; apply an extreme value optimization algorithm to solve the optimal threshold σ_i(t); construct a threshold smoothing mechanism based on an exponential sliding 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; apply the convex optimization method to solve the minimum trigger interval τ_min_new; 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 chaos expansion method; analyze the influence of parameter disturbance based on Wiener-Askey expansion; design a stability criterion based on the probability ellipsoid; apply the Monte Carlo method to evaluate the system robustness; 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 the dual decomposition theory; design a distributed coordination mechanism based on consistency constraints; apply the proximal point method to solve sub-problems; and output the local optimal control sequence u_i*.

[0122] Firstly, the triggering time of the control signal is optimized by designing the state-dependent trigger function and updating the adaptive threshold. Secondly, the Zeno behavior exclusion mechanism is used to ensure the minimum triggering interval and avoid frequent triggering. Thirdly, the parameter uncertainty is analyzed by the polynomial chaos expansion method to expand the stability region of the system. Finally, the distributed model predictive control is used to achieve local optimal control. This not only reduces the communication burden of the system, but also ensures the real-time and reliability of the control. The efficient use of communication resources is achieved through the dynamic adjustment of the adaptive trigger threshold. At the same time, the stable operation of the system under parameter uncertainty is ensured through robustness analysis. Especially when dealing with the distributed control of large-scale hydrogen production systems, it can effectively balance the control performance and communication overhead, providing reliable technical support for the large-scale deployment of the system. The combination of event triggering mechanism and distributed optimization significantly improves the control efficiency and scalability of the system, 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 improving the overall performance of the system.

[0123] According to one aspect of the present application, step S53 is specifically:

[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 component {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 hierarchical analysis method; 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 of each scale;

[0126] Step S533, input each scale weight w_i,k into the index fusion module; read the fusion parameter set F_param from the pre-stored parameter library; use Choquet integral to build a nonlinear index fusion model; design an interactive characteristic evaluation method based on fuzzy measure theory; design a dynamic fusion strategy based on coordination; apply a convex optimization algorithm to solve the fusion weight; and obtain a 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 sorting; design a genetic operation based on adaptive crossover rate; apply the non-dominated sorting method to update the population; output the Pareto optimal solution set P*.

[0128] Firstly, the wavelet packet transform algorithm is used for multi-scale decomposition, and the optimal number of decomposition layers is selected by the Shannon entropy criterion to achieve a refined analysis of performance indicators; secondly, the entropy weight method and the hierarchical analysis method are used to construct a comprehensive evaluation system of objective weights and subjective weights to ensure the scientific nature of weight distribution; thirdly, the Choquet integral is used to achieve effective fusion of nonlinear indicators, and the interaction between indicators is considered; finally, the NSGA-III algorithm is used to solve the multi-objective optimization problem and obtain 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. In particular, when evaluating the frequency modulation performance of large-scale hydrogen production systems, it can accurately capture the performance characteristics at different time scales 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 the present application, step S54 is specifically:

[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 the cluster analysis method to evaluate the representativeness of the sample; and 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 an anisotropic correlation function structure; apply the diagonalization algorithm to solve the correlation matrix; and 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 the exploration-utilization balance strategy based on the confidence limit ratio; design the local search mechanism based on the gradient information; apply the sequential quadratic programming method to solve 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 proxy model; output the final Kriging proxy model M_Kriging.

[0134] By constructing a response surface model and calculating robustness indicators, the reliability evaluation of the optimization results is achieved, ensuring the effectiveness of the optimization scheme in practical applications.

[0135] In this embodiment, the variational modal decomposition method introduces adaptive bandwidth optimization and center frequency dynamic adjustment mechanisms, so that the modal decomposition process can automatically adjust the decomposition parameters according to the signal characteristics. Specifically, an adaptive modal 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 modal decomposition method using fixed parameters.

[0136] The Hammerstein-Wiener model adds a dynamic programming optimal segmentation and composite kernel function mapping mechanism based on the traditional HW model. Specifically, the data is first segmented optimally using the dynamic programming algorithm, and then the composite kernel function is constructed by combining the radial basis kernel and the polynomial kernel, which improves the model's ability to describe the nonlinear characteristics of the electrolytic cell. At the same time, the B-spline function is used to implement smooth constraints to avoid sudden changes in the model output.

[0137] The double recursive least squares method, combined with adaptive moment estimation and sliding time window, realizes real-time parameter update. Specifically, it constructs an augmented state vector, which includes the system state and the parameters to be identified, dynamically adjusts the forgetting factor by calculating high-order statistical moments (first to fourth order moments), and uses the Metropolis-Hastings algorithm to generate parameter samples, which improves the accuracy of parameter estimation.

[0138] The multi-time scale control framework divides the control system into a fast response layer (0.1-1s), a medium-speed adjustment layer (1-10s) and a slow optimization layer (10-60s) to achieve coordinated control of different time scales. Specifically, the fast layer uses H∞ control to achieve fast response, the medium-speed layer uses generalized predictive control for regulation, and the slow layer performs economic optimization based on the Pontryagin maximum principle, achieving global performance optimization.

[0139] The adaptive weight Kriging surrogate model 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 to achieve optimal selection of sampling points.

[0140] Multi-layer deep kernel network adds multi-layer feature extraction and attention mechanism on the basis of traditional kernel network. Specifically, it designs multi-layer nonlinear transformation and residual connection structure, highlights key features through attention mechanism, and realizes more effective feature learning and parameter prediction.

[0141] The event-triggered distributed control mechanism, combined with the adaptive trigger threshold and Zeno behavior exclusion mechanism, optimizes the trigger strategy of the control signal. Specifically, it designs a state-dependent trigger function, dynamically adjusts the trigger condition through an adaptive threshold, and ensures the minimum trigger interval to avoid frequent communication problems.

[0142] The distributed multi-objective performance evaluation method integrates wavelet packet transform and Choquet integral to achieve comprehensive evaluation of multi-dimensional performance indicators. Specifically, the optimal number of decomposition levels is selected through the Shannon entropy criterion, the weights are determined using the entropy weight method and the hierarchical analysis method, and the interaction between indicators is considered through the Choquet integral to achieve a comprehensive performance evaluation.

[0143] The NSGA-III algorithm adds reference point-based selection operations and adaptive crossover rate design. Specifically, the next generation population is selected through the reference point association mechanism, and the crossover and mutation probabilities are dynamically adjusted, which improves the convergence and diversity of multi-objective optimization.

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

[0145] In another embodiment of the present application, the structure of the hardware part is:

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

[0147] The 24-pulse rectifier circuit is composed of four 6-pulse rectifier circuits connected in parallel. The 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 the three-phase voltage source A, the anode of the third diode is connected to the three-phase voltage source B, the anode of the fifth diode is connected to the three-phase voltage source C, the cathode of the second diode is connected to the three-phase voltage source A, the cathode of the fourth diode is connected to the three-phase voltage source B, and the cathode of the sixth diode is connected to the three-phase voltage source C. The cathodes of the first diode, the third diode and the fifth diode are connected to serve as the positive output terminal of the 6-pulse rectifier circuit, and the anodes of the second diode, the fourth diode and the sixth diode are connected to serve as the negative output terminal of the 6-pulse rectifier circuit.

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

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

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

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

[0152] Among them, i Lavg is the average value of the 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 value of the inductor current can be obtained, which is expressed by the following formula:

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

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

[0156] The voltage expression of the electrolytic cell terminal can be obtained, which is expressed by the following formula:

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

[0158] According to one aspect of the present application, the reversible potential, overvoltage, and activation voltage expressions can be obtained, which are 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), P v，KOH is the vapor pressure of the electrolyte, α H is the water activity in the electrolyte, r is a parameter related to the cell temperature, A is the cross-sectional area of ​​the region through which current flows in the cell, I0 is the equilibrium current at which the oxidation reaction rate is equal to the reduction reaction rate, and α is the charge transfer coefficient.

[0161] According to one aspect of the present application, an electrolytic cell current expression can be obtained, which is expressed by the following formula:

[0162]

[0163] According to one aspect of the present application, a complex domain expression of the current flowing through the electrolytic cell and the voltage at the electrolytic cell terminal can be obtained, which is 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: Collect the system frequency value to obtain the system frequency deviation and the inductor current I of the 20 submodules L1 ,I L2 ,I L2 …I L20 ;

[0167] Step 2: Calculate the average value of the submodule inductor current i according to 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 i of the submodule obtained in step 2 Lavg The error signal is input into the outer loop subtractor 10 for difference, and then input into the outer loop proportional integral controller 11;

[0169] Step 4: input the average value of the submodule inductor current obtained in step 2 and the submodule inductor current obtained in step 1 into the current sharing loop subtractor 12 to make a difference, and input the error signal into the current sharing loop proportional integral controller 13;

[0170] Step 5: 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 are respectively used as inputs of the adder 14;

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

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

[0173] The preferred embodiments of the present invention are described in detail above; however, the present invention is not limited to the specific details in the above embodiments. Within 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 belong to the protection scope of the present invention.

Claims

1. A method for controlling the primary frequency regulation of a power system based on a modular large-scale hydrogen production power source, characterized in that: The steps include: Step S1, obtaining the original PMU data stream P(t) from the high-precision synchronized phasor measurement unit, obtaining the original electrolyzer data stream E(t) from the electrolyzer, and obtaining the original submodule current data stream I(t) from the submodule; inputting P(t), E(t) and I(t) into the precise time protocol module for time synchronization, obtaining the time-aligned data stream D(t) and inputting it into the robust principal component analysis module for anomaly detection and repair, obtaining the repaired data stream D'(t); Apply biorthogonal wavelet transform to the repaired data stream D'(t) for multi-scale decomposition and denoising, and output the denoised wavelet coefficients W' and eigenvector F; Step S2, input the denoised wavelet coefficient W' and the eigenvector F into the 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 processing to obtain the frequency disturbance characteristic matrix F_dist; according to the frequency disturbance characteristic matrix F_dist and the system state vector x_est, apply the improved area method to calculate the power grid inertia estimation value M_est; Step S3, inputting the electrolytic cell data stream E(t) into the improved Hammerstein-Wiener model to obtain the electrolytic cell dynamic model structure HW; applying the recursive least squares method to the model structure HW for parameter identification to obtain the model parameter vector θ_est; inputting the model parameter vector θ_est into the model verification module for parameter optimization to obtain the final model parameter vector θ_final; Step S4, input the system state vector x_est, the frequency disturbance characteristic matrix F_dist, the grid inertia estimation value M_est, and the 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 modulation strategy L_opt.

2. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 1, characterized in that: The step S5 is also included: 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 consistency 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.

3. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: Step S11, collecting basic data including voltage, current and phase angle from the synchronous phasor measurement unit, preprocessing through the adaptive sampling rate adjustment module, and obtaining the original PMU data stream P(t) with a sampling frequency of a predetermined frequency; collecting parameters such as voltage, current, temperature, pressure, electrolyte concentration, etc. from the electrolyzer, and performing preliminary noise reduction through a multivariable composite Kalman prefilter to obtain the original electrolyzer data stream E(t) with a sampling frequency of 10Hz; collecting inductor current data from N submodules, compressing the data through a distributed data compression sampling module, and obtaining the original submodule current data stream I(t) with a sampling frequency of 10kHz; inputting P(t), E(t) and I(t) into the data quality assessment module, calculating the data quality index in combination with the pre-stored historical data, and obtaining the data stream set D_q(t) with quality scores; N is a predetermined natural number; Step S12, input the data stream set D_q(t) into the precise time protocol module, use the pre-stored time synchronization parameter set T to perform timestamp alignment, and obtain the time-aligned data stream D(t); input D(t) into the tensor decomposition anomaly detection module, obtain the low-rank tensor L and the 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 data set R; merge R with the normal data in D(t), and obtain the repaired data stream D'(t) after local manifold alignment; apply adaptive multi-scale filtering to D'(t) for preliminary denoising to obtain the filtered data stream D_f(t); Step S13, input the filtered data stream D_f(t) into the variational mode decomposition module to obtain multiple intrinsic mode functions IMF_k; apply biorthogonal wavelet transform to each intrinsic mode function IMF_k for multi-scale decomposition to obtain approximate coefficients a_J and detail coefficients d_j; input all coefficients into the adaptive threshold calculation module, generate the optimal threshold λ_j according to the standard deviation of each layer of noise, and perform soft threshold processing on the coefficients to obtain the denoised wavelet coefficients W'; input W' into the multi-dimensional feature extraction module, calculate the statistical features, entropy features and time-frequency features, and obtain the feature vector F; wherein the statistical features include mean, variance, skewness and kurtosis; the entropy features include sample entropy, fuzzy entropy and permutation entropy; the time-frequency features include Hilbert spectrum and wavelet packet energy.

4. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: The step S2 is specifically as follows: Step S21, input the denoised wavelet coefficients W' and the eigenvector F into the adaptive state dimension expansion module, construct the dimension expansion state vector x, including the node voltage amplitude, phase angle and system frequency; apply the adaptive weight iterative traceless transformation to the dimension expansion state vector x, generate 2n+1 sigma point sets Xi and corresponding weights Wi; input Xi into the nonlinear state prediction equation, and calculate the predicted state x(k+1|k) in combination with 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) in combination with the observation noise covariance matrix R; update the state estimate through the sigma point inverse transformation and weighted summation to obtain the system state vector x_est; Step S22, extracting the frequency component f_est from the system state vector x_est; inputting f_est into the variational mode decomposition module, and obtaining a series of intrinsic mode functions IMF_k by iteratively optimizing the center frequency ω_k and the bandwidth constraint τ_k; applying adaptive multi-scale permutation entropy analysis to each intrinsic mode function IMF_k, calculating the permutation entropy PE_k(τ) at different scales, and extracting entropy value features; inputting IMF_k into the time-varying fuzzy information entropy module, constructing the time-varying fuzzy membership function μ_k(t), and calculating the fuzzy information entropy FIE_k(t); applying Hilbert transform to each intrinsic mode function IMF_k, calculating the instantaneous frequency IF_k(t) and the instantaneous amplitude IA_k(t); integrating all features using a multi-layer feature fusion network to obtain a frequency disturbance feature matrix F_dist; Step S23, input the frequency disturbance characteristic matrix F_dist and the system state vector x_est into the power balance calculation module, and calculate the power imbalance value ΔP in combination with the pre-stored power generation power P_gen and load power P_load data; extract the frequency change rate data from the frequency disturbance characteristic matrix F_dist, apply the adaptive time window selection algorithm, and determine the optimal analysis window [t1, t2]; in the selected window, use the inertia estimator to calculate the system inertia to obtain the preliminary inertia estimation value M_raw; The preliminary inertia estimation value M_raw is input into the dynamic correction module, and combined with the historical inertia estimation results and the system operation status, the recursive Bayesian estimation method is used to correct it to obtain the grid inertia estimation value M_est.

5. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: 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 a segmented data set E_seg; apply polynomial kernel function mapping to each segmented data in the segmented data set 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; Step S32, input the model structure HW and the electrolytic cell data stream E(t) into the dual recursive least squares module; construct an augmented state vector x_aug, including the system state and the parameters to be identified; generate multiple groups 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 estimate from x_update to obtain the model parameter vector θ_est; Step S33, input the parameter vector θ_est into the adaptive moment estimation module to construct a sliding time window W(k); Calculate the high-order statistical moments M_1 to M_4 of the data in the window; According to the changing trend of the moment, adaptively adjust the forgetting factor λ(k); update the moment estimate M_update and calculate the parameter correction Δθ; input the corrected parameters into the multi-model adaptive estimator and 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 weight w_i(k); use particle filtering to optimize the parameters and obtain the final model parameter vector θ_final; Step S34, input the final parameter vector θ_final and the electrolytic cell verification data set E_val(t) into the model evaluation module; construct a parameter constraint set θ, and design a projection-based constrained optimization problem; use a quadratic programming solver to obtain a parameter estimate θ_constrained that satisfies the constraints; input θ_constrained into the online model verifier to calculate the predicted output y_pred(t); compare y_pred(t) with the actual output y(t), and evaluate the model performance indicator J_perf; according to the performance evaluation results, use an iterative optimization algorithm to fine-tune the parameters and output the final optimized model parameter vector θ_opt.

6. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: Step S41, input the system state vector x_est, the frequency disturbance characteristic matrix F_dist and the power grid inertia estimation value M_est into the hierarchical control structure design module; construct a fast response layer controller, output a fast control instruction u_f; design a medium-speed regulation layer controller, output a medium-speed control instruction u_m; configure a slow optimization layer controller, output a slow control instruction u_s; input the three-layer control instructions into the coordinated control module, apply the time-sharing and hierarchical optimization algorithm, and obtain a 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 optimization layer controller is 10-60s; Step S42, read the historical VSG parameter data set D_hist from the pre-stored database; construct the grid inertia estimation value M_est and the frequency disturbance characteristic matrix F_dist as the input characteristic vector x; Apply deep kernel network to transform the feature vector x to obtain 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 by maximizing the marginal likelihood criterion; calculate the covariance matrix K_opt and the feature decomposition result {U, Λ}; use sparse Gaussian process regression to predict parameters, select the induced point set Z and calculate the cross-covariance matrix K_xz; update the VSG parameters through the adaptive particle swarm optimization algorithm to obtain the virtual moment of inertia J_vsg and virtual damping coefficient D_vsg; Step S43, input the control framework MSC, VSG parameters [J_vsg, D_vsg] and model parameter vector θ_final into the dynamic scene tree generation module; use the kernel density estimation method to construct the frequency disturbance probability distribution p(w); apply the stratified sampling strategy to generate the scene tree S_tree; execute the scene reduction algorithm on S_tree to obtain the simplified scene set S_reduced and the probability set Π_reduced; construct the robust stochastic model predictive control problem, design the objective function J and the constraint condition g; apply the distributed alternating direction multiplier method to solve the optimization problem and obtain the optimal control sequence u*; Step S44, input the optimal control sequence u into the adaptive learning module; construct a value function approximator V(x, θ_V) and a policy function approximator π(x, θ_π); use the deep deterministic policy gradient algorithm to update the value function and policy function parameters; establish an experience replay buffer D to store transfer samples; sample mini-batch from D for policy update; perform strategy fusion on the updated policy function π(x, θ_π) and u, and design an adaptive weight function w(x); calculate the fusion control input u_fusion; apply input constraint processing to u_fusion to obtain the final control input u_final; Step S45, read the historical control effect data H_ctrl from the pre-storage database; input u_final output from 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 use the online gradient descent method to update the decision variables; 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; The strategy is dynamically adjusted based on J_long, and the final long-term frequency modulation strategy L_opt is output.

7. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: The step S5 is specifically as follows: Step S51, input the control strategy MSC, J_vsg, D_vsg, u*, L_opt, and submodule current data I(t) into the distributed coordinated control module; read the communication topology information G and the initial weight matrix W from the pre-stored database; construct a local state vector x_i, including 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 weight, 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; Step S52, input the local control law u_i into the event trigger 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 chaos expansion method to process the parameter uncertainty, 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 use the alternating direction multiplier method to solve and obtain the local optimal control sequence u_i*; 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 Ji 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 frontier generation module; use the NSGA-III algorithm to initialize the population P(0), and design the reference point-based selection operation; 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*; Step S54, input the Pareto optimal solution set P* into the adaptive weight Kriging proxy model; use Latin hypercube sampling to generate the initial sample point set X; evaluate the objective function value Y of the sample point; construct a Kriging model to predict the function value and prediction variance of the unknown point; design the expected improvement criterion EI(x); use the adaptive weight strategy to dynamically adjust the EI weight w_EI; solve the weighted EI maximization problem to obtain a new sampling point x_new; update the sample set and retrain the model to obtain the final Kriging proxy model M_Kriging; Step S55, input the proxy model M_Kriging into the multi-objective Bayesian optimization module; design the acquisition function α(x), and adopt the Thompson sampling strategy to sample from the posterior distribution; solve the sub-optimization 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 multi-objective evaluation results, Pareto optimal solution, sensitivity analysis and optimization suggestions, and output the optimized control parameter set P_opt.

8. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 3, characterized in that: The step S11 is specifically as follows: Step S111, obtaining 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; inputting these data into the information entropy evaluator to calculate the information entropy H(t) of different time windows; using sliding variance analysis to calculate the signal fluctuation index σ(t); constructing an adaptive sampling rate adjustment function f_s(t) based on H(t) and σ(t); applying a variable sampling rate reconstruction algorithm to the original data to obtain the resampled data V_s(t), I_s(t), and θ_s(t); Step S112, construct the V_s(t), I_s(t), θ_s(t) output by step S111 into the observation vector y(t); read the initial state estimate x(0) and the 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 untraceable transformation to generate the sigma point set χ_i and the weight W_i; perform state prediction and observation prediction for each sigma point; and update the gain K(t) by cross-covariance calculation; Obtain the pre-filtered PMU data stream P(t); Step S113, collecting original 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); applying the wavelet dictionary learning algorithm to each parameter to construct a sparse representation basis D_i; using the orthogonal matching pursuit algorithm to calculate the compression coefficient α_i; designing an adaptive compression rate controller to dynamically adjust the compression rate r_i(t) according to the signal importance index; performing distributed compression sampling to obtain a compressed electrolytic cell data stream E(t); Step S114, inputting the pre-filtered PMU data stream P(t) and the compressed electrolyzer data stream E(t) into a multi-dimensional data quality assessment module; The inductor current I_L(t) is collected from 20 sub-modules to construct the current data matrix I(t); the integrity index C_i, accuracy index A_i and consistency index U_i are calculated for P(t), E(t) and I(t) respectively; the fuzzy comprehensive evaluation method is applied, combined with the pre-stored evaluation rule set R, to calculate the comprehensive quality score Q_i; the data is weighted according to Q_i, and a data stream set D_q(t) with quality score is output.

9. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: The step S12 is specifically as follows: Step S121, decomposing the data stream set D_q(t) output from step S11 into PMU data stream P_q(t), electrolyzer data stream E_q(t) and submodule current data stream I_q(t); reading the reference clock source T_ref and the delay compensation parameter set Δt from the pre-stored time synchronization parameter library; applying the bidirectional delay measurement algorithm to each data stream to calculate the transmission delay τ_i and the clock deviation δ_i; using the adaptive Kalman filter to estimate and correct τ_i and δ_i in real time; performing time alignment interpolation on the three types of data streams based on the corrected time parameters to obtain the preliminary aligned data stream D_sync(t); Step S122, reconstruct D_sync(t) outputted 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 abnormality judgment threshold λ(t) using the adaptive threshold algorithm; Reconstruct the tensor decomposition problem through sparse constrained optimization to obtain a low-rank tensor L and a sparse tensor S; Extract the abnormal data point set A from S, and output the normal data tensor L and the abnormal data set A; Step S123, input the abnormal data set A output from step S122 into the multi-layer 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); use the attention mechanism to select similar patterns from the normal data tensor L as a reference; combine the reference pattern and the hidden layer representation, and generate repair data R = f_d(h;θ_d) through the decoder; Apply residual learning to optimize the restoration results and obtain the final restoration dataset R_opt; Step S124, taking the normal data tensor L output by step S122 and the repaired data set R_opt output by step S123 as input; constructing a k-nearest neighbor graph G_k, and calculating the local geometric structure of each data point; designing an adaptive weight matrix W_l, and quantifying the local relationship between data points; using the manifold regularization algorithm to construct a global alignment objective function J = tr(X^TLX), where L is the graph Laplacian matrix; solving the alignment problem by alternating optimization to obtain the aligned data stream D'(t); applying adaptive multi-scale filtering to D'(t), and outputting the filtered data stream D_f(t).

10. The method for controlling the primary frequency regulation of a power system based on modular large-scale hydrogen production power source according to claim 2, characterized in that: The step S22 is specifically as follows: Step S221, input the frequency component f_est in the system state vector x_est into the adaptive bandwidth optimization module; construct a variational mode decomposition objective function, including 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, and dynamically adjust the decomposition layer number 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 corresponding center frequencies {ω_1,ω_2,...,ω_K}; Step S222, input the modal function set {u_k} output by step S221 into the multi-scale analysis module; construct a coarse-grained time series {x_τ} of different scales for each mode u_k; design an adaptive embedding dimension selector, and determine the optimal embedding dimension m_k and 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 the curvature c_k; output the permutation entropy feature vector PE_features; Step S223, input the modal function {u_k} of step S221 and the PE_features of step S222 into the time-frequency analysis module; apply the generalized S transform to each mode 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; Step S224, input the fused feature vector F_fused outputted from step S223 into a deep learning module; construct a deep neural network including multiple residual blocks; load pre-trained parameters θ_pre from a pre-stored model library; process feature data using local response normalization; extract multi-level features through residual connections and skip connections; Design feature selection gating units to dynamically screen important features; apply dropout regularization to prevent overfitting; perform principal component analysis on network output to reduce dimensionality and retain the main components; organize the final frequency perturbation feature matrix F_dist; The step S32 is specifically as follows: 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 to include physical quantities such as voltage and current, and construct the parameter vector θ to include nonlinear mapping coefficients and dynamic parameters; combine x and θ to form an 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 transfer equation f(X_i) to each Sigma point to predict the state at the next moment; 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; Step S322, taking the predicted state x_pred outputted in step S321 as input; constructing a sliding time window W(t), wherein the window length L is adaptively adjusted according to the signal change rate; calculating the first-order moment M_1 to the fourth-order moment M_4 of the data in the window; Design a forgetting factor update rule based on the moment change rate to obtain the adaptive forgetting factor λ(t); Apply the exponentially weighted recursive algorithm to update the moment estimate, M_new = λ(t)M_old + (1-λ(t))M_current; calculate the parameter correction Δθ based on the high-order moment characteristics to obtain the corrected parameter estimate θ_corr; Step S323, inputting the θ_corr outputted from step S322 into the multi-model library; constructing N candidate models {M_1, M_2, ..., M_N} according to historical data and expert knowledge; calculating the likelihood function p(y(k)|M_i) for each model; designing a weight update mechanism based on minimizing the prediction error; using the softmax function to calculate the model weight w_i(k); constructing an interactive multi-model estimation framework to handle the complementarity between models; fusing the parameter estimates of each model to obtain the weighted parameter vector θ_weighted; Step S324, using the θ_weighted outputted in step S323 as the initial parameter; initializing Np particles {θ_p} and their weights {w_p}; applying the state space model to each particle for prediction; Calculate the particle likelihood based on the latest observation data; design an adaptive resampling strategy based on the number of effective particles; construct a moving resampling algorithm to avoid particle degradation; update particle weights and normalize them; calculate the mean and variance of the posterior distribution; obtain the final parameter vector θ_final through iterative optimization; The step S42 is specifically as follows: Step S421, read the historical VSG parameter data set D_hist from the pre-stored database; combine the power grid inertia estimation value 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 an initial kernel feature map φ_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 a representation x_trans that highlights the key features; use the back propagation algorithm to optimize the network parameters and output the transformed feature vector x_opt; 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 induction points Z={z_1,...,z_m}; calculate the cross-covariance matrix K_xz and K_zz; construct a sparse approximate prediction equation to obtain the VSG parameter prediction value [J_pred,D_pred] and prediction variance [v_J,v_D]; Step S423, taking the parameter prediction value [J_pred, D_pred] output by step S422 as the initial solution; initializing the particle swarm P={p_1,..., p_N}; designing an adaptive inertia weight update strategy w(t); constructing a comprehensive performance evaluation function f(J, D), taking into account frequency deviation, rate of change and kinetic energy loss; evaluating the fitness value f(p_i) for each particle; updating the individual optimal position p_best_i and the global optimal position g_best; dynamically adjusting the learning factors c1, c2 according to the performance gradient; applying the speed update formula to avoid local optimality; iteratively optimizing to obtain the parameter optimization value [J_opt, D_opt]; Step S424, input the optimized parameters [J_opt, D_opt] output from 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 the structural singular value analysis method 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 to generate compensation [Δ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; The step S51 is specifically as follows: 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 the adjacency matrix A; calculate the network algebraic connectivity λ_2 as an evaluation index; design a topology optimization objective function J_topo based on graph analysis; use a discrete particle swarm algorithm to search for the optimal topology structure; consider the communication delay τ_ij and the packet loss rate p_loss to construct the communication quality constraint; apply the branch and bound method to solve the discrete optimization problem; obtain the optimized topology structure G_opt and the initial weight matrix W_init; Step S512, take the initial weight matrix W_init outputted in step S511 as input; construct a local state vector x_i, including the inductor current I_i, the output voltage V_i and the 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 projected 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; Step S513, input the Laplacian matrix L outputted from 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 unknown states ẑ_i; construct an 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 state feedback and observer design to obtain the complete controller structure C_i; Step S514, input the controller structure C_i output from step S513 into the distributed optimization module; construct the local cost function J_i and the coupling constraint g_i; design a distributed solution framework based on ADMM; introduce the augmented Lagrangian to process the coupling constraint; construct a prediction model for each subsystem to predict the state trajectory x_i(k+j|k); update local variables, global variables and Lagrange multipliers through iterative alternating optimization; 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.

Images