A method for modeling and parameter optimization of reaction kinetics of an electrocatalytic denitrification process

Abstract

The present application relates to the technical field of computational chemistry and chemical reaction engineering, and particularly relates to a reaction kinetics modeling and parameter optimization method for an electrocatalytic denitrification process, which obtains electrode material crystal structure graph data, nitrogen-containing intermediate species thermodynamic and kinetic benchmark data, and macroscopic experimental data; constructs a potential-dependent physical information graph neural network micro-kinetics model to predict multi-path elementary reaction activation energy and generate a rate equation set; constructs a reaction-mass transfer bidirectional coupling cross-scale kinetics model, which is reduced by concomitant sensitivity analysis; extracts intrinsic catalytic descriptors by manifold learning to construct a physically constrained Gaussian process proxy model; performs multi-objective manifold Bayesian optimization to output an optimal combination of electrode material structure characteristics and operating parameters. The present application realizes the collaborative optimization of material and operating parameters, and takes into account both prediction accuracy and physical feasibility.

Description

Technical Field

[0001] This invention relates to the fields of computational chemistry and chemical reaction engineering, specifically to a method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process. Background Technology

[0002] In terms of kinetic modeling, traditional methods mainly rely on density functional theory (DFT) calculations to obtain the energy barriers of elementary reactions, and then combine them with transition state theory to construct microscopic kinetic models. However, these methods have the following shortcomings: First, DFT calculations are usually performed in a vacuum or static solvation model, making it difficult to accurately capture the modulation effect of dynamic changes in electrode potential on activation energy, resulting in large deviations in kinetic predictions at different potentials; Second, microscopic kinetic models usually assume that the electrode surface is a uniform constant potential surface, ignoring the spatial non-uniformity of the electric field distribution within the electric double layer and the influence of local pH evolution on the reaction rate, which significantly reduces the prediction accuracy of the model at high current densities; Third, existing methods for coupling microscopic kinetics and macroscopic mass transfer mostly adopt a one-way coupling strategy, that is, solving the microscopic kinetics first and then substituting it into the mass transfer equation, lacking a two-way feedback mechanism between reaction and mass transfer, and failing to accurately describe the dynamic regulation of reaction selectivity by concentration polarization and local pH effects.

[0003] In terms of surrogate model construction and optimization, traditional response surface methodology and Kriging surrogate models directly use the original operating parameter space as input. When the operating parameters have high dimensionality and complex nonlinear coupling relationships, the prediction accuracy and generalization ability of the surrogate model drop sharply. Furthermore, existing Bayesian optimization methods do not consider the geometric manifold structure of candidate solutions in the material feature space during sampling, resulting in low sampling efficiency in the material-operating parameter joint design space and a tendency to get trapped in local optima. Simultaneously, purely data-driven surrogate models lack physical constraints and may output invalid solutions that violate thermodynamic conservation laws, reducing the reliability and reproducibility of the optimization results.

[0004] In terms of parameter optimization, existing methods typically separate electrode material selection from operating parameter optimization, first using trial and error to select materials and then optimizing operating conditions. This makes it difficult to achieve synergistic optimization of material structural characteristics and operating parameters, and fails to fully realize the performance improvement potential brought about by material-operating parameter matching. In addition, existing optimization methods lack quantitative analysis of the sensitivity to elementary reactions, and model all elementary reactions with the same precision, resulting in a huge computational burden that is difficult to meet the needs of efficient iterative optimization.

[0005] Therefore, there is an urgent need for a systematic approach that can achieve accurate prediction of potential-dependent microdynamics, bidirectional coupling simulation of reaction-mass transfer, efficient construction of physical constraint surrogate models, and synergistic optimization of material-operation parameters, in order to overcome the above-mentioned shortcomings of existing technologies. Summary of the Invention

[0006] The purpose of this invention is to provide a method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] A method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process includes the following steps:

[0009] S1: Obtain basic data of the electrocatalytic denitrification reaction system. The basic data includes crystal structure diagram data of the target electrode material, thermodynamic and kinetic baseline data of nitrogen-containing intermediate species at different electrode potentials, and macroscopic reaction kinetics and mass transfer experimental data.

[0010] S2: Construct a potential-dependent physical information graph neural network microdynamic model. Input crystal structure diagram data and double-layer electric field physical constraints into the physical information graph neural network to predict the activation energy of elementary reactions under multi-path competition. Combine the transition state theory to generate a set of potential-dependent microdynamic rate equations. The multi-path competition includes the nitrogen denitrification to nitrogen generation path, the ammonia reduction to nitrogen generation path, and the hydrogen evolution side reaction path.

[0011] S3: Construct a reaction-mass transfer bidirectional coupled cross-scale kinetic model, and couple the microscopic kinetic rate equations as source terms into the one-dimensional Nernst-Planck-Poisson mass transfer equation to simulate the local pH dynamics and concentration polarization effect within the electric double layer; calculate the global sensitivity index of each elementary reaction to the macroscopic denitrification index based on adjoint sensitivity analysis, and reduce the order of insensitive elementary reactions to obtain a reduced-order cross-scale kinetic model.

[0012] S4: Construct a physically constrained manifold surrogate model. Based on the reduced-order cross-scale dynamic model, generate a multi-dimensional sample set containing material diagram feature descriptors and macroscopic operating parameters. Use the manifold learning algorithm to extract intrinsic catalytic descriptors. Construct a physically constrained Gaussian process surrogate model with intrinsic catalytic descriptors and macroscopic operating parameters as inputs and denitrification rate, nitrogen selectivity and energy consumption as outputs.

[0013] S5: Perform multi-objective manifold Bayesian optimization, and use the acquisition function based on the Pareto front expectation hypervolume improvement to iteratively optimize the physical constraint Gaussian process surrogate model. In the optimization process, a thermodynamic consistency penalty function is introduced to eliminate invalid solutions that violate physical conservation laws, and the optimal electrode material structure characteristics and macroscopic operating parameter combination are output.

[0014] As can be seen from the technical solution provided by the present invention above, the reaction kinetics modeling and parameter optimization method for an electrocatalytic denitrification process provided by the present invention has the following beneficial effects:

[0015] By constructing a potential-dependent physical information graph neural network microdynamic model, the physical constraints of the Gouy-Chapman-Stern double-layer electric field are embedded into the message passing layer of the graph neural network. Adaptive feature injection of electrode potential is achieved through a gating attention mechanism, enabling the model to accurately capture the nonlinear modulation effect of electrode potential changes on the activation energy of each elementary reaction. Compared with the traditional DFT linear scaling relationship correction method, it significantly improves the prediction accuracy of activation energy of multi-path competitive reactions under different electrode potentials. At the same time, thermodynamic consistency correction is performed on the predicted activation energy based on thermodynamic cycle verification conditions to ensure the thermodynamic self-consistency of the kinetic parameters.

[0016] By constructing a reaction-mass transfer bidirectional coupled cross-scale kinetic model, the microscopic kinetic rate equations are coupled as interface reaction source terms into the one-dimensional Nernst-Planck-Poisson mass transfer equation, realizing bidirectional real-time feedback between microscopic surface reactions and macroscopic mass transfer processes. This model can accurately simulate the dynamic evolution of local pH and concentration polarization effects within the electric double layer. Based on this, the global sensitivity index of each elementary reaction to macroscopic denitrification indicators is calculated using the adjoint sensitivity analysis method. Quasi-steady-state approximation is performed on insensitive elementary reactions to reduce their order. This significantly reduces the dimensionality of the state variables in the cross-scale model while ensuring the model's prediction accuracy, providing a high-fidelity simulation framework that balances accuracy and efficiency for subsequent iterative optimization.

[0017] By employing a manifold learning algorithm, low-dimensional intrinsic catalytic descriptors are extracted from high-dimensional multi-source data containing material graph feature descriptors and macroscopic operating parameters. This effectively eliminates redundant information and noise interference in the original feature space, achieving a compact representation of material features and operating conditions in the low-dimensional manifold space. Based on this, a physically constrained Gaussian process surrogate model is constructed. By introducing a material similarity kernel, an operating parameter kernel, and a cross-coupling kernel into the kernel function, and applying thermodynamic consistency and mass conservation constraints as regularization penalty terms, the surrogate model can strictly adhere to physical conservation laws while achieving high prediction accuracy, avoiding physically infeasible solutions that may occur in purely data-driven models.

[0018] By performing multi-objective manifold Bayesian optimization and employing a collection function based on the product of the Pareto front expected hypervolume improvement and the manifold distance metric, an efficient global search is conducted in the low-dimensional intrinsic catalytic descriptor space. Simultaneously, a thermodynamic consistency penalty function is introduced to automatically eliminate invalid candidate solutions that violate physical conservation laws during the optimization process. This achieves synergistic optimization of electrode material structural features and macroscopic operating parameters, which can simultaneously maximize steady-state denitrification rate and nitrogen Faraday selectivity, minimize denitrification per unit energy consumption, and reconstruct the optimal intrinsic catalytic descriptor into the crystal structure features of the electrode material that can guide experimental synthesis through inverse mapping. Attached Figure Description

[0019] Figure 1 This is a schematic diagram of the overall structure of the reaction kinetics modeling and parameter optimization method for an electrocatalytic denitrification process according to the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0021] To better understand the above technical solutions, the following will provide a detailed description of the technical solutions in conjunction with the accompanying drawings and specific embodiments.

[0022] like Figure 1 As shown in the figure, this invention provides a method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process, including the following steps:

[0023] S1: Obtain basic data of the electrocatalytic denitrification reaction system. The basic data includes crystal structure diagram data of the target electrode material, thermodynamic and kinetic baseline data of nitrogen-containing intermediate species at different electrode potentials, and macroscopic reaction kinetics and mass transfer experimental data.

[0024] S2: Construct a potential-dependent physical information graph neural network microdynamic model. Input crystal structure diagram data and double-layer electric field physical constraints into the physical information graph neural network to predict the activation energy of elementary reactions under multi-path competition. Combine the transition state theory to generate a set of potential-dependent microdynamic rate equations. The multi-path competition includes the nitrogen denitrification to nitrogen generation path, the ammonia reduction to nitrogen generation path, and the hydrogen evolution side reaction path.

[0025] S3: Construct a reaction-mass transfer bidirectional coupled cross-scale kinetic model, and couple the microscopic kinetic rate equations as source terms into the one-dimensional Nernst-Planck-Poisson mass transfer equation to simulate the local pH dynamics and concentration polarization effect within the electric double layer; calculate the global sensitivity index of each elementary reaction to the macroscopic denitrification index based on adjoint sensitivity analysis, and reduce the order of insensitive elementary reactions to obtain a reduced-order cross-scale kinetic model.

[0026] S4: Construct a physically constrained manifold surrogate model. Based on the reduced-order cross-scale dynamic model, generate a multi-dimensional sample set containing material diagram feature descriptors and macroscopic operating parameters. Use the manifold learning algorithm to extract intrinsic catalytic descriptors. Construct a physically constrained Gaussian process surrogate model with intrinsic catalytic descriptors and macroscopic operating parameters as inputs and denitrification rate, nitrogen selectivity and energy consumption as outputs.

[0027] S5: Perform multi-objective manifold Bayesian optimization, and use the acquisition function based on the Pareto front expectation hypervolume improvement to iteratively optimize the physical constraint Gaussian process surrogate model. In the optimization process, a thermodynamic consistency penalty function is introduced to eliminate invalid solutions that violate physical conservation laws, and the optimal electrode material structure characteristics and macroscopic operating parameter combination are output.

[0028] In this embodiment, the core function of step S1 is to systematically collect three types of basic data from the electrocatalytic denitrification reaction system. This provides a complete technical path for constructing crystal structure diagram data of the target electrode material, obtaining thermodynamic and kinetic benchmark data of nitrogen-containing intermediate species at different electrode potentials, and collecting macroscopic reaction kinetics and mass transfer experimental data. This lays the data foundation for the subsequent training of the physical information graph neural network micro-dynamic model and the construction of the cross-scale kinetic model. The detailed steps are as follows:

[0029] Step S1-1: Construction of crystal structure data for the target electrode material:

[0030] Based on the chemical composition and crystal phase information of the target electrode material, an initial crystal model of the target electrode material is constructed. The initial crystal model is constructed based on lattice parameters, space group symmetry, and atomic occupancy information from experimental characterization data or material databases. Density functional theory is used to optimize the geometric structure of the initial crystal model, with the goal of minimizing the total ground state energy of the crystal model. During the optimization process, the conjugate gradient algorithm is used to iteratively relax the atomic coordinates and lattice parameters. The optimization is considered complete when the residual force on each atom is less than the preset force convergence threshold and the total energy change of the system is less than the preset energy convergence threshold.

[0031] Atomic coordinates and bonding information are extracted from the optimized crystal model. Atomic coordinates are the three-dimensional Cartesian space positions of each atom in the unit cell reference system in the optimized crystal model. Bonding information is determined by the ratio of the distance between atoms to the sum of the covalent radii of the corresponding atoms. When the distance between two atoms is less than a preset multiple of the sum of their covalent radii, it is determined that there is a chemical bond between the atom pair.

[0032] Based on the extracted atomic coordinates and bonding information, crystal structure diagram data is constructed. Each atom in the optimized crystal model is defined as a graph node, and a node feature vector is assigned to each node. The node feature vector contains atomic number features and electronegativity features. The atomic number features represent elemental information using one-hot encoding or continuous encoding, and the electronegativity features represent the electron affinity of the atom using the Pauling electronegativity scale. Each chemical bond in the optimized crystal model is defined as a graph edge, and an edge feature vector is assigned to each edge. The edge feature vector contains bond length features and bond order features. The bond length feature is the Euclidean distance between the atoms at the two ends of the bond, and the bond order feature is represented by the chemical bond strength calculated by bond valence theory. The set of node feature vectors, the set of edge feature vectors, and the graph topological connections are integrated to generate crystal structure diagram data of the target electrode material.

[0033] Step S1-2: Acquisition of thermodynamic and kinetic baseline data for nitrogen-containing intermediate species at different electrode potentials:

[0034] Based on the crystal structure data of the target electrode material, a solid-liquid interface double layer calculation model was constructed, which includes an implicit solvation model and an explicit water molecule layer. The implicit solvation model uses a conductor-like shielding model to describe the continuous medium polarization effect of the bulk solvent water. The explicit water molecule layer constructs one or two layers of ordered water molecule arrangement above the adsorption sites on the electrode surface to describe the hydrogen bond network structure and local solvation effect of the interfacial water. The solid-liquid interface double layer calculation model is constructed by spatial superposition of an atomic-level surface model and the solvation environment. The surface model contains a vacuum layer of sufficient thickness along the normal direction perpendicular to the electrode surface to accommodate the solvation environment.

[0035] Different electrode potential conditions are introduced into the solid-liquid interface double-layer calculation model. The electrode potential is adjusted by changing the total charge of the system or by applying an external electric field to the double-layer region. The electrode potential range covers the target potential range of the electrocatalytic denitrification reaction. Under each electrode potential condition, the adsorption free energy of nitrogen-containing intermediate species at each highly symmetric adsorption site on the surface of the target electrode material is calculated. The formula for calculating the adsorption free energy is as follows: ,in, The adsorption free energy of nitrogen-containing intermediate species, The total energy of the system after nitrogen-containing intermediate species are adsorbed onto the surface of the target electrode material is denoted as . The total energy of the clean surface of the target electrode material. This serves as a reference energy for isolated nitrogen-containing intermediates in the gas or liquid phase. This is the zero-point vibration energy correction term. The reaction temperature. The entropy change before and after adsorption is calculated; the adsorption free energy of each nitrogen-containing intermediate species at each adsorption site under each electrode potential is summarized to obtain the thermodynamic baseline data of nitrogen-containing intermediate species under different electrode potentials.

[0036] Simultaneously, the micro-motion elastic band method is used to search for the lowest energy path through which nitrogen-containing intermediate species transform into the next intermediate species or transition state along the reaction path. This method inserts multiple interpolated configurations between the reactant and product states as initial reaction paths, constrains the spacing between adjacent configurations using spring forces, and minimizes the energy of each configuration along the potential energy surface, eventually converging to the lowest energy path. The configuration corresponding to the highest potential energy point on the converged lowest energy path is extracted as the transition state, and the energy difference between the transition state and the reactant state is the activation energy barrier for each elementary reaction. Combined with the electrode potential's shift correction of the reaction free energy diagram, the correction of the activation energy barrier by the electrode potential adopts a linear scaling relationship based on the electron transfer coefficient. ,in, Electrode potential The modified activation barrier below, As the uncorrected activation barrier at the reference electrode potential, The electron transfer coefficient, The basic charge is used; the activation energy barriers of each elementary reaction under each electrode potential are summarized to obtain the kinetic baseline data of nitrogen-containing intermediate species under different electrode potentials; the thermodynamic baseline data and the kinetic baseline data are integrated into thermodynamic and kinetic baseline data of nitrogen-containing intermediate species under different electrode potentials.

[0037] Step S1-3: Data acquisition for macroscopic reaction kinetics and mass transfer experiments:

[0038] Based on the thermodynamic and kinetic baseline data of nitrogen-containing intermediates at different electrode potentials, the thermodynamically feasible potential window and the kinetic rate-determining step of the electrocatalytic denitrification reaction were determined. The thermodynamically feasible potential window is the electrode potential range where the adsorption free energy of all elementary reactions is negative or close to zero, and the kinetic rate-determining step is the elementary reaction step with the highest activation energy barrier in all elementary reaction sequences. The potential ranges corresponding to the thermodynamically feasible potential window and the kinetic rate-determining step are used as the basis for setting the operating parameters of the high-throughput microfluidic electrochemical testing array.

[0039] Under set operating parameters and fluid shear force conditions, the dynamic response curves of product concentration versus time during the electrocatalytic denitrification reaction were continuously acquired using a high-throughput microfluidic electrochemical testing array. The high-throughput microfluidic electrochemical testing array contains multiple parallel microreaction channels, each channel integrating a three-electrode system consisting of a working electrode, a reference electrode, and a counter electrode. The working electrode is a microelectrode with a surface-modified target electrode material. The fluid shear force is controlled by adjusting the volumetric flow rate of the electrolyte in the microchannel, and the volumetric flow rate range is set according to the sensitivity of the kinetic rate-determining step to the mass transfer effect.

[0040] In-situ analysis of the dynamic response curves was performed using online mass spectrometry and ion chromatography. Online mass spectrometry introduced volatile products from the reaction solution at the microchannel outlet into the mass ionization source in real time via a membrane injection interface, detecting the change in the mass-to-charge ratio of nitrogen products and dinitrogen intermediates over time. Ion chromatography periodically collected reaction solution samples from the microchannel outlet using an automatic injection valve, separating and quantitatively detecting the concentrations of nitrate ions, nitrite ions, ammonium ions, and hydroxylamine intermediates in the liquid phase. The macroscopic apparent reaction rate constant, limiting current density, and diffusion layer thickness were extracted from the analyzed dynamic response curves. The macroscopic apparent reaction rate constant was obtained by fitting the slope of the initial linear segment of the product concentration versus time curve; the limiting current density was determined by the average value of the plateau region of the current density versus electrode potential curve; and the diffusion layer thickness was calculated using the linear relationship between the limiting current density and the reactant phase concentration and diffusion coefficient. The macroscopic apparent reaction rate constant, limiting current density, and diffusion layer thickness were then combined to generate macroscopic reaction kinetics and mass transfer experimental data.

[0041] In this embodiment, the core function of step S2 is to use the crystal structure diagram data of the target electrode material obtained in step S1 and the physical constraints of the double-layer electric field as input to construct a potential-dependent physical information graph neural network microdynamic model. This model predicts the activation energies of elementary reactions under three competing pathways: nitrogen generation through denitrification, ammonia generation through reduction, and hydrogen evolution as a side reaction. It also combines transition state theory to generate a potential-dependent set of microdynamic rate equations, providing microdynamic source terms for the subsequent construction of a cross-scale dynamic model. The detailed steps are as follows:

[0042] Step S2-1: Construction of the physical information graph neural network skeleton:

[0043] Based on the crystal structure map data of the target electrode material, a physical information graph neural network skeleton including a node feature embedding layer, an edge feature embedding layer, and a multi-physics message passing layer is constructed.

[0044] The node feature embedding layer maps the node feature vectors of graph nodes to high-dimensional atomic representation vectors; let the node... The node feature vector is Its dimensions are If the atomic number and electronegativity characteristics are included, then the mapping formula for the high-dimensional atomic representation vector is: ,in, For graph nodes The initial high-dimensional atomic representation vector after mapping through the node feature embedding layer, The learnable weight matrix for the node feature embedding layer. This is the bias vector for the node feature embedding layer. It is a non-linear activation function;

[0045] The edge feature embedding layer maps the edge feature vectors of the graph edges to high-dimensional chemical bond representation vectors; let the graph nodes be connected. With graph nodes The edge eigenvectors of the graph are Its dimensions are If the bond length and bond order features are included, then the mapping formula for the high-dimensional chemical bond characterization vector is: ,in, This is the initial high-dimensional chemical bond representation vector after the graph edges are mapped through the edge feature embedding layer. The learnable weight matrix for the edge feature embedding layer. The bias vector of the edge feature embedding layer;

[0046] In a multiphysics message passing layer, graph nodes In the The layer's representation vector is updated through message passing with neighboring graph nodes. The message passing formula is: ,in, For graph nodes In the The high-dimensional atomic representation vector of the layer, For graph nodes In the The high-dimensional atomic representation vector of the layer, For the first The high-dimensional chemical bond characterization vector of the layer, For graph nodes The set of nodes in the adjacency graph. For the first Layer node self-transformation weight matrix, For the first Layer edge message transformation weight matrix, This is an element-wise multiplication operation. For graph nodes In the The updated high-dimensional atomic representation vector after layer update;

[0047] Step S2-2: Embedding the physical constraints of the double-layer electric field:

[0048] A double-layer electric field physical constraint module is constructed in the multiphysics message passing layer; the double-layer electric field physical constraint module calculates the local electric field intensity distribution and potential decay curve on the surface of the target electrode material based on the Gouy-Chapman-Stern double-layer theory;

[0049] Let the normal direction of the electrode surface be... The shaft, the position of the electrode surface is The position of the bulk solution is Then the potential decay curve satisfy: ,in, Distance from electrode surface The potential value at that point, The potential of the electrode surface. This is the reciprocal of the Debye shield length. The calculation formula is: , Let Avogadro's constant be 1. Ionic strength, The vacuum permittivity, The relative permittivity of the solvent, Boltzmann's constant;

[0050] Local electric field intensity distribution The following is obtained by taking the negative gradient of the potential decay curve: ,in, It is the unit normal vector perpendicular to the electrode surface and pointing towards the bulk solution;

[0051] The local electric field intensity distribution is superimposed onto the high-dimensional atomic representation vector as the external field perturbation vector of each graph node; let the graph node... The corresponding atom has the following coordinates above the electrode surface: Then the external field perturbation vector of the node in the graph for: ,in, The superimposed enhanced graph node features are learnable electric field response vectors. for: ;

[0052] The potential decay curve is used as a weighting correction factor for each graph edge to correct the conventional chemical bond characterization vector; for connecting graph nodes... Graph Nodes Edges, weight correction factor for: ;

[0053] in, For graph nodes The electric potential value at the corresponding atomic position, For graph nodes The electric potential value at the corresponding atomic position, The reference potential normalization constant; the modified enhanced graph edge features for: ;

[0054] Enhanced graph node features and enhanced graph edge features Replace the corresponding terms in the original message passing formula to generate enhanced graph node features and enhanced graph edge features that include electric field physical constraints;

[0055] Step S2-3: Potential-dependent feature injection and atomic-level reactivity characterization:

[0056] A potential-dependent feature injection mechanism is constructed to encode different electrode potential conditions into potential embedding vectors. The potential embedding vectors are then adaptively fused with enhanced graph node features containing electric field physical constraints through a gated attention mechanism.

[0057] First, set the electrode potential Position encoding is mapped to potential embedding vectors. : ,in, For position encoding functions, For potential embedding vector, Let be the dimension of the potential embedding vector. For dimension indexing;

[0058] Gated attention mechanisms achieve adaptive fusion by calculating the attention weights between the potential embedding vector and the features of the enhanced graph nodes; assuming graph nodes... The enhanced graph node features are Then the gating attention weight for: ,in, This is a learnable attention query vector. The weight matrix is ​​the feature transformation matrix for graph nodes. The potential embedding transformation weight matrix, The set of all graph nodes. For graph nodes The corresponding enhanced graph node features include physical constraints on the double-layer electric field. For graph nodes The corresponding enhanced graph node features include physical constraints on the double-layer electric field;

[0059] Potential-dependent atomic-level reactivity characterization The result obtained through gating fusion is: ,in, For graph nodes The corresponding potential-dependent atomic-level reactivity characterization vector;

[0060] Step S2-4: Prediction of activation energies of elementary reactions under multi-path competition:

[0061] Characterizing potential-dependent atomic-level reactivity Input multi-task pre-contact; the multi-task pre-contact includes a reaction pathway classification subnetwork and an activation energy regression subnetwork;

[0062] The reaction path classification subnetwork first aggregates the atomic-level reaction activity representations of all graph nodes through graph-level pooling operations to obtain graph-level representation vectors. : ;

[0063] Input the graph-level representation vector into a classification multilayer perceptron, and output the probability distribution of each response path: ,in, It is a three-element probability vector. The probability of the nitrogen generation pathway occurring during denitrification. To reduce the probability of the ammonia generation pathway occurring, This represents the probability of the hydrogen evolution side reaction pathway occurring. The normalized exponential function maps the input vector to a probability distribution where each element takes a value between 0 and 1, and the sum is 1. and The learnable weight matrix for the classification subnetwork. and This is the bias vector for the classification subnetwork;

[0064] The activation energy regression subnetwork, for each elementary reaction in each reaction pathway, takes the atomic-level reaction activity characterization of the corresponding active site as input and outputs the predicted forward activation energy of that elementary reaction; let the elementary reaction... The set of nodes corresponding to the active sites in the graph is The predicted activation energy of this elementary reaction is... for: ,in, The learnable weight vector of the regression subnetwork. For the bias scalar of the regression subnetwork;

[0065] Similarly, the predicted value of reverse activation energy The elementary reaction sequences of the denitrification to ammonia generation pathway, the reduction to nitrogen generation pathway, and the hydrogen evolution side reaction pathway were identified by the reaction pathway classification subnetwork, and the forward and reverse activation energies of each elementary reaction sequence at the corresponding electrode potential were predicted by the activation energy regression subnetwork.

[0066] Step S2-5: Thermodynamic Cycle Verification and Activation Energy Correction:

[0067] Based on the forward activation energy and reverse activation energy, and combined with the adsorption free energy of the nitrogen-containing intermediate species on the surface of the target electrode material obtained in step S1, thermodynamic cycle verification conditions for each elementary reaction are constructed.

[0068] For elementary reactions Thermodynamic cycle consistency requires positive activation energy. With reverse activation energy The difference is equal to the change in reaction free energy. : ,in, Elementary reaction The change in reaction free energy is calculated from the adsorption free energy of each nitrogen-containing intermediate species in step S1;

[0069] Introducing a thermodynamic consistency correction term Symmetrical corrections are made to the forward and reverse activation energies: ; ;

[0070] in, This is the error in the closure of the thermodynamic cycle. This is the corrected positive activation energy. This is the corrected reverse activation energy;

[0071] Step S2-6: Generation of potential-dependent microscopic kinetic rate equations:

[0072] Based on the modified positive activation energy and the corrected reverse activation energy The microscopic reaction rate constants of each elementary reaction are calculated using transition state theory.

[0073] For elementary reactions Its forward microscopic reaction rate constant and the reverse microscopic reaction rate constant They are respectively: ; ;

[0074] in, Let be Planck's constant. The molar gas constant, This is the positive activation energy after thermodynamic cycle verification and correction. The reverse activation energy is corrected by thermodynamic cycling; the microscopic reaction rate constant is a function of electrode potential, local electric field strength, and temperature.

[0075] Based on the microscopic reaction rate constants of each elementary reaction, microscopic kinetic rate equations were constructed for the denitrification to nitrogen generation pathway, the reduction to ammonia pathway, and the hydrogen evolution side reaction pathway, respectively; for nitrogen-containing intermediate species... Its surface coverage change rate is: ,in, Nitrogen-containing intermediate species Coverage on the surface of the target electrode material It is the set of all elementary reactions. Nitrogen-containing intermediate species In elementary reactions The stoichiometric coefficients in the figure are represented by positive values ​​for products and negative values ​​for reactants. Elementary reaction The net reaction rate;

[0076] Elementary reactions Net reaction rate for: ,in, Elementary reaction The set of reactants, Elementary reaction The collection of products;

[0077] The microkinetic rate equations for the denitrification to nitrogen generation pathway, the reduction to ammonia pathway, and the hydrogen evolution side reaction pathway are coupled to construct a potential-dependent set of microkinetic rate equations under multi-path competition: ,in, Let this be the coverage vector for all nitrogen-containing intermediate species. This represents the total number of nitrogen-containing intermediate species. This is the vector-valued function of the microdynamic rate equations; thus, the construction of the potential-dependent physical information graph neural network microdynamic model is complete.

[0078] In this embodiment, the core function of step S3 is to use the potential-dependent micro-kinetic rate equations generated in step S2 as the interface reaction source terms, bidirectionally couple them to the one-dimensional Nernst-Planck-Poisson mass transfer equation, simulate the dynamic evolution of local pH value and concentration polarization effect within the electric double layer, and calculate the global sensitivity index of each elementary reaction to macroscopic denitrification indicators based on adjoint sensitivity analysis. Insensitive elementary reactions are then reduced in order, ultimately yielding a reduced-order cross-scale kinetic model. This provides a cross-scale simulation framework that balances accuracy and computational efficiency for the subsequent construction of surrogate models. The detailed steps are as follows:

[0079] Step S3-1: Construction of the interface reaction source term vector:

[0080] Based on the reactant consumption rate and product formation rate of each elementary reaction in the micro-kinetic rate equation set generated in step S2, an interfacial reaction source term vector containing nitrogen-containing intermediate species, protons, hydroxide ions and electrolyte ions is constructed.

[0081] Suppose that the system of microscopic mathematical rate equations contains Nitrogen-containing intermediate species, protons hydroxide ions as well as The total number of electrolyte ions is defined as follows: ; Interface reaction source term vector The dimension is , its first The components are: ,in, For species Net formation rate at the electrode-solution interface, in units of , For species In elementary reactions The stoichiometric coefficients in the figure are represented by positive values ​​for products and negative values ​​for reactants. Elementary reaction Net positive reaction rate Let be the set of all elementary reactions; the interfacial reaction source term vector is represented as an integral table. ;

[0082] Step S3-2: Construction of the one-dimensional spatial discretized mass transfer computational domain:

[0083] A one-dimensional spatial discretization mesh is constructed along the direction perpendicular to the normal to the electrode surface; let the position of the electrode surface be... The outer boundary of the diffusion layer is located at ,in The thickness of the diffusion layer is determined by the macroscopic reaction kinetics and mass transfer experimental data from step S1; the interval is... Uniformly discretized There are 1 grid layer, with a grid layer spacing of . , No. The center coordinates of each grid layer are ,in ;

[0084] Defining concentration field variables in each grid layer and electric potential field variables Concentration field variables Indicates the first In each grid layer Concentration of each species, in units of Electric potential field variables Indicates the first The potential value at the center of each grid layer, in units of The concentration field variables and electric potential field variables of each grid layer are integrated to obtain a one-dimensional spatial discretized mass transfer computation domain.

[0085] Step S3-3: Construction of a cross-scale dynamic model for bidirectional coupling of reaction and mass transfer:

[0086] Interface reaction source term vector As a boundary flux source term, it is coupled into the Nernst-Planck mass transfer equation of the boundary grid layer near the electrode surface in the one-dimensional spatial discretized mass transfer computation domain;

[0087] For the Species in each grid layer Its Nernst-Planck mass transfer equation is: ,in, For the first Species in each grid layer The concentration, in units of , For species The diffusion coefficient, in units of , For species The number of charges, It is Faraday's constant. The molar gas constant, The reaction temperature;

[0088] In the boundary grid layer near the electrode surface At this point, an interface reaction source term is introduced as a boundary flux condition: ;

[0089] This condition indicates the species. The diffusion flux at the electrode surface is equal to the net consumption or net generation rate at the interface due to the electrocatalytic reaction.

[0090] The charge transfer rate involved in the electron transfer steps of each elementary reaction is coupled as a source term to the Poisson potential equation in the one-dimensional discretized mass transfer calculation domain; the Poisson potential equation is: ,in, The vacuum permittivity, The relative permittivity of the solvent; in the boundary mesh layer At this point, the surface charge density of the electrode Determined by the net charge transfer rate through the electron transfer step: ,in, For a subset of elementary reactions involving electron transfer steps, Elementary reaction The number of electrons transferred;

[0091] By integrating the above coupling equations, an initial cross-scale dynamic model of bidirectional reaction-mass transfer coupling is obtained;

[0092] Step S3-4: Implicit time integration and alternating iterative solution:

[0093] The initial multiscale kinetic model of the reaction-mass transfer bidirectional coupling is solved by time-stepping using an implicit time integration method; in the first... There are 1 time steps, and the time step size is denoted as . , from concentration field variable Time advances to time;

[0094] Within each time step, the concentration field variable and the potential field variable are updated alternately based on the Newton-Raphson iterative method; let the th time step be... The approximate solutions for the concentration field variable and the potential field variable in the next iteration are respectively and Then the first The update formula for the next iteration is: ; ;

[0095] in, Let be the residual vector function of the concentration field. Let be the residual vector function of the electric potential field. for For the Jacobian matrix of the concentration field, for Jacobian matrix for the electric potential field;

[0096] The time step is considered convergent when the residual norm of two consecutive iterations satisfies the following convergence criterion: ,in, The preset convergence threshold is used;

[0097] Step S3-5: Simulation of local pH dynamics and concentration polarization effect within the electric double layer:

[0098] Based on the proton concentration of each grid layer in the converged concentration field distribution Calculate the dynamic evolution curve of local pH value within the electric double layer over time; Each grid layer at time The local pH value is: ,in, For the first Each grid layer at time The proton concentration, The standard concentration is set to a value of [value missing]. The local pH values ​​of each grid layer at each time point are summarized and plotted to obtain the dynamic evolution curve of the local pH value within the double layer over time.

[0099] Calculate the concentration polarization overpotential based on the concentration gradient of nitrogen-containing intermediate species in each grid layer of the converged concentration field distribution. : ,in, The total number of electrons transferred. Nitrogen-containing intermediate species The concentration of the grid layer on the electrode surface, Nitrogen-containing intermediate species The concentration in the bulk solution; summarizing the dynamic evolution curves and concentration polarization overpotentials, and completing the construction of a cross-scale kinetic model of bidirectional coupling of reaction and mass transfer;

[0100] Step S3-6: Accompanying sensitivity analysis and global sensitivity index calculation:

[0101] Steady-state denitrification rate extracted from a reaction-mass transfer bidirectional coupling cross-scale kinetic model. Nitrogen Faraday selectivity and nitrogen removal per unit energy consumption As a macroscopic nitrogen removal indicator ; using the microscopic reaction rate constants of each elementary reaction For the perturbation variable, construct the adjoint equation;

[0102] Set macroscopic nitrogen removal index For microscopic reaction rate constant The gradient can be calculated using the adjoint method, introducing an adjoint variable vector. The adjoint equation is: ,in, The system of microscopic dynamic rate equations under steady-state conditions vector-valued functions, Let be the coverage vector of nitrogen-containing intermediate species. Let be the Jacobian matrix of the steady-state equation with respect to the coverage vector;

[0103] Solving the adjoint equation yields the microscopic reaction rate constants for each elementary reaction. The corresponding distribution of adjoint variables; and the reaction flux of each elementary reaction in the microkinetic rate equations. Perform inner product calculations to obtain the effect of each elementary reaction on the macroscopic denitrification index. Local sensitivity coefficient: ,in, Elementary reaction Macroscopic denitrification indicators The local sensitivity coefficient;

[0104] Within the preset electrode potential variation range and temperature variation range The activation energies of each elementary reaction were determined using the Latin hypercube sampling method. and pre-exponential factors Perform random sampling to generate A set of kinetic parameter samples was obtained; for each set of kinetic parameter samples, a reaction-mass transfer bidirectional coupled cross-scale kinetic model was run to obtain the corresponding macroscopic denitrification index sample set. ;

[0105] The local sensitivity coefficients are weighted and integrated with the macroscopic denitrification index sample set to calculate the global sensitivity index of each elementary reaction to the steady-state denitrification rate, ammonia Faraday selectivity, and denitrification per unit energy consumption; elementary reactions Macroscopic denitrification indicators Global sensitivity index for: ,in, These correspond to three indicators: steady-state denitrification rate, ammonia Faraday selectivity, and denitrification per unit energy consumption. For the first Elementary reactions under group samples For indicators The local sensitivity coefficient;

[0106] Step S3-7: Insensitive elementary reaction order reduction and generation of reduced-order cross-scale kinetic model:

[0107] Global sensitivity index Below the preset sensitivity threshold The elementary reactions were labeled as insensitive elementary reactions; the insensitive elementary reactions were subjected to a quasi-steady-state approximation, and the change rate of surface coverage of the nitrogen-containing intermediate species corresponding to the insensitive elementary reactions was set to zero.

[0108] Let the elementary reaction be... The rate of change in the coverage of nitrogen-containing intermediates for reactions labeled as insensitive elementary reactions is: ,in, For the first Nitrogen-containing intermediate species in the first Stoichiometric coefficients in elementary reactions It is a set of insensitive elementary reactions; the equation is transformed into an algebraic constraint equation, directly eliminating the corresponding surface coverage differential equation from the differential equation system;

[0109] After eliminating the surface coverage differential equation corresponding to the insensitive elementary reaction, the dimension of the original microscopic kinetic rate equations changes from... Down to , received a reduction

[0110] The set of microscopic dynamic rate equations: ,in, This is the reduced-order coverage vector of sensitive nitrogen-containing intermediate species. For the reduced-order microdynamic rate equations;

[0111] The reduced-order microdynamic rate equations are recoupled to the one-dimensional Nernst-Planck-Poisson mass transfer equations, and the interfacial reaction source term vector is updated as follows: It only includes the contribution of sensitive elementary reactions retained after the order reduction; repeat the coupling and solution process from step S3-3 to step S3-5 to verify that the relative deviations between the reduced-order model and the original transscale kinetic model in the three indicators of steady-state denitrification rate, nitrogen Faraday selectivity and denitrification amount per unit energy consumption are all less than the preset allowable deviation threshold; if the deviation requirements are met, the reduced-order model passes the verification and the reduced-order transscale kinetic model is obtained.

[0112] In this embodiment, the core function is to use the reduced-order, multi-scale dynamic model constructed in step S3 as a high-fidelity data generator. Through manifold learning, low-dimensional intrinsic catalytic descriptors are extracted from the high-dimensional material-operating parameter joint feature space. This leads to the construction of a Gaussian process surrogate model embedded with physical conservation constraints, providing an efficient and physically interpretable surrogate response surface for the multi-objective manifold Bayesian optimization in step S5. The detailed steps are as follows:

[0113] Step S4-1: Extraction of material map feature descriptors:

[0114] Starting from the target electrode material upon which the reduced-order cross-scale dynamics model constructed in step S3 depends, the embedded representation of the intermediate layer of the physical information graph neural network microdynamics model in step S2 is extracted. Specifically, the crystal structure graph data of the target electrode material is re-input into the node feature embedding layer trained in step S2-1 to obtain the high-dimensional atomic representation vectors of each graph node. (in This represents the total number of message passing layers. , (The total number of atoms in the crystal structure); input the same crystal structure diagram data into the edge feature embedding layer trained in step S2-1 to obtain the high-dimensional chemical bond representation vector of each edge. (in , (the set of all objects with edges).

[0115] Global average pooling is performed on the high-dimensional atomic representation vectors to obtain node-level graph feature vectors: ,in, For node-level graph feature vectors;

[0116] Performing global max pooling on the high-dimensional chemical bond characterization vector yields the edge-level graph feature vector: ,in, This indicates taking the maximum value element by element. These are the feature vectors of the edge-level graph;

[0117] After concatenating the eigenvectors of the node-level graph and the edge-level graph, principal component analysis (PCA) is performed for dimensionality reduction to obtain the material graph feature descriptor: ; ;

[0118] in, PCA transformation matrix (preserving the original value) One principal component, such that the cumulative variance contribution rate is not lower than a preset threshold, for example, 95%. This is the mean vector of the concatenated feature vectors of all candidate materials (see below). For material map feature descriptors, dimensions ;

[0119] When multiple candidate target electrode materials exist (e.g., electrode materials doped with different transition metals or with different crystal planes), the above extraction process is performed on each candidate material separately to obtain the corresponding set of material map feature descriptors. ,in The number of candidate materials;

[0120] Step S4-2: Macroscopic operating parameter sampling and multidimensional sample set generation:

[0121] With electrode potential (Unit: Vvs. RHE) Electrolyte flow rate (Unit: m·s) -1 Initial concentration of reactants (Unit is mol·m) -3 Electrolyte temperature (Unit: K) and electrolyte pH value (Dimensionless) serves as a macroscopic operational parameter, forming a five-dimensional macroscopic operational parameter vector. The preset variation range of each macroscopic operating parameter is: electrode potential (Determined by the thermodynamic feasibility potential window determined in steps S1-3), electrolyte flow rate (For example m·s -1 to m·s -1 ), initial concentration of reactants (e.g., 1 mol·m -3 Up to 100 mol·m -3 Electrolyte temperature (e.g., 273K to 353K), electrolyte pH value (e.g., 3 to 11);

[0122] Within the preset range of five-dimensional operation parameters, the Latin hypercube sampling method is used to generate... Group of macroscopic operating parameter samples Latin hypercube sampling divides the range of values ​​for each dimension into equal parts. In each interval, a value is randomly selected and then randomly permuted and combined to ensure that the sample is uniformly distributed in the parameter space.

[0123] For each candidate target electrode material ( ), and its material map feature descriptor With each set of macroscopic operating parameter samples The reduced-order turbulent-scale dynamic model, constructed in step S3, is run group by group as input conditions. When running the reduced-order turbulent-scale dynamic model, the electrolyte flow rate is... Converting microchannel geometry parameters into fluid shear force As boundary conditions for the mass transfer equation: (in The dynamic viscosity of the electrolyte. (for microchannel height); converting electrolyte pH to bulk proton concentration boundary conditions: mol·m -3 ; Initial concentration of reactants Set the bulk nitrate concentration boundary condition;

[0124] After each set of input conditions reaches steady state, the corresponding steady-state denitrification rate is extracted. Nitrogen deLaday selectivity and nitrogen removal per unit energy consumption As a macroscopic performance label; by summarizing the input-output data of all candidate materials and all operating parameter samples, a multidimensional sample set containing material diagram feature descriptors and macroscopic operating parameters is obtained: ;

[0125] The total number of samples in the multidimensional sample set is ;

[0126] Step S4-3: Construction of high-dimensional spatial proximity graph:

[0127] For multidimensional sample sets Constructing a high-dimensional spatial proximity graph: First, represent each sample point as a high-dimensional joint feature vector: ,in, For the global index of the sample points, ;

[0128] Each sample point in the multidimensional sample set is used as a graph node. Nearest neighbor weighted edges are constructed based on the Euclidean distance between sample points; for sample points and sample points The Euclidean distance between them is: ;

[0129] For each sample point Select the one with the smallest Euclidean distance. Each sample point is considered its nearest neighbor. To preset the number of nearest neighbors, for example ), and at sample points With each of its nearest neighbor sample points Construct a weighted edge between them, with the edge weight being: ,in, For sample points To its first The distance between nearest neighbor sample points (local scale parameter); this yields a high-dimensional sample neighborhood map. ,in for A set of graph nodes, For the nearest neighbor set, This is the edge weight matrix (the weights between non-nearest neighbor sample points are set to zero).

[0130] Step S4-4: Manifold learning to extract intrinsic descriptors:

[0131] Neighborhood map of high-dimensional samples Construct a manifold learning objective function, which includes a local geometry preservation term and a global distance preservation term:

[0132] Local geometry preservation term The linear reconstruction weights of adjacent sample points in the low-dimensional embedding space are consistent with those in the high-dimensional space; for sample points Its linear reconstruction weights in high-dimensional space are obtained by minimizing the following local reconstruction error: ;

[0133] Satisfy normalization constraints Solving this constrained least squares problem yields the optimal reconstructed weight vector. The local geometry preservation term is then defined as: ,in, Let the coordinate matrix of all sample points be in the low-dimensional embedding space. For sample points Coordinate vectors in a low-dimensional embedding space Preset dimensions for intrinsic descriptors (e.g.) or );

[0134] Global distance preservation The geodesic distance between sample points in the low-dimensional embedding space is proportional to the geodesic distance between sample points in the high-dimensional space; firstly, in the high-dimensional sample neighborhood map... Above, Dijkstra's shortest path algorithm is used to calculate the sample between any two points. and Map geodesic distance between (i.e., the sum of the shortest path weights along the graph edges); then the global distance preservation term is defined as: ;

[0135] The complete manifold learning objective function is a weighted combination of local geometry preservation terms and global distance preservation terms: ,in, Regularization weights for the global distance preservation term (e.g.) );

[0136] Iterative optimization algorithms (such as gradient-based Adam optimizers or quasi-Newton-LBFGS algorithms) are used to learn the objective function of the manifold. The solution is performed using random initialization or PCA initialization as the initial values ​​for the low-dimensional embedding coordinates, iteratively optimizing until the change in the objective function value is less than a preset convergence threshold (e.g., ...). ), to obtain the optimal low-dimensional embedding coordinate matrix ; for each sample point The corresponding optimal low-dimensional embedding coordinate vector As its intrinsic catalytic descriptor: ,in, For sample points Intrinsic catalytic descriptor, dimension Much smaller than the dimension of the original high-dimensional joint feature vector This completes the extraction of the intrinsic catalytic descriptor.

[0137] Step S4-5: Construction of the physical constraint Gaussian process surrogate model:

[0138] intrinsic catalytic descriptor and macroscopic operating parameters Concatenate them into a joint input feature vector: ,in, ; steady-state denitrification rate Ammonia / Lafarge Selectivity and nitrogen removal per unit energy consumption Construct as a multi-task output vector: ;

[0139] Using the joint input feature vector as the independent variable and the multi-task output vector as the dependent variable, a multi-output Gaussian process regression framework is constructed. Assuming that the three output tasks share a latent hidden Gaussian process through the linear model common regionization (LMC) method, the prior assumptions for multi-output Gaussian process regression are: ,in, For input The corresponding multi-task output vector, mean function ( (The linear mean weight matrix), covariance matrix Defined by physical constraints;

[0140] Physically constrained composite kernel function Material similarity kernel function Operation parameters kernel function and cross-coupled kernel functions It consists of three parts: material similarity kernel function. Material structural similarity is calculated based on the cosine similarity between intrinsic catalytic descriptors: ,in, For the first The intrinsic catalytic descriptor vector of each sample. For the first The intrinsic catalytic descriptor vector of each sample. The signal variance hyperparameter of the material similarity kernel;

[0141] Operational parameters kernel function Calculate the similarity of operating conditions based on the radial basis functions (RBF) between macroscopic operating parameters: ,in, The signal variance hyperparameter is the kernel of the operating parameters. For the first The macroscopic operating parameter vector of each sample For the first The macroscopic operating parameter vector of each sample An anisotropic length-scale diagonal matrix (the five diagonal elements correspond to the length scales of electrode potential, electrolyte flow rate, initial reactant concentration, electrolyte temperature, and electrolyte pH, respectively).

[0142] Cross-coupled kernel function Multiply the material similarity kernel function and the operating parameter kernel function to capture the synergistic effect of materials and operating conditions: ,in, The signal variance hyperparameter of the cross-coupled core;

[0143] The complete physical constraint composite kernel function is a weighted combination of the three kernel functions mentioned above: ;

[0144] Construct physical conservation constraints; these constraints include thermodynamic consistency constraints and mass conservation constraints.

[0145] Thermodynamic consistency constraints require nitrogen Faraday selectivity Nitrogen removal per unit energy consumption The product is within the pre-defined thermodynamic feasible region: ,in, and These are the lower and upper bounds of the thermodynamic feasible region, respectively, determined by the thermodynamic reference data in step S1-2;

[0146] The mass conservation constraint requires a steady-state denitrification rate With nitrogen Faraday selectivity The product of these factors is not greater than the theoretical maximum denitrification rate corresponding to the initial reactant concentrations. ,in, is the diffusion coefficient of nitrate ions. The thickness of the diffusion layer;

[0147] The physical conservation constraints are transformed into regularization penalty terms; for each training sample point... Define the physical constraint violation penalty function: ,in, and These are the penalty coefficients for the thermodynamic consistency constraint and the mass conservation constraint, respectively.

[0148] By adding the regularization penalty term to the negative log-likelihood function of the multi-output Gaussian process regression framework, we obtain the physical constraint loss function: ,in, The joint input feature matrix for all training samples, The multi-task output matrix for all training samples. The set of all hyperparameters (including the linear mean weight matrix) for a multi-output Gaussian process regression framework. The signal variance of the three kernel functions , , Anisotropic length scale matrix (observation noise variance and LVMC coefficient matrix).

[0149] The bounded memory quasi-Newton optimization algorithm (L-BFGS) is used to apply the physical constraint loss function. Hyperparameter optimization is performed; the L-BFGS algorithm approximates the inverse of the Hessian matrix by maintaining finite-step historical gradient information, and updates the hyperparameter vector along the quasi-Newton direction in each iteration. until the gradient norm of the physical constraint loss function is less than a preset convergence threshold (e.g., ...). Or reach the preset maximum number of iterations (e.g., 500 times); after optimization, obtain the optimal set of hyperparameters. and the trained physical constraint Gaussian process surrogate model;

[0150] The trained physical-constrained Gaussian process surrogate model can handle any given joint input feature vector. Predicted output denitrification rate Ammonia selective prediction value and energy consumption forecast And the corresponding estimate of prediction uncertainty (prediction variance). , and .

[0151] In this embodiment, the core function is to utilize the physically constrained Gaussian process surrogate model and intrinsic catalytic descriptor low-dimensional manifold structure constructed in step S4. Through multi-objective manifold Bayesian optimization, it efficiently searches within the joint design space of intrinsic catalytic descriptor-operating parameters, automatically balancing the three conflicting objectives of maximizing the denitrification rate, maximizing ammonia selectivity, and minimizing energy consumption. Simultaneously, a thermodynamic consistency penalty function is used to ensure the physical feasibility of the optimization process, ultimately outputting the optimal electrode material structure characteristics and macroscopic operating parameter combination. The detailed steps are as follows:

[0152] Step S5-1: Joint Design Space and Multi-Objective Optimization Problem Construction:

[0153] intrinsic catalytic descriptor The range of values ​​and macroscopic operating parameters The pre-defined operating parameter variation range constructs the joint design space: Among them, the value range of the intrinsic catalytic descriptor The intrinsic catalytic descriptors of all training sample points in step S4-4 The definition of the hyperrectangular region determined by the minimum and maximum values ​​in each dimension: ;

[0154] Preset operating parameter variation range for macroscopic operating parameters Consistent with the range defined in step S4-2;

[0155] Predicted denitrification rate Maximize the ammonia selective prediction value Maximize and predict energy consumption values Minimization is defined as three conflicting optimization objectives within a joint design space. Construct a multi-objective optimization problem: Among them, taking a negative value for the third objective transforms the energy consumption minimization problem into a maximization problem to unify the optimization direction;

[0156] Step S5-2: Generation of the initial Pareto front based on the lower confidence bound response surface:

[0157] Based on the predicted mean and prediction uncertainty estimate of the physical constraint segregation process surrogate model output after training in steps S4-5, a lower confidence bound (LCB) response surface for the multi-output Gaussian process is constructed; for any point in the joint design space... The LCB values ​​of the three objective functions are as follows: ; ; ;

[0158] in, To explore and utilize trade-off factors (e.g.) (corresponding to approximately 97.7% confidence level); The LCB response surface penalizes regions with high prediction uncertainty, encouraging optimization algorithms to strike a balance between exploring unknown regions and utilizing known high-performance regions;

[0159] On the lower confidence bound response surface, a non-dominated sorting genetic algorithm (NSGA-II) is used for multi-objective optimization search to generate a uniformly distributed initial Pareto front approximate solution set; the parameters of NSGA-II are set as follows: population size Maximum number of generations Crossover probability Probability of mutation The tournament selection size is 2; NSGA-II preserves elite solutions and maintains the diversity of the Pareto front in each generation of evolution through fast non-dominated sorting and crowding distance calculation; after the evolution terminates, all individuals with a non-dominated level of 1 are extracted from the final population to form the initial Pareto front approximate solution set. ,in The initial number of individuals at the Pareto front;

[0160] Step S5-3: Constructing the Bayesian optimization acquisition function for a multi-objective manifold:

[0161] Approximate solution set with initial Pareto front As a set of reference points, the Expected Hypervolume Improvement (EHVI) of the Pareto front is calculated on the set of reference points; let the current iteration step be... The approximate solution set of the Pareto front at time is The corresponding reference point is (The reference point is usually the worst value of each objective function on the current Pareto front, or a conservative value specified by the user), then any candidate point The expected overvolume improvement at this location is: ,in, For hypervolume function (with reference point) (lower boundary) Let the three objective function vectors be: For a physically constrained Gaussian process surrogate model at point The multivariate Gaussian prediction distribution at the location; the expected hypervolume improvement is numerically calculated using the Monte Carlo sampling method, with the number of sampling times... ;

[0162] At the same time, based on candidate points Intrinsic Catalysis Descriptor For each location in the low-dimensional embedding space, the manifold distance metric is computed; the manifold distance metric is defined as the minimum geodesic distance in the intrinsic catalytic descriptor space between the candidate point and each point in the current Pareto front approximation solution set. ,in, This represents the geodesic distance between two points calculated on the low-dimensional embedding manifold learned in step S4-4; when candidate points When a sample is not in the existing training samples, its geodesic distance is approximated by the following method: First, a high-dimensional sample neighborhood map is used. The graph geodesic distance between existing sample points is used to train a local linear regression model from the intrinsic catalytic descriptor space to the high-dimensional graph geodesic distance. Then, this regression model is used to predict the approximate geodesic distance between candidate points and each existing sample point.

[0163] By multiplying the expected hypervolume improvement at the Pareto front with the manifold distance metric, a multi-objective manifold Bayesian optimization acquisition function is constructed: ,in, Weight exponents for manifold distance metrics (e.g.) This function is used to regulate the degree of emphasis on manifold structure exploration; it also encourages candidate points with high expected hypervolume improvement in the target space and appropriate distance from the known Pareto front in the material manifold space, thereby achieving balanced exploration in the target space and the material manifold space.

[0164] Step S5-4: Candidate sampling point generation and physical feasibility determination:

[0165] Based on acquisition function In the joint design space A maximization search is performed to generate a set of candidate sampling points; specifically, a multi-starting-point local optimization strategy is adopted: in the joint design space, low-difference sequences of Sobol are generated. For each initial search point, the LBFGS-B algorithm (a quasi-Newton algorithm with boundary constraints) is used to perform local maximization. After convergence, a set of candidate sampling points is obtained. ( (Remove duplicate convergence points);

[0166] For each candidate sampling point in the candidate sampling point set Feasibility is determined by applying physical conservation constraints; these constraints include thermodynamic consistency constraints and mass conservation constraints, consistent with the constraints defined in step S4-5.

[0167] Determination of thermodynamic consistency constraints: ,in, This is a Faraday-selective predictor of nitrogen gas. This is the predicted value of nitrogen removal per unit of energy consumption;

[0168] Determination of mass conservation constraints: ;

[0169] Candidate sampling points that violate thermodynamic consistency constraints or mass conservation constraints are marked as invalid solutions. The sampling function values ​​of invalid solutions are then negatively penalized and corrected using a thermodynamic consistency penalty function to obtain the corrected sampling function values. ,in, For example, the penalty intensity coefficient (e.g.) ), The physical constraint violation penalty function defined in step S4-5 is calculated as follows: ;

[0170] Based on the corrected acquisition function value The optimal candidate sampling point with the largest corrected sampling function value is selected from the candidate sampling point set. ;

[0171] Step S5-5: High-fidelity numerical simulation evaluation and dataset update:

[0172] Optimal candidate sampling points The reduced-order, multi-scale kinetic model constructed in step S3 is used for high-fidelity numerical simulation evaluation; due to the intrinsic catalytic descriptor It may not correspond to any existing candidate electrode materials, so it is necessary to first obtain the corresponding material features through the inverse mapping from intrinsic catalytic descriptor to material map feature descriptor; the inverse mapping is achieved in the following way:

[0173] The high-dimensional joint features established using manifold learning in step S4-4 To low-dimensional eigenvalue descriptors The positive mapping relationship is used to construct an inverse regression model; specifically, the eigenvalues ​​of all training samples are used. The input is the material map feature descriptor part of the corresponding high-dimensional joint feature vector. For the output, train a multilayer perceptron (MLP) regression model. The MLP contains two hidden layers (64 neurons per layer, ReLU activation function), and is trained by minimizing the mean squared error; after training, the optimal eigenvalue descriptor is... Inputting an MLP regression model yields a material map feature descriptor that is reconstructed inversely. ;

[0174] Material map feature descriptor reconstructed in reverse With optimal macroscopic operating parameters In combination, the corresponding input conditions are set in the reduced-order span dynamics model in step S3: using the same parameter transformation method as in step S4-2, the input conditions are... Electrolyte flow rate Converted to fluid shear force pH value Convert to volume ion concentration boundary conditions For material map feature descriptors This is then matched with the closest candidate material to the physical information neural network trained in step S2—selecting the material map feature descriptor and... The candidate material with the smallest Euclidean distance between them is used as the approximate material, and the crystal structure diagram data of the candidate material is used to drive the reduced-order span dynamics model.

[0175] The reaction-mass transfer bidirectional coupling simulation was run to steady state to obtain the true steady-state denitrification rate corresponding to the optimal candidate sampling point. True ammonia Faraday selectivity and actual nitrogen removal per unit energy consumption , as a true macroscopic performance label;

[0176] Optimal candidate sampling points and their corresponding real macroscopic performance labels and the material map feature descriptor for reverse reconstruction As a new sample, it is added to the multidimensional sample set in step S4. In the middle, update the high-dimensional sample neighborhood graph. Repeat step S4-4 of the manifold learning algorithm to update the high-dimensional sample neighborhood graph and extract the updated intrinsic descriptors. The finite-memory quasi-Newton optimization algorithm (L-BFGS) in steps S4-5 is used to update the hyperparameters of the physical constraint loss function. The physical constraint Gaussian process surrogate model is then retrained with the complete dataset after the addition of new samples to obtain the updated physical constraint Gaussian process surrogate model.

[0177] Steps S5-6: Convergence Judgment and Iterative Optimization:

[0178] Calculate the Pareto front hypervolume improvement corresponding to the updated physically constrained Gaussian process surrogate model. That is, the difference in hypervolume between the approximate solution set of the Pareto front before and after the update: ;

[0179] Determine if the overvolume improvement is less than the preset convergence threshold. (For example ):

[0180] like If so, return to step S5-3, recalculate the acquisition function based on the updated surrogate model and the updated Pareto front approximate solution set, and perform the next round of iterative optimization;

[0181] like If the optimization has converged, then the iterative optimization should be stopped.

[0182] Step S5-7: Optimal Solution Selection and Output:

[0183] After optimization and convergence, the approximate solution set from the current Pareto front is obtained. In the middle, according to the preset decision preference weight vector (satisfy For example, in applications that prioritize nitrogen removal... ), calculate the weighted composite score for each Pareto front solution: ,in, , and These are the optimal values ​​for each objective on the Pareto front (maximum nitrogen removal rate, maximum selectivity, and minimum energy consumption), used for normalization; a weighted comprehensive score is selected. The highest Pareto front solution is taken as the optimal solution;

[0184] The intrinsic catalytic descriptor corresponding to the optimal solution The optimal material map feature descriptor is reconstructed using the MLP inverse mapping model described in step S5-5. Furthermore, from the set of candidate electrode materials, material map feature descriptors are selected and... The candidate material with the smallest Euclidean distance is selected as the optimal electrode material, and its crystal structure characteristics (including lattice parameters, atomic occupancy, exposed crystal planes, and coordination environment of active sites) represent the optimal electrode material structure characteristics. The macroscopic operating parameters corresponding to the optimal solution are then used. As the optimal combination of macroscopic operating parameters;

[0185] The optimal combination of electrode material structure characteristics and macroscopic operating parameters is output to complete the reaction kinetics modeling and parameter optimization of the electrocatalytic denitrification process.

[0186] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process, characterized in that: Includes the following steps: S1: Obtain basic data of the electrocatalytic denitrification reaction system. The basic data includes crystal structure diagram data of the target electrode material, thermodynamic and kinetic baseline data of nitrogen-containing intermediate species at different electrode potentials, and macroscopic reaction kinetics and mass transfer experimental data. S2: Construct a potential-dependent physical information graph neural network microdynamic model. Input crystal structure diagram data and double-layer electric field physical constraints into the physical information graph neural network to predict the activation energy of elementary reactions under multi-path competition. Combine the transition state theory to generate a set of potential-dependent microdynamic rate equations. The multi-path competition includes the nitrogen denitrification to nitrogen generation path, the ammonia reduction to nitrogen generation path, and the hydrogen evolution side reaction path. S3: Construct a reaction-mass transfer bidirectional coupled cross-scale kinetic model, and couple the microscopic kinetic rate equations as source terms into the one-dimensional Nernst-Planck-Poisson mass transfer equation to simulate the local pH dynamics and concentration polarization effect within the electric double layer; calculate the global sensitivity index of each elementary reaction to the macroscopic denitrification index based on adjoint sensitivity analysis, and reduce the order of insensitive elementary reactions to obtain a reduced-order cross-scale kinetic model. S4: Construct a physically constrained manifold surrogate model. Based on the reduced-order cross-scale dynamic model, generate a multi-dimensional sample set containing material diagram feature descriptors and macroscopic operating parameters. Use the manifold learning algorithm to extract intrinsic catalytic descriptors. Construct a physically constrained Gaussian process surrogate model with intrinsic catalytic descriptors and macroscopic operating parameters as inputs and denitrification rate, nitrogen selectivity and energy consumption as outputs. S5: Perform multi-objective manifold Bayesian optimization, and use the acquisition function based on the Pareto front expectation hypervolume improvement to iteratively optimize the physical constraint Gaussian process surrogate model. In the optimization process, a thermodynamic consistency penalty function is introduced to eliminate invalid solutions that violate physical conservation laws, and the optimal electrode material structure characteristics and macroscopic operating parameter combination are output.

2. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 1, characterized in that: Acquire fundamental data for the electrocatalytic denitrification reaction system. This fundamental data includes crystal structure diagrams of the target electrode material, baseline thermodynamic and kinetic data of nitrogen-containing intermediates at different electrode potentials, and macroscopic reaction kinetics and mass transfer experimental data; specifically including: An initial crystal model of the target electrode material is constructed, and the geometric structure of the initial crystal model is optimized using density functional theory. Atomic coordinates and bonding information are extracted from the optimized crystal model. Atoms are defined as graph nodes, chemical bonds are defined as graph edges, and node atomic number and electronegativity characteristics are assigned, while edge bond length and bond order characteristics are assigned, generating crystal structure diagram data of the target electrode material. Based on the crystal structure data of the target electrode material, a solid-liquid interface double-layer calculation model including an implicit solvation model and an explicit water molecule layer is constructed. Different electrode potential conditions are introduced into the solid-liquid interface double-layer calculation model to calculate the adsorption free energy of nitrogen-containing intermediate species on the surface of the target electrode material, obtaining thermodynamic baseline data of nitrogen-containing intermediate species under different electrode potentials. At the same time, the micro-motion elastic band method is used to search for the lowest energy path for nitrogen-containing intermediate species to transform into the transition state, and the activation energy barrier of each elementary reaction is calculated. Combined with the offset correction of the reaction free energy diagram by the electrode potential, kinetic baseline data of nitrogen-containing intermediate species under different electrode potentials are obtained. The thermodynamic baseline data and the kinetic baseline data are integrated into thermodynamic and kinetic baseline data of nitrogen-containing intermediate species under different electrode potentials. Based on the thermodynamic and kinetic baseline data of nitrogen-containing intermediates at different electrode potentials, the thermodynamic feasibility potential window and kinetic rate-determining step of the electrocatalytic denitrification reaction were determined, which served as the basis for setting the operating parameters of the high-throughput microfluidic electrochemical testing array. Under the set operating parameters and fluid shear force conditions, the dynamic response curves of product concentration versus time during the electrocatalytic denitrification reaction were continuously acquired. The dynamic response curves were analyzed in situ using online mass spectrometry and ion chromatography to extract the macroscopic apparent reaction rate constant, limiting current density, and diffusion layer thickness, and the macroscopic reaction kinetics and mass transfer experimental data were compiled.

3. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 1, characterized in that: The construction of a potential-dependent physical information graph neural network microdynamic model specifically includes: Based on the crystal structure diagram data of the target electrode material, a physical information graph neural network skeleton is constructed, which includes a node feature embedding layer, an edge feature embedding layer, and a multi-physics message passing layer. The node feature embedding layer maps the node atomic number and electronegativity features into high-dimensional atomic representation vectors, and the edge feature embedding layer maps the edge bond length and bond order features into high-dimensional chemical bond representation vectors. A double-layer electric field physical constraint module is constructed in the multiphysics message passing layer. The double-layer electric field physical constraint module calculates the local electric field intensity distribution and potential decay curve on the surface of the target electrode material based on the Gouy-Chapman-Stern double-layer theory. The local electric field intensity distribution is superimposed on the high-dimensional atomic characterization vector as the external field perturbation vector of each graph node. The potential decay curve is used as the weight correction factor of each graph edge to correct the high-dimensional chemical bond characterization vector, generating enhanced graph node features and enhanced graph edge features containing electric field physical constraints. A potential-dependent feature injection mechanism is constructed to encode different electrode potential conditions into potential embedding vectors. The potential embedding vectors are adaptively fused with enhanced graph node features containing electric field physical constraints through a gated attention mechanism to obtain a potential-dependent atomic-level reactive activity characterization. The potential-dependent atomic-level reactivity characterization is input into a multi-task prediction head, which includes a reaction path classification subnetwork and an activation energy regression subnetwork. The multi-task prediction head outputs the occurrence probability distribution of each reaction path and the corresponding predicted value of the activation energy of the elementary reaction, thus completing the construction of a potential-dependent physical information graph neural network microdynamic model.

4. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 3, characterized in that: Crystal structure data and double-layer electric field physical constraints are input into a physical information graph neural network to predict the activation energy of elementary reactions under multi-path competition. This is combined with transition state theory to generate a potential-dependent set of microscopic kinetic rate equations. The multi-path competition includes the nitrogen denitrification pathway, the ammonia reduction pathway, and the hydrogen evolution side reaction pathway; specifically: Based on the potential-dependent physical information graph neural network microdynamic model, the crystal structure diagram data of the target electrode material and the physical constraints of the double electric layer electric field are simultaneously input. The reaction path classification sub-network in the multi-task prediction head identifies the elementary reaction sequences of the denitrification to nitrogen generation path, the reduction to ammonia generation path, and the hydrogen evolution side reaction path, respectively. The activation energy regression sub-network predicts the forward activation energy and reverse activation energy of each elementary reaction sequence at the corresponding electrode potential. Based on the forward and reverse activation energies, and combined with the adsorption free energy of nitrogen-containing intermediate species on the surface of the target electrode material, thermodynamic cycle verification conditions for each elementary reaction are constructed. The predicted forward and reverse activation energies are then corrected for thermodynamic consistency using these thermodynamic cycle verification conditions, resulting in the corrected forward and reverse activation energies. Based on the corrected forward activation energy and the corrected reverse activation energy, the microscopic reaction rate constants of each elementary reaction are calculated using transition state theory. The microscopic reaction rate constants are functions of electrode potential, local electric field strength, and temperature. Based on the microscopic reaction rate constants of each elementary reaction, microscopic kinetic rate equations were constructed for the denitrification to nitrogen generation pathway, the reduction to ammonia pathway, and the hydrogen evolution side reaction pathway, respectively. The microscopic kinetic rate equations include the surface coverage change rate of each nitrogen-containing intermediate species, the reactant consumption rate, and the product formation rate. The micro-kinetic rate equations for the denitrification to nitrogen generation pathway, the reduction to ammonia pathway, and the hydrogen evolution side reaction pathway are coupled to construct a potential-dependent set of micro-kinetic rate equations under multi-path competition.

5. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 1, characterized in that: A bidirectional, reaction-mass transfer coupled cross-scale kinetic model was constructed, with the microscopic kinetic rate equations coupled as source terms to the one-dimensional Nernst-Planck-Poisson mass transfer equation to simulate the dynamic evolution of local pH and concentration polarization effects within the electric double layer; specifically including: Based on the reactant consumption rate and product formation rate of each elementary reaction in the microdynamic rate equation set, an interfacial reaction source term vector containing nitrogen-containing intermediate species, protons, hydroxide ions and electrolyte ions is constructed. A one-dimensional spatial discretization grid is constructed along the normal direction perpendicular to the electrode surface. The diffusion layer is divided into multiple grid layers. Concentration field variables and potential field variables are defined in each grid layer to obtain a one-dimensional spatial discretization mass transfer calculation domain. The interface reaction source term vector is coupled as the boundary flux source term to the Nernst-Planck mass transfer equation of the boundary grid layer near the electrode surface in the one-dimensional spatial discretized mass transfer computation domain. The charge transfer rate involving electron transfer steps in each elementary reaction is coupled as the source term to the Poisson potential equation in the one-dimensional spatial discretized mass transfer computation domain, thus obtaining the initial cross-scale dynamic model of reaction-mass transfer bidirectional coupling. The initial multi-scale dynamic model of reaction-mass transfer bidirectional coupling is solved by time stepping using the implicit time integration method. In each time step, the concentration field variables and electric potential field variables are updated alternately based on the Newton-Raphson iteration method. When the residual norm of two adjacent iterations is less than the preset convergence threshold, it is determined to be converged, and the converged concentration field distribution and electric potential field distribution are obtained. The dynamic evolution curve of local pH value in the double layer over time is calculated based on the proton concentration of each grid layer in the converged concentration field distribution. The concentration polarization overpotential is calculated based on the concentration gradient of nitrogen-containing intermediate species in each grid layer in the converged concentration field distribution. The dynamic evolution curve and concentration polarization overpotential are summarized to complete the construction of a cross-scale kinetic model of bidirectional coupling between reaction and mass transfer.

6. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 5, characterized in that: Based on the adjoint sensitivity analysis, the global sensitivity index of each elementary reaction to macroscopic denitrification indicators was calculated. Insensitive elementary reactions were then reduced in order to obtain a reduced-order cross-scale kinetic model; specifically including: Steady-state denitrification rate, nitrogen Faraday selectivity, and denitrification per unit energy consumption are extracted from the reaction-mass transfer bidirectional coupling cross-scale kinetic model as macroscopic denitrification indicators, and the adjoint equation is constructed with the microscopic reaction rate constants of each elementary reaction as perturbation variables. Solving the adjoint equation yields the distribution field of the adjoint variables corresponding to the microscopic reaction rate constants of each elementary reaction. The inner product operation is performed between the distribution field of the adjoint variables and the reaction flux of each elementary reaction in the microscopic kinetic rate equation set to obtain the local sensitivity coefficient of each elementary reaction to the macroscopic denitrification index. Within the preset range of electrode potential variation and temperature variation, the activation energy and pre-exponential factor of each elementary reaction are randomly sampled using the Latin hypercube sampling method to generate multiple sets of kinetic parameter samples. For each set of kinetic parameter samples, a reaction-mass transfer bidirectional coupling cross-scale kinetic model is run to obtain the corresponding macroscopic denitrification index sample set. The local sensitivity coefficient is weighted and integrated with the macroscopic denitrification index sample set to calculate the global sensitivity index of each elementary reaction to the steady-state denitrification rate, nitrogen Faraday selectivity and denitrification per unit energy consumption. Elementary reactions with a global sensitivity index below a preset sensitivity threshold are marked as insensitive elementary reactions. Quasi-steady-state approximations are applied to insensitive elementary reactions. The surface coverage change rate of nitrogen-containing intermediates corresponding to insensitive elementary reactions is set to zero. The surface coverage differential equations corresponding to insensitive elementary reactions are eliminated. The eliminated microscopic kinetic rate equations are recoupled to the one-dimensional Nernst-Planck-Poisson mass transfer equation to obtain a reduced-order cross-scale kinetic model.

7. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 1, characterized in that: A physically constrained manifold surrogate model is constructed, and a multidimensional sample set containing material diagram feature descriptors and macroscopic operating parameters is generated based on a reduced-order multi-scale dynamics model. Intrinsic catalytic descriptors are extracted using a manifold learning algorithm. Specifically, this includes: Crystal structure data of the target electrode material is extracted from the reduced-order multi-scale dynamic model. The high-dimensional atomic representation vectors output by the node feature embedding layer and the high-dimensional chemical bond representation vectors output by the edge feature embedding layer in the physical information graph neural network microdynamic model are then processed. Global average pooling is performed on the high-dimensional atomic representation vectors to obtain node-level graph feature vectors, and global max pooling is performed on the high-dimensional chemical bond representation vectors to obtain edge-level graph feature vectors. The node-level graph feature vectors and edge-level graph feature vectors are then concatenated and subjected to principal component analysis for dimensionality reduction to obtain the material graph feature descriptor. Using electrode potential, electrolyte flow rate, initial reactant concentration, electrolyte temperature and electrolyte pH as macroscopic operating parameters, multiple sets of macroscopic operating parameter samples are generated using the Latin hypercube sampling method within the preset range of operating parameter variations. The material graph feature descriptor is combined with each set of macroscopic operating parameter samples as input conditions. The reduced-order cross-scale dynamic model is run group by group. The steady-state denitrification rate, nitrogen Faraday selectivity and denitrification per unit energy consumption corresponding to each set of input conditions are extracted as macroscopic performance labels. The material graph feature descriptor, macroscopic operating parameter samples and macroscopic performance labels are summarized to obtain a multidimensional sample set containing the material graph feature descriptor and macroscopic operating parameters. A high-dimensional spatial proximity graph is constructed for the multidimensional sample set. Each sample point in the multidimensional sample set is used as a graph node, and nearest neighbor weighted edges are constructed based on the Euclidean distance between sample points to obtain the high-dimensional sample proximity graph. A manifold learning objective function is constructed on the high-dimensional sample proximity graph. The manifold learning objective function includes a local geometry preservation term and a global distance preservation term. The local geometry preservation term constrains the linear reconstruction weights of adjacent sample points in the low-dimensional embedding space to be consistent with those in the high-dimensional space. The global distance preservation term constrains the geodesic distance between sample points in the low-dimensional embedding space to be proportional to the geodesic distance between sample points in the high-dimensional space. An iterative optimization algorithm is used to solve the manifold learning objective function to obtain the coordinate vector of each sample point in the low-dimensional embedding space in the multidimensional sample set. The coordinate vector is then used as the intrinsic catalytic descriptor to complete the extraction of the intrinsic catalytic descriptor.

8. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 7, characterized in that: A physically constrained Gaussian process surrogate model is constructed, taking intrinsic catalytic descriptors and macroscopic operating parameters as inputs and denitrification rate, nitrogen selectivity, and energy consumption as outputs; specifically including: The intrinsic catalytic descriptor and macroscopic operating parameters are concatenated into a joint input feature vector. The steady-state denitrification rate, nitrogen Faraday selectivity and denitrification per unit energy consumption are constructed into a multi-task output vector. A multi-output Gaussian process regression framework is constructed with the joint input feature vector as the independent variable and the multi-task output vector as the dependent variable. In the multi-output Gaussian process regression framework, a physically constrained composite kernel function is constructed. The physically constrained composite kernel function includes a material similarity kernel function, an operating parameter kernel function, and a cross-coupling kernel function. The material similarity kernel function calculates the material structure similarity based on the cosine similarity between intrinsic catalytic descriptors. The operating parameter kernel function calculates the operating condition similarity based on the radial basis function between macroscopic operating parameters. The cross-coupling kernel function performs a product operation on the material similarity kernel function and the operating parameter kernel function to capture the synergistic effect between materials and operating conditions. Physical conservation constraints are constructed, which include thermodynamic consistency constraints and mass conservation constraints. The thermodynamic consistency constraints require that the product of nitrogen Faraday selectivity and nitrogen removal per unit energy consumption be within the preset thermodynamic feasible region. The mass conservation constraints require that the product of steady-state denitrification rate and nitrogen Faraday selectivity be no greater than the theoretical maximum denitrification rate corresponding to the initial concentration of reactants. The physical conservation constraints are transformed into regularization penalty terms, which are then superimposed onto the negative log marginal likelihood function of the multi-output Gaussian process regression framework to obtain the physical constraint loss function. A finite-memory quasi-Newton optimization algorithm is used to optimize the hyperparameters of the physical constraint loss function, resulting in a trained physical constraint Gaussian process surrogate model. The trained physical constraint Gaussian process surrogate model can output the predicted values ​​of denitrification rate, nitrogen selectivity, and energy consumption, as well as the corresponding prediction uncertainty estimates.

9. The method for reaction kinetic modeling and parameter optimization of an electrocatalytic denitrification process according to claim 1, characterized in that: Multi-objective manifold Bayesian optimization is performed, using a data acquisition function based on the Pareto front-based expectation hypervolume improvement to iteratively optimize a physically constrained Gaussian process surrogate model. During the optimization process, a thermodynamic consistency penalty function is introduced to eliminate invalid solutions that violate physical conservation laws, outputting the optimal electrode material structural characteristics and macroscopic operating parameter combination; specifically including: A joint design space is constructed using the value range of the intrinsic catalytic descriptor and the preset operating parameter variation range of the macroscopic operating parameters. The three conflicting optimization objectives of maximizing the predicted value of denitrification rate, maximizing the predicted value of nitrogen selectivity, and minimizing the predicted value of energy consumption are set as three mutually conflicting optimization objectives. A multi-objective optimization problem is constructed in the joint design space. Based on the prediction uncertainty estimation of the output of the trained physical constraint Gaussian process surrogate model, a lower confidence bound response surface of the multi-output Gaussian process is constructed. On the lower confidence bound response surface, a non-dominated sorting genetic algorithm is used to generate a uniformly distributed initial Pareto front approximate solution set. Using the initial Pareto front approximation solution set as the reference point set, the expected hypervolume improvement of the Pareto front is calculated on the reference point set. The expected hypervolume improvement of the Pareto front is multiplied by the manifold distance metric calculated based on the geodesic distance of the intrinsic catalytic descriptor in the low-dimensional embedding space to construct a multi-objective manifold Bayesian optimization acquisition function. Based on the acquisition function, a set of candidate sampling points is obtained by maximizing the search in the joint design space. Physical conservation constraints are applied to each candidate sampling point in the set to determine its feasibility. The physical conservation constraints include thermodynamic consistency constraints and mass conservation constraints. Candidate sampling points that violate thermodynamic consistency constraints or mass conservation constraints are marked as invalid solutions. The acquisition function values ​​of invalid solutions are negatively penalized and corrected using a thermodynamic consistency penalty function to obtain the corrected acquisition function values. Based on the corrected acquisition function value, the optimal candidate sampling point with the largest acquisition function value is selected from the candidate sampling point set. The optimal candidate sampling point is input into the reduced-order multi-scale dynamic model for high-fidelity numerical simulation evaluation. The true steady-state denitrification rate, true nitrogen Faraday selectivity, and true denitrification amount per unit energy consumption corresponding to the optimal candidate sampling point are obtained as true macroscopic performance labels. The optimal candidate sampling point and its corresponding true macroscopic performance labels are added as new samples to the multidimensional sample set. The manifold learning algorithm is re-executed to update the high-dimensional sample proximity graph and extract the updated intrinsic catalytic descriptor. The finite-memory quasi-Newton optimization algorithm is used to update the hyperparameters of the physical constraint loss function to obtain the updated physical constraint Gaussian process surrogate model. The algorithm determines whether the Pareto front hypervolume improvement corresponding to the updated physical constraint Gaussian process surrogate model is less than a preset convergence threshold. If it is greater than the preset convergence threshold, it returns to the step of maximizing the search in the joint design space based on the acquisition function for the next round of iterative optimization. If it is less than the preset convergence threshold, it stops iterative optimization. The optimal solution is selected from the current Pareto front approximate solution set according to the preset decision preference weight. The intrinsic catalytic descriptor corresponding to the optimal solution is reconstructed into the crystal structure features of the target electrode material through inverse mapping. The macroscopic operating parameters corresponding to the optimal solution are used as the optimal macroscopic operating parameter combination, and the optimal electrode material structure features and macroscopic operating parameter combination are output.

Images