Coordination control method of photovoltaic off-grid hydrogen production system

Abstract

The invention provides a coordination control method for a photovoltaic off-grid hydrogen production system, which comprises the following steps of: acquiring voltage and current time domain signals of a photovoltaic array, constructing an impedance frequency characteristic mapping model through tensor decomposition and a radial basis function network, and calculating converter control parameters; based on converter control parameters and environment data, temperature-sensing-free self-correction photovoltaic characteristic prediction is achieved through electrical parameter indirect temperature estimation, and an MPPT reference track and a power prediction track are generated; multi-time-scale data coordination processing of millisecond-level electrical control, second-level pressure regulation and minute-level temperature optimization is executed, potential conflicts are solved, and an execution instruction set is generated; and generating and issuing a system comprehensive control instruction by adopting a double-layer deep reinforcement learning network. According to the method, the system efficiency and stability are improved, and the performance is excellent especially under the illumination fluctuation working condition.

Description

Technical Field

[0001] The present invention relates to off-grid photovoltaic hydrogen production, and in particular to a coordinated control method for an off-grid photovoltaic hydrogen production system. Background Art

[0002] Solar photovoltaic hydrogen production, as a zero-carbon energy conversion method, enables long-term energy storage and holds significant potential for application in remote areas, microgrids, and energy interconnected networks. Off-grid photovoltaic hydrogen production systems convert volatile solar energy into stable, storable hydrogen energy through water electrolysis, effectively addressing the intermittent and random nature of photovoltaic power generation and representing a key technological path for achieving efficient renewable energy utilization.

[0003] Currently, control methods for photovoltaic hydrogen production systems primarily focus on single-domain, single-timescale control strategies, such as traditional MPPT control based on the perturbation and observation (P&O) method and the incremental conductance method (INC), electrolyzer temperature control based on PID, and pressure regulation technology based on model predictive control (MPC). These methods typically operate independently and lack system-level coordination. In photovoltaic modeling, one- or two-dimensional lookup tables are primarily used to establish a mapping relationship between illumination, temperature, and output power, or equivalent circuit models are used to describe the electrical characteristics of photovoltaic arrays. For the study of system dynamic characteristics, small-signal analysis methods are often used to construct transfer function models or state-space models, but these models are generally only applicable to specific operating conditions and have limited generalization capabilities.

[0004] However, the core technical bottleneck faced by existing photovoltaic off-grid hydrogen production systems is the problem of coordinated control in multiple physical domains and multiple time scales. Specifically, the dynamic coupling mechanism between photovoltaic electrical characteristics and electrolytic chemical processes in existing systems has not been effectively analyzed, resulting in a prominent power mismatch between the photovoltaic array and the electrolytic load under conditions of rapid light fluctuations, and a significant decrease in hydrogen production efficiency. At the same time, the traditional control strategy lacks information interaction between the subsystems and is unable to coordinate and deal with conflicts between millisecond-level electrical response, second-level pressure regulation, and minute-level temperature optimization. Control command contradictions are very likely to occur in the dynamic process of the system, seriously affecting system stability and energy efficiency. These key issues restrict the promotion and application of photovoltaic off-grid hydrogen production systems in actual application scenarios, especially in areas with poor resource conditions. Summary of the Invention

[0005] The purpose of the invention is to provide a coordinated control method for a photovoltaic off-grid hydrogen production system, in order to solve one of the above-mentioned problems existing in the prior art.

[0006] The technical solution is a coordinated control method for a photovoltaic off-grid hydrogen production system, including: Collect the time domain signals of PV array voltage and current, build impedance frequency characteristic mapping model, and calculate converter control parameters; Collect ambient temperature data and light intensity data, perform temperature sensor-free self-correction photovoltaic characteristic prediction, generate MPPT reference trajectory and short-term power prediction trajectory, and combine them with converter control parameters to perform multi-time scale data coordination and processing to generate the final execution instruction set; Based on the final execution instruction set, the constructed multi-domain coupling model and two-layer deep reinforcement learning network are used to generate and issue system comprehensive control instructions.

[0007] Beneficial effects: The present invention reduces the system's response time to changes in impedance characteristics, effectively avoids system instability caused by impedance mismatch, and improves the power output of the photovoltaic array; effectively solves the problem of large differences in dynamic characteristics of different physical processes, improves system efficiency and control stability, and performs particularly well under conditions of fluctuating light. BRIEF DESCRIPTION OF THE DRAWINGS

[0008] Figure 1 This is a flowchart of the steps of a coordinated control method for a photovoltaic off-grid hydrogen production system provided in an embodiment of the present application.

[0009] Figure 2 This is a flowchart of the steps for calculating converter control parameters provided in an embodiment of the present application.

[0010] Figure 3 This is a flowchart of the steps for forming a global mapping model provided in an embodiment of the present application.

[0011] Figure 4 This is a flowchart of the steps provided in an embodiment of the present application for obtaining the current impedance parameters and calculating the converter control parameters based on them. DETAILED DESCRIPTION

[0012] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0013] It should be noted that to clearly illustrate the steps of this application, serial numbers are assigned to each step in the specification. These serial numbers are for illustrative purposes only and do not limit the order in which the steps must be executed. In actual operation, depending on the technical requirements of the specific implementation scenario, the steps may be executed in a different order than shown in the specification, and in some cases, parallel processing between steps may be implemented.

[0014] The study found that there are many technical problems in the existing control methods of photovoltaic off-grid hydrogen production systems: First, the traditional impedance characteristic characterization method relies on complex equivalent circuit models, which has high computational overhead and limited accuracy, and it is difficult to capture the dynamic impedance changes under different light and temperature conditions; second, the existing MPPT technology is heavily dependent on temperature sensors. Sensor failure causes a sharp drop in system efficiency and cannot solve the problem of distributed temperature unevenness; third, the dynamic characteristics of different physical processes vary significantly (electrical milliseconds, pressure seconds, temperature minutes), and the independent control of each subsystem leads to poor coordination, which is prone to system oscillation or even loss of control when the operating conditions suddenly change; fourth, a single control algorithm is difficult to balance the conflicting needs of system efficiency and stability. Especially under conditions of rapid light fluctuations, the existing control strategy cannot achieve stable output of hydrogen production rate, which seriously restricts the practical application of off-grid hydrogen energy systems.

[0015] like Figure 1 As shown, a coordinated control method for a photovoltaic off-grid hydrogen production system is provided, comprising: S1. Collect the time domain signals of the photovoltaic array voltage and current, build the impedance frequency characteristic mapping model, and calculate the converter control parameters; Specifically, the system reads the time-domain voltage and current signals of the photovoltaic array, collects impedance characteristic data sets under various lighting and temperature conditions, and constructs a database of raw impedance spectrum. The high-dimensional data is then subjected to tensor decomposition, and a radial basis function network is used to construct a four-dimensional mapping relationship, forming a global mapping model. A compressed mapping model is constructed using model sparsification and quantization techniques, and an indexed mapping model is created by designing a quadtree or octree index structure. This enables efficient retrieval of current impedance parameters and calculates the optimal converter control parameters for subsequent electrical control.

[0016] S2, collects ambient temperature data and light intensity data, performs temperature-sensor-free self-correction photovoltaic characteristic prediction, and generates MPPT reference trajectory and short-term power prediction trajectory; Specifically, ambient temperature data, light intensity data, and photovoltaic output power are collected to establish a thermal balance model for the photovoltaic modules. Model parameters are identified using an adaptive recursive least squares algorithm, and a set of calibration parameters is output. Partial array sampling is used to obtain equivalent open-circuit voltage and equivalent short-circuit current. A fused temperature estimate is calculated using parameter credibility assessment metrics. The complete PV curve characteristics are predicted based on the estimated temperature and light intensity. An error correction mapping function is applied to correct the curve, generating the corrected maximum power point parameters, maximum power point tracking (MPPT) reference trajectory, and short-term power prediction trajectory, which serve as inputs for multi-timescale control.

[0017] S3, combining the MPPT reference trajectory, short-term power prediction trajectory and converter control parameters, performs multi-time scale data coordination and processing, and generates the final execution instruction set; Specifically, at the millisecond-level electrical control layer, high-frequency electrical signals are collected and combined with the MPPT reference trajectory and converter control parameters to generate electrical fast-response control instructions and event characteristics. At the second-level pressure control layer, based on event characteristics and short-term power forecast trajectories, the optimal pressure set point is calculated using the transformer regulation prediction model, and pressure regulation control instructions and regulation status are generated. At the minute-level temperature control layer, external weather forecasts and system historical data are combined to generate energy forecast data, calculate the long-term optimal temperature set point, and output temperature optimization control instructions and temperature control status data. A cross-timescale data coordination algorithm is used to process control instructions across three timescales, resolve potential conflicts, and generate a final execution instruction set for system execution and reinforcement learning training.

[0018] S4. Based on the final execution instruction set, the constructed multi-domain coupling model and two-layer deep reinforcement learning network are used to generate and issue system comprehensive control instructions.

[0019] Specifically, an impedance-frequency characteristic mapping model, a combined thermoelectric model, and an electrolysis temperature-pressure characteristic model are integrated to construct a photovoltaic-electrolysis multi-domain coupling model and encapsulate it as a reinforcement learning environment interface. The global state space and action space are defined, a hierarchical reward function is constructed, and a knowledge-enhanced policy network is trained through the reinforcement learning environment interface in combination with a posterior sampling exploration strategy. Based on the decision output by the global policy, the local state space and action space are designed, and an adaptive execution network is trained in combination with a gradient-guided exploration strategy. The current system state is input into a two-layer network to generate reinforcement learning control instructions, which are then integrated with the final execution instruction set into a comprehensive system control instruction. These instructions are distributed to each execution unit to achieve system coordinated control, and execution feedback data is collected and stored in the system operation database for continuous optimization.

[0020] like Figure 2 As shown, according to one aspect of the present application, calculating the converter control parameters includes: S11, collecting the time domain signals of the photovoltaic array voltage and current under various light intensity and temperature conditions, and processing them through fast Fourier transform (FFT) to obtain the impedance modulus and phase angle to form the original impedance spectrum; S12, performing tensor decomposition on the original impedance spectrum, and based on the tensor decomposition results, using a radial basis function network to construct a four-dimensional mapping relationship of light-temperature-frequency-impedance to form a global mapping model; S13, performing sparseness and quantization on the global mapping model, configuring an index structure, and obtaining an indexed mapping model; S14. Input the current operating condition parameters into the index mapping model to obtain the current impedance parameters; and calculate the converter control parameters based on the current impedance parameters.

[0021] Specifically, the time-domain voltage and current signals of the photovoltaic array are read and processed via fast Fourier transform to obtain frequency-domain impedance data. Under varying light intensities S and temperatures T, the frequency ω is swept to acquire an impedance characteristic dataset {ZPV(jω, S, T)}. For each operating point, the impedance modulus and phase angle are extracted to construct a raw impedance spectrum database. A tensor decomposition algorithm is applied to the raw impedance spectrum database, decomposing the four-dimensional data (jω, S, T, ZPV) into a combination of a low-dimensional core tensor and a factor matrix, reducing the number of parameters. In this decomposed low-dimensional space, a local radial basis function (RBF) network is used to construct the mapping: ZPV(jω, S, T) = Σ wi•φi(||X-Xi||), where X = [jω, S, T] is the input feature vector, wi is the weight coefficient, φi is the radial basis function, and Xi is the center point. An adaptive cross-validation method is used to determine the optimal number and distribution of RBFs, balancing model complexity and accuracy. The trained RBF network model is saved as a global mapping model to implement a four-dimensional mapping relationship between light, temperature, frequency, and impedance. Sparsification and quantization techniques are applied to the global mapping model, retaining key RBF nodes and constructing a compressed mapping model, reducing data volume by over 80%. A quadtree or octree index structure is designed to organize the compressed mapping data into an efficient query format and generate an indexed mapping model. Based on the current operating parameters (ωcurrent, Scurrent, Tcurrent), the nearest neighbor region is retrieved from the indexed mapping model. The current impedance parameter ZPV,current is calculated through fast interpolation and substituted into the converter control model to calculate the optimal converter control parameters for subsequent electrical fast response control.

[0022] like Figure 3 As shown, according to one aspect of the present application, a global mapping model is formed, including: The original impedance spectrum is organized into a fifth-order tensor, and Tucker decomposition is applied to obtain the core tensor and factor matrix, generating a decomposed data tensor, and extracting the characteristic distribution spectrum from the decomposed data tensor; Based on the feature distribution map, a three-layer radial basis function network structure is constructed and the initial RBF center point is determined to form an initial RBF network structure model and train it. The number of RBF nodes is dynamically adjusted during training, and a global mapping model is obtained through Bayesian hyperparameter optimization.

[0023] Specifically, the high-dimensional data tensor decomposition preprocessing is performed. The original impedance spectrum database is read and the four-dimensional data {jω, S, T, ZPV} is organized into a fifth-order tensor X∈R (N_ω×N_S×N_T×2), where N_ω, N_S, and N_T represent the number of sampling points for frequency, light intensity, and temperature, respectively, and 2 represents the modulus and phase of the impedance. The Tucker decomposition algorithm is applied to decompose the tensor X into a core tensor and a factor matrix. An adaptive rank selection algorithm is used to determine the optimal rank, balancing data compression rate and reconstruction error, to generate a decomposed data tensor. Statistical analysis is performed on the decomposed data tensor, calculating the importance weights of each element in the core tensor and extracting feature distribution maps, which serve as a reference for subsequent RBF network design.

[0024] Construct a radial basis function network structure. Read the decomposed data tensor and feature distribution map, and determine the complexity requirement of the RBF network based on the feature distribution. Design a three-layer radial basis function network structure: input layer, hidden layer, and output layer. Based on the feature distribution map, use the adaptive K-means++ algorithm to determine the initial RBF center point {c_i} so that its distribution matches the data density: c_i = argmax_x{min_{j <i} ||x-c_j|| 2 • prob(x)}, where prob(x) is the sampling probability determined by the feature distribution map. Initialize the RBF width {σ_i}. Construct the initial RBF network structure model, including the center point {c_i}, width {σ_i}, and initial weights {w_i}.

[0025] Adaptive mapping model training and validation. Read the initial RBF network structure model and training dataset (80% random samples from the original impedance spectrum database), set the learning rate η, regularization coefficient λ, and maximum number of iterations N_iter. Use weighted least squares to train the output layer weights {w_i}, solving the optimization problem: min_{w} Σ||y_k - Σw_i•φ_i(x_k)|| 2 + λ•||w|| 2 , where φ_i(x)=exp(-||x-c_i|| 2 / (2σ_i 2 )) is the activation function of the i-th RBF node, w is the weight vector; y_k is the target output value (impedance modulus and phase); x_k is the input feature vector (frequency, light, temperature), ||w|| 2 is the L2 norm of the weight. Implement an adaptive growing-pruning algorithm to dynamically adjust the number of RBF nodes: calculate the prediction error of each sample. If the regional error exceeds the threshold ε_add, add a new RBF node; evaluate the contribution of each RBF node. If it is lower than the threshold ε_prune, delete the node; re-optimize the parameters {c_i, σ_i, w_i} of the retained nodes; evaluate the model performance on the validation set (20% of the samples from the original impedance spectrum database), and calculate RMSE, MAE, and R 2The performance evaluation report is generated based on the metrics. Bayesian hyperparameter optimization is applied to adjust the regularization coefficient λ and RBF parameters to minimize the validation error and obtain the final global mapping model. The global mapping model is saved to the system model library for subsequent model compression and real-time query.

[0026] like Figure 4 As shown, according to one aspect of the present application, obtaining the current impedance parameter and calculating the converter control parameter accordingly includes: Calculate the contribution of each RBF node in the global mapping model and arrange them in descending order. Filter the nodes whose contribution ranking exceeds the preset threshold as a sparse node set and quantize them at fixed points to obtain a quantization parameter set. Combine the sparse node set and the quantization parameter set to form a compressed mapping model. Based on the RBF center point and coverage radius in the compression mapping model, a tree index structure is constructed, the RBF node contribution is calculated for the leaf nodes to generate a node influence table, and a spatial hash index is constructed to form an indexed mapping model. The indexed mapping model is used to locate leaf nodes according to the current operating parameters, and the most relevant RBF nodes are extracted to calculate the current impedance parameters and generate the converter control parameters.

[0027] Specifically, the global mapping model is read and the set of all RBF nodes {Node_i} and their corresponding weights {w_i} are extracted. The contribution index CI of each RBF node is calculated as: CI_i = |w_i| • ∫_Ωφ_i(x)dx, where |w_i| is the modulus of the weight vector; φ_i(x) is the i-th radial basis function; Ω is the valid region of the input space (jω, S, T), and ∫_Ωφ_i(x)dx represents the integral of the radial basis function over the entire input space. All nodes are sorted in descending order by their contribution CI, and a node importance ranking table T_rank is constructed. An iterative greedy sparsification algorithm is applied: an empty model M_sparse = {} is initialized; nodes are added to M_sparse one by one in the order of T_rank; the performance of the current model on the validation set is evaluated every time K nodes are added; additions are stopped when the performance improvement falls below a threshold Δ. The top N_top nodes in the node importance ranking table T_rank are retained to form a sparse node set, reducing the number of model parameters by 60-85%. The RBF parameters in the sparse node set are quantized using 8-bit fixed-point quantization to further compress storage space and generate a quantized parameter set. The sparse node set and quantized parameter set are combined into a compressed mapping model, reducing the data size by over 80% compared to the original model.

[0028] Read the compressed mapping model and extract the RBF center points {c_i} and coverage radius {r_i=3σ_i}. Determine the boundaries B_min and B_max of the input space (jω, S, T) and construct the root node Root that contains the entire space. Implement the quadtree or octree recursive partitioning algorithm: function BuildQuadTree(Node, Points, Depth): if |Points| ≤ N_thresh or Depth ≥ D_max: Node.Points = Points; return; Divide the Node space evenly into 2 3 = 8 subspaces; for each subspace S_j: Child_j = a new child node; Child_j.Boundary = the boundary of S_j; P_j = an RBF node that falls within or intersects S_j; if P_j is not null: BuildQuadTree(Child_j, P_j, Depth+1); Node.Children.add(Child_j); Execute the quadtree or octree construction algorithm BuildQuadTree(Root, {c_i}, 0) to generate a complete quadtree or octree structure. For each quadtree or octree leaf node, precompute the contribution of the contained RBF nodes to the region and sort them by contribution to generate a node influence table. Construct a spatial hash table to quickly map the input coordinates (jω, S, T) to quadtree or octree leaf nodes, forming a spatial hash index. The quadtree or octree structure, node influence table, and spatial hash index are combined into an indexed mapping model to support millisecond-level query response.

[0029] Read the current operating parameters and indexed mapping model in real time. Use spatial hash indexing to locate the quadtree or octree leaf node L containing the query point with O(1) complexity. Extract the top-K most relevant RBF nodes {Node_top} from the node influence table of L nodes. Calculate the response value of these K nodes to the query point: R_i = w_i • exp(-||x-c_i|| 2 / (2σ_i 2), i∈{Node_top}; quickly sum to obtain the impedance estimate: Z_est = [|Z|_est, φ_est] = Σ R_i; convert the complex impedance value into equivalent circuit parameters: R_s = Re(Z_est); X_L = Im(Z_est); combine the current operating parameters with the current impedance parameter ZPV, current = [R_s, X_L], and apply the photovoltaic inverter control algorithm: D = f_D(ZPV, current, ωcurrent); K_p = f_Kp(ZPV, current, ωcurrent); K_i = f_Ki(ZPV, current, ωcurrent); where D is the duty cycle, and K_p and K_i are the PI controller parameters. The optimal converter control parameters C_converter = {D, K_p, K_i, R_s, X_L} are generated.

[0030] According to one aspect of the present application, generating an MPPT reference trajectory and a short-term power prediction trajectory includes: S21. Collect ambient temperature data, light intensity data, and photovoltaic output power, establish a photovoltaic module thermal balance model and identify model parameters, and output a calibration parameter set and a combined thermal and power model; S22, obtaining the equivalent open-circuit voltage and the equivalent short-circuit current, and calculating the fusion temperature estimate value in combination with the calibration parameter set; S23. Combine the fused temperature estimate with the current light intensity, predict the PV curve characteristics and perform error correction to generate a corrected maximum power point, and form an MPPT reference trajectory and a short-term power prediction trajectory.

[0031] Specifically, the ambient temperature data T_ambient, light intensity data S and photovoltaic output power PPV are collected to establish the thermal balance model of the photovoltaic module: T_module = T_ambient + (S / S0)•[NOCT - 20] + k•PPV; where NOCT is the nominal operating temperature, S0 is the standard light condition (1000W / m 2), where k is the thermal coefficient. A sliding time window adaptive recursive least squares algorithm is applied to process historical operating data to identify model parameters k and NOCT. A self-adjusting forgetting factor mechanism is designed to dynamically adjust the forgetting factor λ based on the rate of environmental change, achieving rapid tracking and stabilization of model parameters. A parameter credibility evaluation metric τ is established, combining the parameter change rate and data fitting error to output the calibration parameter set θ = {k, NOCT, λ, τ}. Based on the PV thermal balance model and parameter identification results, a combined heat and power model is constructed: T_module = f_T(T_ambient, S, PPV, θ); PPV = f_P(VPV, IPV, T_module, S); where f_T is the forward temperature prediction function and f_P is the PV power characteristic function. The thermal balance model is coupled with the electrical characteristic model to form a complete combined heat and power model: M_thermo_elec = {f_T, f_P, θ, thermal_params, electrical_params}; where thermal_params is the thermal characteristic parameter set and electrical_params is the electrical characteristic parameter set. Evaluate the model's prediction accuracy under different operating conditions and generate a model performance report. Output the combined heat and power model.

[0032] Through a unique circuit topology design and partial array sampling technology, the equivalent open-circuit voltage (VOC_eq) and equivalent short-circuit current (ISC_eq) of some components are periodically acquired without interrupting normal system operation. Leveraging the temperature-voltage sensitivity relationship and combining it with the parameter credibility assessment metric τ, a temperature estimate, T_est_V, is calculated as: T_STC + (VOC_eq - VOC_STC) / (βVOC) • W_V(τ), where T_STC is the standard test temperature, VOC_STC is the open-circuit voltage under standard conditions, βVOC is the voltage temperature coefficient, and W_V(τ) is a credibility-based weighting function. In parallel, utilizing the temperature-current sensitivity relationship, a second temperature estimate, T_est_I, is calculated as: T_STC + (ISC_eq / S - ISC_STC / S_STC) / (αISC / S_STC) • W_I(τ), where αISC is the current temperature coefficient and W_I(τ) is a credibility-based weighting function. The confidence-weighted fusion algorithm is applied to combine the VOC and ISC dual-path estimation results to calculate the fused temperature estimate T_module_est.

[0033] Based on the estimated temperature and current light intensity S, single / double exponential photovoltaic models are applied to predict the complete PV curve characteristics {(V, P)}. A non-uniform sampling feature point extraction algorithm is designed to densely sample key areas of the PV curve and extract a set of key feature points {(Vi, Pi)}, including the maximum power point MPP (VMPP, PMPP). A historical prediction error database E_history is constructed, and the deviation between the predicted MPP and the actual tracked MPP is compared to establish an error correction mapping function. The error correction mapping function is applied to correct the current prediction and output the corrected maximum power point parameters (VMPP_corrected, PMPP_corrected). Combining the corrected maximum power point parameters with system operating constraints, a smooth transition trajectory is planned to generate the MPPT reference trajectory R_MPPT and the short-term power prediction trajectory as inputs for subsequent control.

[0034] According to one aspect of the present application, calculating a fusion temperature estimate includes: Read the photovoltaic system topology data and sampling schedule, divide the photovoltaic array into predetermined sub-array areas, configure the bidirectional switch and sampling circuit, generate the sampling control sequence and sub-array configuration parameters, and execute the measurement program on the sub-arrays in sequence based on the sampling control sequence and sub-array configuration parameters; During measurement, the original sampling data set is obtained and statistical filtering is performed to form a filtered sampling data set. The equivalent open-circuit voltage and equivalent short-circuit current are calculated taking into account the mismatch effect between sub-arrays. Based on the equivalent open-circuit voltage, equivalent short-circuit current, current light intensity and calibration parameter set, the temperature estimates of the voltage path and current path are calculated respectively, and then weighted fused and filtered to generate a fused temperature estimate.

[0035] Specifically, the PV system topology data and sampling schedule are read to determine a partial sampling strategy. A partial array isolation sampling circuit is constructed: the N×M array is divided into K subarray regions {SA_1, SA_2, ..., SA_K}; each subarray SA_i is configured with a bidirectional switch SW_i and a sampling circuit SC_i. The main controller MCU schedules the sampling process according to the sampling schedule: a rotating sampling protocol is implemented, sampling each of the K subarrays sequentially within each sampling period T: for i = 1 to K: SW_j = connect to the main power line (j ≠ i); SW_i = connect to the sampling circuit; SC_i.measure(); SW_i = connect to the main power line. An energy buffer circuit is constructed to ensure system output stability during sampling and eliminate the impact of sampling on system energy output. The sampling control sequence and subarray configuration parameters are generated for actual sampling operations.

[0036] Read the sampling control sequence and subarray configuration parameters and perform the sampling operation. When the sampling time arrives, perform the measurement procedure on the selected subarray SA_i: quickly disconnect the load of SA_i and measure the stable voltage value; compensate for the transient response and extrapolate to obtain the true open-circuit voltage VOC_i; short-circuit SA_i with a small resistor and measure the current value; compensate for the resistance drop and extrapolate to obtain the true short-circuit current ISC_i: collect the raw data {VOC_i, ISC_i} of K subarrays i=1 K , forming a raw sampling dataset. A statistical filtering algorithm is applied to detect and remove outliers, generating a filtered sampling dataset. Considering the mismatch effect between subarrays, an array equivalent model is applied for compensation: VOC_eq = f_V({VOC_i}, M_mismatch); ISC_eq = f_I({ISC_i}, M_mismatch), where M_mismatch is the array mismatch model. Combined with the current light intensity S, the normalized short-circuit current ISC_eq / S is calculated, and the equivalent open-circuit voltage VOC_eq and equivalent short-circuit current ISC_eq are output for subsequent temperature estimation.

[0037] Read the equivalent open circuit voltage VOC_eq, equivalent short circuit current ISC_eq, current light intensity S and parameter credibility evaluation index τ. Read the standard test condition data from the system parameter library: T_STC (25℃), VOC_STC, ISC_STC (short circuit current under standard conditions), S_STC (standard light conditions, 1000W / m 2), βVOC (voltage temperature coefficient) and αISC (current temperature coefficient). The first temperature estimate is calculated using the temperature-voltage sensitivity relationship: T_est_V = T_STC + (VOC_eq -VOC_STC) / (βVOC). In parallel, the second temperature estimate is calculated using the temperature-current sensitivity relationship: T_est_I = T_STC + (ISC_eq / S - ISC_STC / S_STC) / (αISC / S_STC). Based on the parameter credibility evaluation index τ, the dynamic weight coefficient is calculated: the weight of the voltage path W_V(τ) = exp(-|τ-τ_V_opt| / σ_V) / (exp(-|τ-τ_V_opt| / σ_V) + exp(-|τ-τ_I_opt| / σ_I)); the weight of the current path W_I(τ) = 1 - W_V(τ). τ_V_opt and τ_I_opt are the optimal credibility points of the voltage method and current method, respectively, and σ_V and σ_I are sensitivity parameters. A soft voting fusion algorithm is applied to calculate the weighted fused temperature: T_module_est = W_V(τ)•T_est_V + W_I(τ)•T_est_I. A Kalman filter is implemented to combine historical temperature estimates to generate a smoothed fused temperature estimate, T_module_est: x*_k = F•x*_k-1 + K•(z_k - H•F•x*_k-1), where z_k is the current fused temperature, x*_k is the filtered temperature estimate, F is the state transition matrix, K is the Kalman gain, and H is the observation matrix. A confidence index, CI_T, is calculated for the temperature estimate, which is used to quantify uncertainty in subsequent PV curve prediction. The fused temperature estimate and confidence index are then output for PV characteristic prediction.

[0038] According to one aspect of the present application, forming an MPPT reference trajectory and a short-term power prediction trajectory includes: Calculate the operating parameters of the photovoltaic module based on the fusion temperature estimate and the current light intensity, solve the photovoltaic equation within the voltage range, and generate the PV curve characteristics; Search for the maximum power point in the PV curve characteristics, filter similar operating condition records from the historical prediction error database, calculate the statistical characteristics of the prediction deviation, and generate the corrected maximum power point through correction function correction; The corrected maximum power point is combined with the pre-stored current system operating point and system operation constraints, and a trajectory planner based on model predictive control is used to solve the optimization problem considering the system dynamic constraints, generating the MPPT reference trajectory and short-term power prediction trajectory from the current operating point to the target maximum power point.

[0039] Specifically, the fusion temperature estimate, confidence index, and current light intensity are read. The standard PV module parameters are read from the system parameter library: I_L, STC (photogenerated current), I_0, STC (reverse saturation current), R_s, STC (series resistance), R_sh, STC (shunt resistance), n (ideality factor), and N_s (number of cells in series). Based on T_module_est and S, the model parameters for the current operating conditions are calculated: I_L = (S / S_STC)•[I_L, STC + α•(T_module_est - T_STC)]; I_0 = I_0, STC•(T_module_est / T_STC). 3 •exp[E_g / (k•T_STC) - E_g / (k•T_module_est)]; R_s = R_s, STC•[1 + β_R•(T_module_est - T_STC)]; R_sh = R_sh, STC•[1 + γ_R•(T_module_est - T_STC)]; where E_g is the silicon bandgap energy and k is the Boltzmann constant. N_V voltage points {V_i} are generated using nonuniform sampling within the voltage range [0, 1.2•VOC_est], with increasing density near the MPP. For each voltage point V_i, the implicit photovoltaic equation is solved to calculate the corresponding current I_i: I = I_L - I_0•[exp((V+I•R_s) / (n•V_t)) - 1] - (V+I•R_s) / R_sh); where I is the output current; V is the output voltage; n is the ideality factor; V_t = k•T / q is the thermal voltage; and q is the electron charge. The solution is iteratively solved using the Newton-Raphson method. The power point sequence P_i = V_i•I_i is calculated to form the predicted PV curve characteristic {(V_i, P_i)}. Taking into account the uncertainty of the temperature estimate, a confidence interval for the PV curve {(V_i, [P_i, min, P_i, max])} is generated based on the confidence index CI_T to quantify the uncertainty range of the prediction.

[0040] Read the predicted PV curve characteristics and PV curve confidence interval. Initial search for the maximum power point: (VMPP, PMPP) = argmax_i{P_i}; Read the historical prediction error database E_history from the system database, which contains the comparison records of the predicted MPP and actual MPP for the past N cycles: E_history = {(VMPP_pred, j, PMPP_pred, j, VMPP_actual, j, PMPP_actual, j)} j=1 N; Based on the current working condition (S, T_module_est), filter the K most similar historical records from E_history to form a similar working condition set S_similar. Calculate the statistical characteristics of the prediction deviation of these K records: μ_V = mean{(VMPP_actual, j - VMPP_pred, j) / (VMPP_pred, j), j∈S_similar}; σ_V = std{(VMPP_actual, j - VMPP_pred, j) / (VMPP_pred, j), j∈S_similar}; μ_P = mean{(PMPP_actual, j - PMPP_pred, j) / (PMPP_pred, j), j∈S_similar}; σ_P = std{(PMPP_actual, j - PMPP_pred, j) / (PMPP_pred, j), j∈S_similar}; where μ_V and μ_P are the mean relative deviations of voltage and power, respectively; σ_V and σ_P are the standard deviations of the relative deviations of voltage and power, respectively; VMPP_pred,j and PMPP_pred,j are the predicted MPP values ​​of the j-th record; VMPP_actual,j and PMPP_actual,j are the actual MPP values ​​of the j-th record; an adaptive Bayesian correction function F_correct is constructed to adjust the current prediction according to the historical deviation pattern: VMPP_corrected =VMPP•(1 + μ_V); PMPP_corrected = PMPP•(1 + μ_P); where VMPP_corrected and PMPP_corrected are the corrected MPP voltage and power, respectively; μ_V and μ_P are the mean relative deviations of voltage and power, respectively. Calculate the confidence interval of the corrected MPP: VMPP_CI = [VMPP_corrected•(1-1.96•σ_V), VMPP_corrected•(1+1.96•σ_V)]; PMPP_CI = [PMPP_corrected•(1-1.96•σ_P), PMPP_corrected•(1+1.96•σ_P)]; output the corrected maximum power point parameters (VMPP_corrected, PMPP_corrected) and the corresponding confidence intervals (VMPP_CI, PMPP_CI) as the targets of MPPT control.

[0041] Read the corrected maximum power point parameters, confidence interval, and current system operating point (V_current, P_current). Read the system operating constraints C_system from the system parameter library, including the maximum voltage change rate (dV / dt)_max, the maximum current change rate (dI / dt)_max, and the power slope limit (dP / dt)_max. Build a trajectory planner based on model predictive control (MPC), considering the system dynamic constraints: min_{u_0, ..., u_N-1} Σ[w_V•(V_k-VMPP_corrected)] 2 + w_P•(P_k-PMPP_corrected) 2 ] + w_u•Σ||u_k|| 2 ;st |V_k+1-V_k| ≤ (dV / dt)_max•Δt; |P_k+1-P_k| ≤ (dP / dt)_max•Δt; V_min ≤ V_k ≤ V_max; where u_k is the control increment, w_V, w_P, and w_u are weight coefficients, N is the prediction time domain length, and Δt is the time step. Solve the MPC problem and obtain the optimal trajectory from the current operating point to the target MPP, forming the MPPT reference trajectory R_MPPT = {(V_k, P_k)} k=0 N Combining the PV curve characteristics and system dynamic characteristics, the power change trend in the future time window is predicted to generate a short-term power prediction trajectory P_predicted = {P_t} t=t_0 t_0+T_pred . Calculate the uncertainty envelope of the power trajectory to form the power prediction confidence interval P_CI = {[P_min, t, P_max, t]} t=t_0 t_0+T_pred Output MPPT reference trajectory R_MPPT and short-term power prediction trajectory P_predicted.

[0042] According to one aspect of the present application, generating a final execution instruction set includes: S31. At the electrical control layer, combining the MPPT reference trajectory and converter control parameters, monitoring the power change rate, and generating electrical fast response control instructions and event characteristics; S32. At the pressure control layer, based on event characteristics and short-term power prediction trajectory, calculate the optimal pressure set point and generate pressure regulation control instructions and regulation status; S33. In the temperature control layer, the long-term optimal temperature set point is calculated by combining the adjustment state with pre-stored historical operation data and weather forecast data to obtain a temperature optimization control instruction; S34. Electrical fast response, pressure regulation and temperature optimization control instructions are stored in a multi-layer data buffer, and a priority allocation rule is adopted to generate a final execution instruction set.

[0043] Specifically, high-frequency sampling data (10kHz) is collected for the photovoltaic output voltage (VPV), current (IPV), and electrolyzer voltage (VEL) and current (IEL). The MPPT reference trajectory is read and the real-time tracking deviation (e_track) is calculated. A wavelet transform denoising algorithm is applied to the raw sampled data to obtain filtered electrical signals (V_filtered and I_filtered). The real-time power (PPV) and power (PEL) are calculated based on converter control parameters and tracking deviation, and the power change rate (dP / dt) is monitored. When a power transient event (E_transient) is detected (dP / dt exceeds a threshold), the fast-response buffer mechanism is triggered, generating an electrical fast-response control command (C_fast). This command (C_fast) is transmitted to the DC-DC converter control unit, while the event signature (F_transient) is uploaded to the medium-speed data processing unit. Medium-frequency sampling data (10Hz) is collected for the electrolyzer pressure (P_EL), hydrogen production (H2_flow), and system power (P_system). The event signature and short-term power prediction trajectory (P_predicted) are received to assess the impact of energy fluctuations. Applying a variable pressure regulation prediction model, the optimal pressure set point, P_EL_set, is calculated based on the current power and predicted power trajectory. A pressure regulation trajectory planner is designed to generate a smooth transition path from the current pressure, P_EL, to the target pressure, P_EL_set, and output a pressure regulation control command, C_pressure. This command is sent to the electrolysis system's pressure control unit, which then uploads the regulation status to a slow data processing unit.

[0044] Low-frequency sampling data (0.1 Hz) is collected for the electrolyzer temperature T_EL, system efficiency η_system, and environmental condition data ENV. The medium-speed layer regulation status and historical operating data are received to assess long-term operating trends. Weather forecast data is obtained through an external meteorological service interface. Combined with historical light patterns, a time series prediction model is applied to generate energy forecast data. A variable temperature optimization decision model is applied, combining weather forecast data and energy forecast data to calculate the long-term optimal temperature set point. A multi-step forward temperature control strategy is designed, accounting for the thermal inertia of the electrolysis system, to generate temperature optimization control instructions. These temperature optimization control instructions are sent to the electrolysis system temperature control unit, and the temperature control status data S_temp is stored for cross-timescale coordination. Three data buffers, B_fast, B_medium, and B_slow, are established to store control instructions and status data for fast response, pressure regulation, and temperature optimization, respectively. A priority assignment rule is designed to dynamically prioritize objectives such as safety constraints and efficiency optimization. A timescale separation and coordination algorithm is applied to detect and resolve conflicts between control instructions from different timescales and generate a coordinated control instruction set. Develop a control translator to convert the coordinated control instruction set into specific control parameters of each execution unit, and output the final execution instruction set for system execution and reinforcement learning training.

[0045] According to one aspect of the present application, in the pressure control layer, the process of generating pressure regulation control instructions and regulation status is specifically as follows: constructing a variable pressure regulation prediction model. Reading the intermediate frequency sampling data of the electrolyzer pressure P_EL, hydrogen production H2_flow, system power P_system, as well as the event feature F_transient and the short-term power prediction trajectory P_predicted. Reading the electrolyzer pressure-efficiency characteristic curve set P_efficiency = {(P_j, η_j(P))} from the system parameter library j=1 N_P, where P_j represents different pressure levels and η_j(P) represents the power-efficiency function at that pressure. Construct a pressure optimization algorithm based on dynamic programming: discretize the prediction time domain [t_0, t_0+T_pred] into M time steps; for each time step k and pressure level j, calculate the state transition cost: Cost(k, j) = -η_j(P_predicted[k]) + α•TransitionCost(P_current, P_j); construct the dynamic programming state transition equation: V(k, j) = min_i{V(k-1, i) + Cost(k, j)}; trace back to obtain the optimal pressure sequence {P_opt}: Considering the uncertainty of power prediction, introduce a robust optimization objective and modify the cost function to: Cost(k, j) = -E[η_j(P_predicted[k]) | P_CI[k]] + α•TransitionCost(P_current, P_j); where Cost(k, j) is the cost of selecting pressure level j at time step k, E[η_j(P_predicted[k]) | P_CI[k]] is the expected efficiency considering power forecast uncertainty; P_CI[k] is the confidence interval of the power forecast at time step k; α is the transition cost weight; and TransitionCost(P_current, P_j) is the cost of changing from the current pressure P_current to the pressure level P_j, which is proportional to the magnitude of the pressure change. Solve the robust optimization problem to obtain the optimal pressure setpoint and future pressure adjustment plan.

[0046] Establish a pressure-efficiency mapping relationship. Read the pressure-power-efficiency data set from the historical operation database, which contains electrolysis efficiency records under different pressure and power conditions. Apply multiple regression analysis to establish a three-dimensional mapping relationship between pressure, power, and efficiency: η(P, Power) = Σ a_ij•P i •Power j, where a_ij is the regression coefficient, determined by solving a least-squares problem. Gaussian process regression (GPR) is applied to the model residuals to establish a nonparametric correction model: Δη(P, Power) ~ GP(μ(P, Power), k((P, Power), (P', Power'))), where k(•,•) is the kernel function to capture complex nonlinear relationships. Combining the parametric model and nonparametric correction terms, a complete pressure-efficiency mapping model M_P_η is constructed: η_total(P, Power) = η(P, Power) + Δη(P, Power). The model prediction accuracy is evaluated for different pressure levels and power ranges, generating a model uncertainty map U_map. Based on the current system power P_system and the short-term power prediction trajectory P_predicted, the pressure-efficiency mapping model M_P_η is queried to obtain the expected efficiency curves at each pressure level, and the efficiency prediction set η_pred = {η_j(t)}_{j, t} is obtained. Based on the pressure-power-efficiency three-dimensional mapping model M_P_η and historical operating data, a complete electrolysis pressure characteristic model is constructed: M_EL_P = {M_P_η, P_efficiency, η_pred, dynamic_params}; where dynamic_params contains the pressure dynamic response characteristic parameters.

[0047] The current pressure P_EL, the optimal pressure setpoint P_EL_set, and system constraints (maximum pressure rate of change (dP / dt)_max and maximum pressure overshoot ΔP_max) are read. A smooth transition trajectory generator based on an S-curve is constructed: the pressure adjustment ΔP is calculated as P_EL_set - P_EL; the time allocation for the acceleration, constant velocity, and deceleration segments is determined; and a continuous and smooth pressure trajectory is generated using cubic spline interpolation. Considering the dynamic response characteristics of the electrolysis system, a feedforward-feedback composite control strategy is introduced: P_command(t) = P_traj(t) + K_p•e(t) + K_i•∫e(t)dt; where e(t) = P_traj(t) - P_actual(t) is the tracking error. The generated trajectory is safety-verified to ensure that system constraints are not violated at any time, generating a safe pressure trajectory P_traj. The safe pressure trajectory P_traj is discretized into a control sequence, which, combined with the control parameters, generates a pressure regulation control command. Calculate the current pressure regulation status, including regulation target, progress and expected completion time, to form a regulation status.

[0048] According to one aspect of the present application, in the temperature control layer, the process of obtaining the temperature optimization control instruction is specifically as follows: Read environmental condition data (ENV) (including temperature, humidity, and cloud cover), historical light data (S_history), and weather forecast data (W_forecast) (for the next 24 hours) provided by an external meteorological service interface. Build a multi-input, multi-output time series prediction model using a long short-term memory (LSTM) network architecture: x_t = [ENV_t, S_t, W_forecast_t, TOD_t, DOW_t]; h_t = LSTM(x_t, h_t-1); S_pred_t = f_S(h_t). Where TOD_t is the time of day and DOW_t is the day of the week, encoding the temporal pattern; x_t is the input feature vector; ENV_t is the environmental condition data; S_t is the current light intensity; W_forecast_t is the weather forecast data; h_t is the LSTM hidden state; S_pred_t is the predicted light intensity; and f_S is the LSTM output function. Use historical data to train the LSTM network, and the optimization goal is to minimize the prediction mean square error: min_θ Σ||S_actual_t - S_pred_t|| 2 θ is a network parameter. Implement an online learning mechanism, continuously updating the model using the latest data from the sliding window to adapt to environmental changes. Convert the predicted light intensity to the expected system power input: P_forecast_t = f_P(S_pred_t, T_ambient_pred_t); where f_P is the PV system power model. Generate energy forecast data for the next 24 hours: E_forecast = {P_forecast_t} t=t_0 t_0+24h , providing long-term energy trends. Calculate the uncertainty estimate of the forecast to form the energy forecast confidence interval E_CI = {[P_min_t, P_max_t]} t=t_0 t_0+24h .

[0049] Read the electrolyzer temperature T_EL, system efficiency η_system, regulation status S_pressure, energy forecast data E_forecast and energy forecast confidence interval E_CI. Obtain the temperature-efficiency characteristic data set D_T_η = {(T_i, η_i(P, Power))} from the system parameter library i=1 N_T , representing the electrolysis efficiency at different temperatures, pressures, and powers. A temperature-power-efficiency three-dimensional mapping model M_T_η is constructed, combining the parameter model and Gaussian process regression: η(T, P, Power) = η_base(T, P, Power) + Δη(T, P, Power); a temperature decision algorithm based on rolling horizon optimization is constructed: max_{T(t)} ∫t0 t0+Th η(T(t), P(t), E_forecast(t))•E_forecast(t)dt; st T_min ≤ T(t) ≤ T_max; |dT / dt| ≤ (dT / dt)_max; T(t0) = T_current; where Th is the optimization time domain (typically 4-8 hours), η is the system efficiency function, E_forecast is the predicted power, and (dT / dt)_max is the maximum temperature change rate constraint. Considering the thermal inertia of the electrolyzer, a temperature dynamic model is established: dT / dt = (Q_in - Q_loss) / C_thermal; where Q_in is the heating power, Q_loss = k•(T-T_ambient) is the heat loss, k is the heat dissipation coefficient, and C_thermal is the system heat capacity. Considering the uncertainty of energy forecasting, a robust optimization framework is introduced, using the conditional value at risk (CVaR) instead of the expected efficiency: max_{T(t)}CVaR_α[∫η(T(t), P(t), E(t)) • E(t)dt]; where α is the confidence level (typically 0.9). Solving the robust optimization problem yields the long-term optimal temperature setpoint and temperature adjustment plan. Based on the three-dimensional temperature-power-efficiency mapping model M_T_η and the system's thermodynamic characteristics, a complete electrolysis temperature characteristic model M_EL_T is constructed: M_EL_T = {M_T_η, thermal_dynamics, heat_transfer_params}; thermal_dynamics describes the dynamic characteristics of the electrolytic cell temperature, and heat_transfer_params is the set of heat transfer parameters. The resulting electrolysis temperature characteristic model M_EL_T is then output for subsequent multi-domain coupled model construction.

[0050] Read the electrolytic cell temperature T_EL, the long-term optimal temperature setpoint and temperature adjustment plan, and the electrolysis system's thermal characteristics. Construct a feedforward controller based on the thermodynamic model, taking into account the system's thermal inertia: dT / dt = (P_heat - k•(T-T_ambient)) / C_thermal; where P_heat is the heating power, k is the heat dissipation coefficient, and C_thermal is the heat capacity. Based on the model and temperature trajectory, calculate the required heating power curve: P_heat(t) = C_thermal•(dT / dt)_desired + k•(T(t)-T_ambient), where (dT / dt)_desired is the desired temperature change rate (°C / s), T(t) is the current temperature, and T_ambient is the ambient temperature. Construct an adaptive model predictive controller (AMPC) to optimize the control sequence: min_{u_0,...,u_N-1} Σ[w_T•(T_k-T_traj_k) 2 + w_u•||u_k|| 2 ]; st T_k+1 = f(T_k, u_k, T_ambient_k); u_min ≤ u_k ≤ u_max; |T_k+1-T_k| ≤ ΔT_max; where u_k is the heating power control variable and T_traj_k is the target temperature trajectory. An adaptive mechanism is incorporated to online identify model parameters and compensate for external disturbances and model errors: θ_k+1 = θ_k + γ•Φ_k•(T_measured_k - T_predicted_k), where θ is the model parameter, Φ is the sensitivity matrix, and γ is the learning rate. Combining feedforward and feedback control strategies, robust temperature optimization control instructions are generated, consisting of a sequence of heating power setpoints and control parameters. The current temperature regulation state, including regulation target, progress, and expected efficiency improvement, is calculated to generate temperature control state data. The temperature optimization control instructions and temperature control state data are output for temperature control and cross-timescale coordination of the electrolysis system.

[0051] According to one aspect of the present application, the step of generating a final execution instruction set includes: Construct a three-layer data buffer structure, including a millisecond-level ring buffer, a second-level linked list structure, and a minute-level array structure, to implement a cross-buffer data synchronization mechanism, generate synchronized data streams and system status time series, and form a multi-scale buffer dataset; Construct the system state space and objective function set, evaluate the importance of each objective and the risk of conflict, implement the priority adaptive adjustment strategy, and dynamically generate priority allocation rules and objective importance evaluation; Identify potential conflicts between control instructions at different time scales, construct a conflict resolution space, implement a control instruction coordinator, generate a coordinated control instruction set and convert it into the final execution instruction set through a control translator.

[0052] Specifically, read the electrolysis pressure characteristic model M_EL_P and the electrolysis temperature characteristic model M_EL_T. Analyze the interaction between pressure and temperature, establish the coupling matrix C_PT, and describe the influence of temperature on pressure effects and the influence of pressure on temperature dynamics. Integrate the pressure and temperature models to construct a unified electrolysis temperature-pressure characteristic model M_EL_PT: M_EL_PT = {M_EL_P, M_EL_T, C_PT, coupling_functions}, where coupling_functions is a set of functions describing the temperature-pressure coupling relationship. Verify the prediction accuracy of the integrated model under different working conditions and generate a model verification report. Output the electrolysis temperature-pressure characteristic model M_EL_PT for use in multi-domain coupling environment model construction.

[0053] Read the electrical fast-response control command C_fast, the pressure regulation control command C_pressure, and the temperature optimization control command C_temperature. A three-layer data buffer structure is constructed: B_fast: a circular buffer with a capacity of 1000 entries and an update frequency of 10kHz; B_medium: a linked list structure with a capacity of 100 entries and an update frequency of 10Hz; and B_slow: an array structure with a capacity of 24 entries and an update frequency of 0.1Hz. A cross-buffer data synchronization mechanism is implemented to ensure that upper-layer decisions are passed to lower-layer execution: control commands in B_slow are downsampled and updated for related items in B_medium; control commands in B_medium are downsampled and updated for related items in B_fast; status data in B_fast is upsampled and updated for status monitoring in B_medium; status statistics in B_medium are upsampled and updated for trend analysis in B_slow. A data time tagging and synchronization mechanism is implemented to solve the time alignment problem of data at different time scales and generate a synchronized data stream. System status history is maintained, supporting status backtracking and trend analysis, and constructing a system status time series. The output is the integrated multi-scale buffer dataset, which is provided to the conflict detection and resolution module.

[0054] Read the multi-scale buffer dataset and system configuration parameters, including basic priority settings and constraints. Define the system state space S and the set of objective functions {f_i}, including multiple objectives such as maximizing efficiency, ensuring safety, and protecting equipment lifespan. Construct a priority learning algorithm based on hierarchical reinforcement learning: Q(s, a) = E[r + γ•max_a' Q(s', a')], where r is the immediate reward, reflecting the weighted combination of the objective functions, γ is the discount factor, and Q(s', a') is the Q value after selecting action a' in the next state s'. Based on the current system state, dynamically evaluate the importance and conflict risk of each objective: I_i(s) = w_i•f_i(s) + β•Conflict_i(s), where w_i is the objective weight, Conflict_i indicates the degree of conflict with other objectives, f_i(s) is the normalized value of the objective function in state s, and β is the conflict risk coefficient. Implement an adaptive priority adjustment strategy, adjusting priorities based on state sensitivity: P_i(s) = softmax(I_i(s) / τ) = exp(I_i(s) / τ) / Σexp(I_j(s) / τ), where τ is a temperature parameter that controls the degree of "softness" in the allocation; Σexp(I_j(s) / τ) is the exponential sum of all objectives, used for normalization. Predefined priority templates are used for different operating conditions: Normal operating condition: efficiency > stability > response speed; Emergency condition: safety > stability > efficiency; Fault recovery: reliability > efficiency > response speed. Based on the current system state, the most suitable operating condition is identified and a priority allocation rule is generated. The priority allocation rule and objective importance assessment are output for subsequent conflict resolution.

[0055] The system reads a multi-scale buffer dataset, assigns priorities, and evaluates target importance. A conflict detection algorithm is constructed to identify potential conflicts between control instructions at three timescales: resource contention conflicts (multiple control instructions compete for the same execution resource); target conflict conflicts (different control objectives constrain each other); and timing dependency conflicts (execution order influences the final outcome). For each detected conflict C_j, a conflict resolution space A_j is constructed, and the utility of each solution is evaluated: U(a) = Σ P_i•f_i(a), where P_i is the priority of target i and f_i(a) is the contribution of solution a to target i. A constraint satisfaction-based conflict resolution algorithm is implemented: a control variable constraint network is constructed; an arc consistency algorithm is applied to streamline the solution space; and a heuristic search is used to find the optimal solution that satisfies the constraints. A control instruction coordinator is constructed to convert the solutions into specific execution instructions: timing coordination (adjusting the execution order of instructions); parameter coordination (modifying instruction parameters to eliminate conflicts); and target coordination (finding a balance between conflicting objectives). Finally, a coordinated control instruction set is generated to ensure harmonious control collaboration across different timescales. Build a control translator to convert the coordinated abstract instructions into control signals for specific execution units: converter control signals: PWM waveform, switching timing; pressure control signals: valve opening, pump speed setting; temperature control signals: heating power, cooling flow; output the final execution instruction set C_final for system execution and reinforcement learning training.

[0056] According to one aspect of the present application, generating a system comprehensive control instruction includes: S41. Based on the impedance frequency characteristic mapping model and the thermoelectric joint model, a photovoltaic-electrolysis multi-domain coupling model is constructed and encapsulated as a reinforcement learning environment interface; S42, based on the reinforcement learning environment interface, train the knowledge-enhanced policy network and output the optimized policy network and global decision; S43. Based on the global decision, define the local state space and action space, and train the adaptive execution network; S44. Read the current system state, generate a control parameter set by inference through the optimization strategy network and the adaptive execution network, fuse it with the final execution instruction set, and output the system comprehensive control instruction.

[0057] Specifically, the following models and data were read: a combined thermoelectric model, an electrolysis temperature-pressure characteristic model, and an anti-frequency characteristic mapping model. Model interfaces were standardized to ensure consistency and compatibility of input and output variables across models. Physical units and coordinate systems were standardized; a variable mapping table between models was established; and a standard format for model data exchange was designed. Preliminary integration testing was performed to verify the correctness of data transfer between models and generate an integration test report. The preprocessed model set and interface specifications were output for subsequent detailed multi-domain coupling model construction.

[0058] The system integrates a combined thermoelectric model, an electrolysis temperature and pressure characteristic model, and system topology data. A photovoltaic-electrolysis multi-domain coupling model (MCM) is constructed. An adaptive physical parameter identification algorithm is applied to dynamically update key parameters in the MCM based on real-time operating data, resulting in a calibrated multi-domain coupling model. A model verification and credibility assessment mechanism is designed to regularly compare model predictions with measured values ​​and calculate model credibility metrics. Based on the model credibility metrics, the different submodules of the MCM are adaptively adjusted to output an optimal multi-domain coupling model. This optimal multi-domain coupling model is encapsulated as a reinforcement learning environment interface, RL_ENV, providing a state-action-reward interaction mechanism for subsequent reinforcement learning training. The state space, S_global, is defined (including global information such as illumination, temperature, system efficiency, and constraint status) and the action space, A_global, is defined (including high-level decisions such as power allocation ratio and operating mode selection). A hierarchical reward function, R_global, is constructed, taking into account system efficiency, equipment lifespan, and constraint satisfaction. The global policy is trained using the PPO algorithm framework through the reinforcement learning environment interface. The policy is updated using a truncated objective function. A neural-symbolic hybrid architecture is introduced to embed physical knowledge as a priori constraints into the policy network, forming a knowledge-enhanced policy network. A batch-efficient experience replay mechanism is designed to select the most informative samples for training, accelerate network convergence, and output an optimized policy network.

[0059] Based on the global decision output by the global policy network, the local state space and local action space (including specific control parameters) are defined. A multi-objective local reward function is designed to precisely evaluate control performance. Through the reinforcement learning environment interface, the SAC algorithm framework is used to train the local execution network (LEN) based on historical control data. Double Q learning is used to reduce the bias in the value function estimation. The current system state is read and the global decision is calculated through inference using the optimized policy network. Based on the global decision and the current local state, the specific control parameter set is generated through inference using the adaptive execution network. A feedforward-feedback hybrid control strategy is applied, combining model predictive control with feedback control to enhance robustness and output reinforcement learning control instructions. A real-time inference optimization accelerator is designed, applying model pruning and quantization techniques to ensure control instruction generation latency of less than 5ms, providing real-time control capabilities. The reinforcement learning control instructions are integrated with the final execution instruction set to output the system's integrated control instructions. These are distributed to each execution unit (DC-DC converter, electrolysis pressure controller, and electrolysis temperature controller) to achieve coordinated control of the photovoltaic off-grid hydrogen production system. Execution feedback data is collected and stored in the system operation database for model optimization and control adjustments in the next cycle.

[0060] According to one aspect of the present application, the process of encapsulating the interface of the reinforcement learning environment is: reading the impedance frequency characteristic mapping model, the thermoelectric combined model, the electrolysis temperature and pressure characteristic model, and the system topology data. Build a modular model integration framework and define a unified data exchange interface: module PV_Model: inputs: jω, S, T; outputs: ZPV, VOC, ISC, VMPP, PMPP; module DC_DC_Model: inputs: ZPV, D, L, C; outputs: ZDC, efficiency; module Electrolyzer_Model: inputs: VDC, IDC, P, T; outputs: H2_flow, efficiency, dT / dt; implement the data flow pipeline between models, handle variable type conversion and unit unification: the output data structure of the photovoltaic model PV_out = PV_Model(ω, S, T); the input data structure of the DC-DC model DC_in = DataAdapter(PV_out); the output data structure of the DC-DC model DC_out = DC_DC_Model(DC_in); the input data structure of the electrolyzer model EL_in = DataAdapter(DC_out); the output data structure of the electrolyzer model EL_out = Electrolyzer_Model (EL_in); DataAdapter handles data structure conversions; integrates all sub-models to build a multi-physics model of the system, capturing the coupled relationships between photovoltaics, power electronics, and electrolysis. Build a model solver using a hybrid solution strategy: an explicit Runge-Kutta method for electrical fast dynamics; an implicit Backward Euler method for electrolysis chemical reactions; and an adaptive step-size method for thermodynamic processes. Implement a parallel computing architecture for sub-models to improve the computational efficiency of multi-domain coupled models and build a parallel solution engine.

[0061] Read real-time operating data and the system's multi-physics model. Build a recursive extended Kalman filter (REKF) and identify model parameters online: x*_k = f(x*_k-1, u_k-1) + K_k•(y_k - h(x*_k)); θ*_k = θ*_k-1 + L_k•(y_k - h(x*_k)), where x is the state variable, θ is the parameter vector, K and L are gain matrices, x*_k is the state estimate at time k, f(•) is the state transfer function, u_k-1 is the control input at time k-1, K_k is the state estimation gain matrix, y_k is the observation value at time k, h(•) is the observation function, θ*_k is the parameter estimate at time k, and L_k is the parameter estimation gain matrix. Implement a multi-timescale parameter identification strategy for fast and slow-changing parameters: electrical parameters (L, C, R): high-frequency identification (~1kHz); electrolytic parameters (k_0, α, E_a): medium-frequency identification (~1Hz); thermal parameters (C_p, h): low-frequency identification (~0.01Hz); build a parameter constraint mechanism to ensure that the identification parameters are within a physically reasonable range: θ*_constrained = min(max(θ*, θ_min), θ_max); implement a parameter change detection algorithm to identify system degradation or operating condition changes: Δθ = |θ*_k- θ*_k-N| / |θ*_k-N|; if Δθ>threshold: trigger_update(); update the parameters in the system multi-physics model to obtain a calibrated multi-domain coupled model.

[0062] Read the calibration multi-domain coupling model and the model validation and credibility assessment results. Based on the model credibility index, adaptively adjust different submodules and output the optimal multi-domain coupling model. Build an OpenAI Gym-compatible environment interface and encapsulate the system model: python; class PVH2Env(gym.Env): def __init__(self, model_path): self.model = load_model(model_path); self.observation_space = spaces.Box(...); self.action_space = spaces.Box(...); ..; def step(self, action): next_state = self.model.simulate(self.state, action); reward = self.calculate_reward(self.state, action, next_state); done = self.is_terminal(next_state); return next_state, reward, done, {}; ..; def reset(self): self.state = self.sample_initial_state(); return self.state; next_state is the next state after executing action; self.model.simulate is the system model simulation function; self.state is the current state; action is the current action; reward is the reward obtained by executing the action; self.calculate_reward is the reward calculation function; done is the termination flag; self.is_terminal is the termination state determination function. Implement an efficient simulation engine that supports batch parallel simulation to accelerate reinforcement learning training: def batch_simulate(states, actions, n_steps): # Parallel simulation of n_batch trajectories, each executing n_steps steps; return parallel_results.Build a differentiable reward function calculator to support multi-objective reinforcement learning: def calculate_reward(state, action, next_state): r_efficiency = w1 * system_efficiency(next_state); r_stability = w2 * stability_measure(state, next_state); r_safety = w3 * safety_margin(next_state); return r_efficiency + r_stability + r_safety; where r_efficiency is the efficiency reward component; w1 is the efficiency weight; system_efficiency is the system efficiency calculation function; r_stability is the stability reward component; w2 is the stability weight; stability_measure is the stability measure function; r_safety is the safety reward component; w3 is the safety weight; and safety_margin is the safety margin calculation function. Develop an environment randomization mechanism to enhance policy generalization: def randomize_env_params(): # Randomize model parameters within a physically reasonable range; return randomized_params; Output the reinforcement learning environment interface.

[0063] According to one aspect of the present application, outputting an optimization strategy network and a global decision includes: Construct a three-layer reward function, including an efficiency reward function to evaluate energy conversion efficiency, a stability reward function to penalize severe fluctuations, and a safety reward function to ensure the safe operation of the system. This implements a dynamic weight adjustment mechanism to form a hierarchical reward function. The policy network is trained based on the PPO algorithm framework, maintaining the posterior distribution of the environment dynamics model. In each training cycle, multiple model instances are sampled from the posterior distribution and the corresponding sub-strategies are optimized. Multiple sub-strategies are mixed to form an exploration strategy, achieving a progressive exploration-exploitation balance mechanism. Construct a neural-symbolic hybrid architecture, combine connectionist and symbolic methods, implement a neural network policy constraint mechanism, build a differentiable logic rule encoder, convert symbolic rules into trainable neural network representations, inject expert knowledge through the knowledge distillation process, and generate an optimized policy network and global decision.

[0064] Specifically, the system optimization objective parameters and constraints are read, and the state space and action space are defined. A hierarchical reward function structure is constructed: R_total(s, a, s') = w_eff • R_efficiency(s, a, s') + w_stab • R_stability(s, a, s') + w_safe • R_safety(s, a, s'); where R_total is the total reward; s is the current state; a is the current action; s' is the next state; w_eff = 0.6 is the efficiency reward weight; R_efficiency is the efficiency reward function; w_stab = 0.3 is the stability reward weight; R_stability is the stability reward function; w_safe = 0.5 is the safety reward weight; and R_safety is the safety reward function. Implement an efficiency reward function to evaluate energy conversion efficiency: R_efficiency(s, a, s') = η_PV_to_H2(s') - η_baseline, where η_PV_to_H2 is the total conversion efficiency from photovoltaic to hydrogen energy, and η_baseline is the baseline efficiency. Construct a stability reward function to penalize severe fluctuations: R_stability(s, a, s') = -α•||s'-s|| 2 - β•||a-a_prev|| 2 , where α and β are weight coefficients, ||s'-s|| 2 is the square norm of state change; ||a-a_prev|| 2 is the square norm of the difference between the current action and the previous action. Construct a safety reward function to ensure that the system operates within the safety boundary: R_safety(s, a, s') = -∞ if the safety constraint is violated else -γ•distance_to_boundary(s'), where γ is the weight coefficient and distance_to_boundary calculates the distance from the state to the safety boundary. Implement a dynamic weight adjustment mechanism to adaptively adjust the weights of each part according to the system state: w_i(s) = base_w_i +adjust_factor(s, i), where w_i(s) is the weight of target i in state s; base_w_i is the base weight; adjust_factor is the state-based adjustment function. Construct a temporal credit allocation mechanism for sparse rewards to solve long-term planning problems: R_discounted(t) = Σ γ k • R(t+k), where γ is the discount factor. Outputs the complete hierarchical reward function for policy network training.

[0065] Read the reinforcement learning environment interface and hierarchical reward function. Implement the policy network training process based on the PPO algorithm framework: use the current policy π_θ to sample N trajectories; calculate the advantage function A π (s, a); construct PPO objective function: J(θ) = E[min(r(θ)•A π (s, a), clip(r(θ), 1-ε, 1+ε)•A π (s, a))], where r(θ) = π_θ(a|s) / π_θ_old(a|s). Stochastic gradient ascent is used to maximize J(θ), constructing an exploration strategy based on posterior sampling to maintain the posterior distribution of the environment dynamics model: p(M|D) ∝ p(D|M)•p(M), where p(M|D) is the posterior distribution of model M given the data D; p(D|M) is the likelihood function; and p(M) is the model prior distribution. In each training cycle, K model instances are sampled from the posterior distribution: {M_i}_{i=1} K ~ p(M|D). For each sampled model Mi, optimize sub-policy πi: πi = argmax_π E_{M_i}[R(τ)|τ~π]; E_{M_i}[R(τ)|τπ] is the expected reward of executing policy π under model Mi. K sub-policies are mixed to form a final exploration policy: π_explore(a|s) = Σ w_i•π_i(a|s), where w_i is the mixing weight, reflecting the posterior probability of the corresponding model. Progressive exploration is implemented by utilizing a balancing mechanism to reduce the exploration ratio as training progresses: the final policy π_final(a|s) = (1-α_t)•π_θ(a|s) + α_t•π_explore(a|s); α_t is the exploration weight, which decays over time. The trained policy network π_θ is output as the basis for global decision making.

[0066] The policy network and the system's physical knowledge base, including constraints and expert experience, are accessed. A neural-symbolic hybrid architecture is constructed, combining connectionist and symbolic approaches: the neural network component handles perception and decision-making; the symbolic reasoning component applies physical constraints and domain knowledge; the interface layer connects the two components and coordinates information flow. The neural network policy constraint mechanism is implemented to ensure that the output complies with physical laws: a_raw = NN_policy(s); a_constrained = ConstraintSatisfaction(a_raw, K_physics), where a_raw is the raw action output by the neural network; NN_policy is the policy network forward function; a_constrained is the action after applying constraints; ConstraintSatisfaction is the constraint satisfaction function; and K_physics is the physical knowledge base. A differentiable logic rule encoder is constructed to convert symbolic rules into a trainable neural network representation: R(x) = sigmoid(β•(f(x) - threshold)), where R(x) is the rule activation and f(x) is the rule condition function. Implement the knowledge distillation process to inject expert knowledge into the neural network: distillation loss L_distill = KL(π_expert(•|s) || π_θ(•|s)), where KL is the KL divergence and π_expert is the rule-based expert policy. A dynamic node expansion mechanism is built to adaptively adjust the network structure based on task requirements: if error_region(s) > threshold: add_new_symbolic_nodes(region(s)). All components are integrated to construct the knowledge-enhanced policy network GPN+, combining data-driven learning with prior physical knowledge. Monte Carlo cross-validation is used to evaluate the performance of GPN+ under different working conditions, generating a performance analysis report. Model pruning and quantization techniques are applied to improve inference efficiency while maintaining performance, resulting in an optimized policy network.

[0067] According to one aspect of the present application, training an adaptive execution network includes: Based on the global decision, define the local state space containing device-level state variables and the local action space containing specific control parameters, establish the mapping relationship between the global decision and the local parameters, build a parameter normalization and scaling mechanism, realize parameter grouping and decoupling, and form a local parameter space definition; Extracting physical gradient information from the multi-domain coupling model, combined with the SAC algorithm framework, constructing a gradient-guided exploration strategy, combining random exploration with physical gradient information, implementing a gradient priority sampling mechanism, and constructing a hybrid exploration strategy of Thompson sampling and gradient guidance to generate an efficient execution network; Build a sliding time window data buffer, implement an adaptive learning rate adjustment mechanism, adopt a multi-model integration architecture to simultaneously maintain multiple execution network versions, build a model selector, knowledge transfer mechanism and incremental learning framework, realize continuous learning and adaptive adjustment, and form an adaptive execution network.

[0068] Specifically, the global decision output from the optimization strategy network and the system control parameter configuration are read. A local state space, S_local, is defined, containing device-level state variables: s_local = [V_PV, I_PV, V_DC, I_DC, V_EL, I_EL, P_EL, T_EL, ...]; and a local action space, A_local, is defined, containing specific control parameters: a_local = [D_converter, f_switching, K_p_V, K_i_V, K_p_I, K_i_I, ...], where D_converter is the converter duty cycle, f_switching is the switching frequency, and K_* is the controller parameter. A mapping relationship, Φ, between global decisions and local parameters is established: D_global → {constraints on A_local}, where {constraints on A_local} is the set of constraints in the local parameter space. Global decisions define the feasible domain of local parameters. Implement parameter normalization and scaling to map physical parameters to the interval [-1, 1]: the normalized parameter value a_norm = 2•(a - a_min) / (a_max - a_min) – 1, where a is the original parameter value; a_min and a_max are the minimum and maximum values ​​of the parameter, respectively. Build a parameter correlation analyzer to identify coupling relationships between parameters: C_ij = correlation(a_i, a_j). Construct a parameter correlation matrix C. Based on the correlation matrix, group and decouple parameters to simplify the optimization problem: {G_k} k=1 M = parameter_grouping(C, threshold), where G_k is the group of related parameters. This outputs the definition of the local parameter space, including the state space, action space, and constraints, which is used to perform network training.

[0069] Read the local parameter space definition and multi-domain coupling model. Build an execution network training framework based on the SAC algorithm: maintain the policy network π_φ and two Q networks Q_θ1, Q_θ2; sample batch data from the replay buffer; update the Q network: L_Q(θ) = E[(Q_θ(s, a) - (r + γ•(min_j Q_θj'(s', a') - α•log π_φ(a'|s')))) 2]; update the policy network: L_π(φ) = E[α•log π_φ(a|s) - min_j Q_θj(s,a)]; extract physical gradient information from the multi-domain coupled model: g_physics(s,a) = ▽_a R(s,a,s'), where s'=f(s,a) is the next state predicted by the model, ▽_a represents the gradient operator for action a, and R is the reward function. Construct a gradient-guided exploration strategy that combines random exploration and physical gradient information: a_explore = a_mean + σ•(ε + β•normalize(g_physics)), where ε~N(0,I) is standard normal noise, β is the physical gradient weight, a_mean is the action mean output by the policy network, and σ is the standard deviation of the exploration noise. Implement a gradient priority sampling mechanism to increase exploration of high-gradient regions: p(s, a) ∝ exp(η•||g_physics(s, a)||), where p(s, a) is the sampling strategy, η is the priority coefficient, and ||g_physics(s, a)|| is the norm of the physical gradient. Construct an adaptive exploration temperature adjustment mechanism: α_t+1 = α_t + λ•(H_target - H_t), where α_t is the temperature parameter at time t; H_t = -E[log π_φ(a|s)] is the policy entropy, and H_target is the target entropy. Implement a hybrid exploration strategy combining Thompson sampling and gradient guidance: p(a|s) = (1-ω)•π_φ(a|s) + ω•Δ(a-a_explore), where ω is the exploration weight and Δ is the Dirac function. Output is the trained efficient execution network LEN+, which has gradient-guided exploration capabilities and improves sample efficiency.

[0070] Read the efficient execution network LEN+ and historical operation data and historical control data. Build a sliding time window data buffer to keep the latest N days of operation data: B_sliding = {(s_t, a_t, r_t, s_t+1)} t=t-N•T t ; where T is the number of data points per day. Implement an adaptive learning rate adjustment mechanism to adjust the update step size according to performance changes: η_t = η_base•exp(-λ•|R_t - R_t-1| / R_t-1), where η_base is the base learning rate; λ=10 is the adjustment coefficient, and R_t is the current performance indicator. Build a multi-model integration architecture and maintain multiple execution network versions at the same time: {LEN_i} i=1 K, each optimized for different operating conditions. Implement a model selector to select the most suitable model based on the current operating conditions: selected model index i* = argmax_isimilarity(current_condition, condition_i); where similarity is the similarity function; current_condition is the current environment condition; and condition_i is the optimization condition corresponding to model i. Build a knowledge transfer mechanism to share and transfer learning experience between models: target model parameters θ_target = θ_target + τ•(θ_source - θ_target), where θ_source is the source model parameter and τ is the migration rate. Implement an anomaly detection and recovery mechanism to trigger retraining when performance degrades significantly: if R_t < γ•mean(R_{tk:t-1}): trigger_retrain(), where R_t is the performance indicator at time t; γ = 0.8 is the threshold coefficient; mean(R_{tk:t-1}) is the mean of the performance indicator over the past k time steps; trigger_retrain() is the function that triggers retraining. An incremental learning framework is constructed to adapt to new situations without forgetting old knowledge: the total loss function L_total = L_current + λ•L_distill, where L_distill = KL(π_old||π_new) is the distillation loss. All adaptive mechanisms are integrated to output an adaptive execution network with continuous learning capabilities, which can continuously optimize control strategies based on environmental changes.

[0071] According to one aspect of the present application, by optimizing the strategy network and the adaptive execution network, outputting system comprehensive control instructions includes: The current system state is input into the optimization strategy network, and the global decision is calculated through forward propagation, including the power allocation ratio and working mode selection; The global decision and the current local state are input into the adaptive execution network, and a specific control parameter set is generated through forward calculation, including the converter duty cycle, controller parameters and equipment operating parameters; Based on a specific set of control parameters, a feedforward-feedback hybrid control strategy is applied, combining model predictive control with feedback control to generate reinforcement learning control instructions. The reinforcement learning control instructions are weighted and fused with the final execution instruction set. The weight of each instruction is dynamically adjusted according to the current system state, and the system comprehensive control instructions are output.

[0072] This embodiment addresses the challenge of multi-domain coupling by employing impedance-frequency mapping technology to establish a four-dimensional mapping relationship between illumination, temperature, frequency, and impedance, effectively capturing the dynamic electrical characteristics of photovoltaic arrays under different operating conditions. Furthermore, a combined thermoelectric self-correcting MPPT technology is designed. Through indirect temperature estimation using electrical parameters, accurate thermodynamic state prediction is achieved without temperature sensors, breaking the information barrier between the photovoltaic electrical and thermodynamic domains. To address the multi-timescale coordination problem, a three-layer collaborative control architecture is constructed: millisecond-level electrical fast response, second-level pressure medium-speed regulation, and minute-level temperature slow optimization. By designing a dynamic priority allocation algorithm and a command conflict detection and coordination mechanism, the system intelligently resolves potential conflicts between control commands at different timescales, ensuring smooth transitions and coordinated operation of various physical processes. A two-layer deep reinforcement learning architecture separates the global policy network from the local execution network, enabling the system to simultaneously address macro-level energy optimization decisions and precise micro-level parameter control. An exploration strategy based on posterior sampling and a gradient-guided exploration mechanism improve the sample efficiency and convergence speed of reinforcement learning in complex dynamic environments.

[0073] In a specific embodiment of this application, efficient coordinated control is achieved for an off-grid hydrogen production system consisting of a 10kW photovoltaic array, a 5kW DC-DC converter, and a 3kW PEM electrolyzer. The system parameters are as follows: 10kW photovoltaic array (40 250W modules), open-circuit voltage VOC = 45.6V, short-circuit current ISC = 8.9A / module; 5kW DC-DC converter, input voltage range 30-60V, output voltage range 30-80V; 3kW PEM electrolyzer, operating voltage range 40-80V, operating temperature range 40-70°C, and operating pressure range 0.1-3.0MPa.

[0074] Step 1: Collect the time domain signals of the photovoltaic array voltage and current and construct an impedance frequency characteristic mapping model.

[0075] 1.1. In this embodiment, the impedance spectrum of the photovoltaic array is measured using the small signal injection method. 2 , step size 100W / m 2 ) and 16 temperatures (0-75°C, step size 5°C), the impedance measurement was performed at 30 frequency points (0.1Hz-10kHz, logarithmic distribution) for each operating point. The specific implementation process is as follows: Under standard test conditions (STC, light intensity 1000W / m 2, temperature 25°C), a small amplitude (5% of the output current) sinusoidal disturbance signal is injected into the PV array, and the frequency is swept from 0.1Hz to 10kHz. The voltage response VPV(t) and current response IPV(t) of the PV array are collected at a sampling rate of 10kHz. A fast Fourier transform (FFT) is performed on the collected time domain signal to obtain the voltage VPV(jω) and current IPV(jω) in the frequency domain. The impedance is calculated: ZPV(jω) = VPV(jω) / IPV(jω), and the impedance modulus |ZPV(jω)| and the phase angle φPV(jω) are extracted. For example, when the light intensity S = 800W / m 2 Under the conditions of temperature T = 30°C and frequency f = 1kHz, the following values ​​are measured: VPV(j2π•1000) = 0.582∠-15.3° V (AC small signal component); IPV(j2π•1000) = 0.126∠-58.7° A (AC small signal component); and ZPV(j2π•1000) = VPV(j2π•1000) / IPV(j2π•1000) = 4.619∠43.4° Ω = 3.361 + j3.171 Ω. This method yields an impedance characteristic dataset {ZPV(jω, S, T)}, totaling 10 × 16 × 30 = 4800 data points. This dataset constitutes the original impedance spectrum database.

[0076] 1.2. Construction of impedance frequency characteristic mapping model.

[0077] 1.2.1. Organize the raw impedance spectrum data into a fifth-order tensor, N_ω=30 (frequency points), N_S=10 (light intensity points), N_T=16 (temperature points), and 2 represents the impedance modulus and phase. Use the Tucker decomposition algorithm to decompose the tensor: X ≈ G × 1U 1 ×2U2×3U 3 ; X is the original data tensor with dimensions of 30×10×16×2; G is the core tensor with dimensions of r1×r2×r3×2; U 1 ∈R N_ω×r1 is the factor matrix in the frequency dimension; U 2 ∈R N_S×r2 is the factor matrix in the illumination dimension; U 3 ∈R N_T×r3 is the factor matrix along the temperature dimension; ×ᵢ is the tensor product operation along the i-th dimension. A cross-validation method was used to determine the optimal rank (r1, r2, r3) = (8, 4, 5), achieving a data compression rate of approximately 93.3% while maintaining a reconstruction error less than 3%. For example, in the original tensor X, the illumination S = 800W / m 2, the impedance value at temperature T = 30°C and frequency f = 1kHz is ZPV = 4.619∠43.4°Ω. After Tucker decomposition, the reconstructed value through the core tensor and factor matrix is ​​ZPV_reconstructed = 4.583∠42.9°Ω, with a relative error of less than 2%.

[0078] 1.2.2 Radial Basis Function Network Structure Design. Based on the decomposed data, a three-layer radial basis function network (RBF network) was constructed: input layer: 3 nodes corresponding to frequency ω, light intensity S, and temperature T; hidden layer: initially configured with 120 RBF nodes; output layer: 2 nodes corresponding to impedance modulus |ZPV| and phase angle φPV. The initial RBF center points ci were determined using the K-means++ algorithm, ensuring that the center point distribution matched the data density. For each center point ci, its RBF width σi is determined by the local density adaptive method: σᵢ = β•mean_{x∈N_k(cᵢ)}||x-cᵢ||; where σᵢ is the width parameter of the i-th RBF node; β is the scaling factor, which is 0.5; N_k(cᵢ) represents the k nearest neighbor sample points of the center point cᵢ, with k being 5; mean_{x∈N_k(cᵢ)}||x-cᵢ|| is the average distance between these k nearest neighbors and the center point cᵢ; ||•|| represents the Euclidean distance. For example, for the center point c1 = [2π•1000, 800, 30] (corresponding to a frequency of 1kHz and an illumination of 800W / m 2 , temperature 30°C), the average distance of its five nearest neighbors to c1 is 76.5 (normalized comprehensive distance), then σ1 = 0.5 × 76.5 = 38.25.

[0079] 1.2.3. Training and validation of the adaptive mapping model. The RBF network was trained using the weighted least squares method combined with an adaptive growing-pruning algorithm. The network size was optimized from the initial 120 RBF nodes to 93 nodes. The performance indicators on the validation set were: RMSE = 0.187Ω (impedance modulus) and 1.53° (phase angle), R 2 = 0.985. The trained global mapping model contains 93 RBF nodes, each characterized by a center point ci, a width σi, and a weight wi. For one of the most important RBF nodes, its parameters are: center point c1 = [2π·1000, 800, 30] (before normalization); width σ1 = 38.25; and weight w1 = [0.873, 0.651] (corresponding to the contributions of the impedance modulus and phase angle, respectively).

[0080] 1.3. Mapping model compression and real-time query optimization.

[0081] 1.3.1. Model Sparsification and Key Node Extraction. To improve the model's real-time performance, the global mapping model is sparsified. The contribution index (CI) of each RBF node is calculated. For example, for node c1 = [2π • 1000, 800, 30], its contribution index (CI1) is |[0.873, 0.651]| • 18732 = 19877.5. All nodes are sorted in descending order of contribution, and the top 35 nodes (approximately 38% of the total number of nodes) are retained to form a sparse node set, N_sparse. The performance of these 35 nodes on the validation set is RMSE = 0.226 Ω and 1.78°, only approximately 20% lower than the full model, but with approximately 62% less computational effort. The RBF parameters (center point, width, and weight) in the sparse node set are quantized using 8-bit fixed-point quantization to further compress storage space and generate a quantized parameter set. Overall, the data size of the compressed mapping model is reduced by about 85% compared to the original model, while the accuracy only decreases by about 20%.

[0082] 1.3.2. Quadtree or Octree Accelerated Index Construction. To achieve millisecond-level query response, an octree index structure is constructed. The input space is divided into multiple subspaces, with each leaf node containing a set of RBF nodes. The octree construction process: Determine the boundaries of the input space (ω, S, T) and construct a root node encompassing the entire space. Implement a recursive partitioning algorithm to evenly divide the space into eight subspaces. For each subspace, if the number of contained RBF nodes is greater than a threshold (set to 5) and the depth is less than a maximum (set to 4), continue partitioning. The resulting octree structure contains 127 nodes. For example, if the root node contains all 35 RBF nodes, after the first partition, the eight subspaces contain 7, 5, 8, 4, 3, 6, 2, and 0 nodes, respectively. For each leaf node of the octree, the contribution of the contained RBF nodes to the region is precomputed to generate a node influence table. Simultaneously, a spatial hash index, I_hash, is constructed to quickly map the input coordinates (ω, S, T) to leaf nodes. The result is an indexed mapping model that supports an average query response time of 0.5ms.

[0083] 1.3.3 Real-time impedance parameter retrieval and converter parameter calculation. In actual operation, given the current operating parameters (ωcurrent, Scurrent, Tcurrent), the index mapping model is used for rapid retrieval. For example, the current operating condition is ωcurrent = 2π•1000 rad / s, Scurrent = 850W / m 2 , Tcurrent=32°C, and locate the leaf node L15 through the index. The three most relevant RBF nodes are c1, c5, and c9, and the corresponding response value is: R1= w1•exp(-||x-c1|| 2 / (2σ1 2)) = [0.873, 0.651]•exp(-52.9 / 2922.1) = [0.819, 0.611];R5= [0.243, 0.387]•exp(-98.6 / 1681.0) = [0.209, 0.333];R9= [0.109, 0.075]•exp(-173.2 / 5476.5) = [0.104, 0.071]. Impedance estimate Z_est = R1+ R5+ R9= [1.132, 1.015] = |Z|_est∠φ_est = 4.479∠41.9° Ω = 3.344 + j2.991 Ω. Based on the estimated impedance value, the converter control parameter is calculated as: D = 0.3 +0.4•(1-e (-R_s / 5) ); K_p = 0.05 + 0.1•|Z|_est / 10; K_i = 20 + 30•atan(X_L / 5) / π; where D is the duty cycle of the DC-DC converter; R_s = Re(Z_est) is the real part of the impedance, which is 3.344 Ω; |Z|_est is the impedance modulus, which is 4.479 Ω; X_L = Im(Z_est) is the imaginary part of the impedance, which is 2.991 Ω; K_p and K_i are the proportional and integral coefficients of the PI controller, respectively. Substituting the numerical values ​​into the calculation: D = 0.3 + 0.4•(1-e -3.344 / 5 ) = 0.496; K_p = 0.05 + 0.1•4.479 / 10 = 0.095; K_i = 20 + 30•atan(2.991 / 5) / π = 25.64; finally, the converter control parameters C_converter = {D=0.496, K_p=0.095, K_i=25.64, R_s=3.344, X_L=2.991} are obtained.

[0084] Step 2: Based on the converter control parameters, pre-stored ambient temperature data and light intensity data, the temperature-sensing-free self-correction photovoltaic characteristics prediction is performed.

[0085] 2.1. Photovoltaic Thermal Balance Model and Parameter Identification. A photovoltaic module thermal balance model is established to describe the relationship between module temperature and environmental factors. Model parameters are identified using a sliding time window adaptive recursive least squares algorithm by analyzing historical operating data. In this example, parameter identification is performed using a dataset covering the last seven days, eight hours per day, and a 10-minute sampling interval (a total of 336 sample points). For example, at a certain moment, the following measurements are made: ambient temperature T_ambient = 25°C; light intensity S = 800W / m 2The photovoltaic output power (PPV) was 8000W. A recursive least squares algorithm was used to optimize the module temperature, minimizing the mean squared error between the measured temperature and the model-predicted temperature. The optimal parameter values ​​were calculated as follows: NOCT = 47.3°C; k = 0.00042°C / W; and the forgetting factor λ = 0.98 (dynamically adjusted based on the rate of environmental change). Substituting this into the thermal balance model, the module temperature was estimated as: T_module = 25 + (800 / 1000)•[47.3 - 20] + 0.00042•8000 = 50.2°C. The parameter confidence index τ = 0.87 (maximum score 1) was calculated based on the parameter change rate and data fitting error. The final output was the calibration parameter set θ = {k = 0.00042, NOCT = 47.3, λ = 0.98, τ = 0.87} and the combined thermal power model.

[0086] 2.2 Indirect temperature estimation based on electrical parameters.

[0087] 2.2.1. Partial Array Sampling Circuit Design. In this example, a 10kW PV array (40 modules) is divided into four subarray regions {SA_1, SA_2, SA_3, SA_4}, each containing 10 modules. A bidirectional switch and sampling circuit are configured to sample each subarray in rotation according to a schedule. The sampling control sequence is designed to have a 4-second period, sampling one subarray per second. The process is as follows: t = 0s: SA_1 is connected to the sampling circuit, and the other subarrays are connected to the main power line; t = 1s: SA_2 is connected to the sampling circuit, and the other subarrays are connected to the main power line; t = 2s: SA_3 is connected to the sampling circuit, and the other subarrays are connected to the main power line; t = 3s: SA_4 is connected to the sampling circuit, and the other subarrays are connected to the main power line. To ensure stable system output during the sampling process, a 1200μF energy buffer circuit is designed to maintain the main power line voltage fluctuation below 2V during measurement.

[0088] 2.2.2 Extracting Equivalent Open-Circuit Voltage and Short-Circuit Current. For each subarray, the following measurements were performed: the subarray's load was quickly disconnected, and the stabilized voltage was measured with a sampling time of 100ms. The subarray was shorted with a small resistor (0.1 Ω), and the current was measured with a sampling time of 100ms. For example, the measurement results for subarray SA_1 were: original open-circuit voltage: 43.2V; actual open-circuit voltage after compensation: VOC_1 = 43.8V; original short-circuit current: 8.5A; actual short-circuit current after compensation: ISC_1 = 8.7A. A statistical filtering algorithm was applied to all subarray measurement results to generate a filtered sampling data set. Then, considering the mismatch effect between sub-arrays, the equivalent open-circuit voltage and equivalent short-circuit current are calculated: VOC_eq = Σ(VOC_i•wi) / Σwi; ISC_eq = Σ(ISC_i•wi) / Σwi; where VOC_eq is the equivalent open-circuit voltage; ISC_eq is the equivalent short-circuit current; VOC_i and ISC_i are the open-circuit voltage and short-circuit current of the i-th sub-array, respectively; wi is the weight coefficient, which is related to the quality of the sub-array and is {1.0, 0.98, 0.97, 0.99} in this example. Substituting the numerical values ​​into the calculation: VOC_eq=(43.8×1.0+44.0×0.98+43.5×0.97+ 43.9×0.99) / (1.0+0.98+0.97+0.99)= 44.43V; ISC_eq = (8.7×1.0+ 8.6×0.98 + 8.5×0.97 + 8.7×0.99) / (1.0+0.98+0.97+0.99)= 8.68A.

[0089] 2.2.3. Calculate the temperature estimates for the two paths based on the equivalent open-circuit voltage, equivalent short-circuit current, and current light intensity. Substituting the values ​​into the equation: T_est_V = 25 + (44.43 - 45.6) / (-0.15) = 32.8°C; T_est_I = 25 + (8.68 / 800 - 8.9 / 1000) / (0.0036 / 1000) = 52.8°C. The temperature estimates for the two paths differ significantly, requiring confidence-weighted fusion. Calculate the weight coefficient based on the parameter confidence index τ = 0.87. Substituting the numerical values ​​into the equations, we obtain: W_V(τ) = exp(-|0.87-0.9| / 0.15) / (exp(-|0.87-0.9| / 0.15) + exp(-|0.87-0.7| / 0.15)) = 0.748; W_I(τ) = 1 - 0.748 = 0.252. Using these weights, we calculate the weighted fused temperature: T_module_est = 0.748 × 32.8 + 0.252 × 52.8 = 37.84°C. Finally, applying a Kalman filter and combining it with the historical temperature estimates yields a smoothed fused temperature estimate, T_module_est = 38.2°C, with a confidence index, CI_T, of 0.83.

[0090] 2.3. Self-correcting PV characteristic prediction and dynamic MPPT.

[0091] 2.3.1 Multi-condition PV curve prediction model. Based on the estimated module temperature T_module_est = 38.2°C and the current light intensity S = 800W / m 2 , calculate the operating parameters of the photovoltaic module. Substitute the numerical calculation: I_L=(800 / 1000)•[9.0+0.0006×(38.2-25)]=7.21A; I_0=1.2×10 -10 •(311.35 / 298.15) 3 •exp[1.12 / (8.617×10 -5 ×298.15)-1.12 / (8.617×10 -5 ×311.35)] = 4.78×10 -10A; R_s = 0.35•[1 + 0.001×(38.2 - 25)] = 0.355Ω; R_sh = 200•[1 + (-0.0005)×(38.2 - 25)] = 198.68Ω. Non-uniform sampling is used to generate 50 voltage points {V_i} within the voltage range [0, 1.2•VOC_est], with a higher density near the MPP (0.7•VOC_est to 0.85•VOC_est). For each voltage point V_i, the implicit photovoltaic equation is solved. Because I appears on both sides of the equation, the Newton-Raphson method is used for iterative solution. For the voltage point V = 36.5V, the iterative process is as follows: initial guess value I_0 = I_L = 7.21A; calculate f(I_0) = I_0 - I_L + I_0•[exp((V+I_0•R_s) / (n•V_t)) - 1]+ (V+I_0•R_s) / R_sh = 0.873A; calculate f'(I_0) = 1 + I_0•exp((V+I_0•R_s) / (n•V_t))•R_s / (n•V_t) + R_s / R_sh = 1.612; update I_1 = I_0 - f(I_0) / f'(I_0) =7.21 - 0.873 / 1.612 = 6.668A; repeat until convergence, and finally obtain I = 6.653A. This method calculates the current corresponding to each voltage point and, in turn, the power, P_i = V_i • I_i. For example, for V = 36.5V, P = 36.5 × 6.653 = 242.83W. By scanning all voltage points, the maximum power point is found: (VMPP, PMPP) = (35.8V, 243.29W). For the entire array (40 modules), the total maximum power point is: (VMPP_array, PMPP_array) = (35.8V, 9731.6W). Due to the uncertainty of the temperature estimate, a confidence interval for the PV curve is generated based on the confidence index CI_T = 0.83. At the MPP, the power confidence interval is [PMPP_min, PMPP_max] = [9512.4W, 9956.3W].

[0092] 2.3.2, Adaptive error correction mechanism. In order to improve the prediction accuracy, the historical prediction error database E_history is selected to match the current working condition (S = 800W / m 2, T_module_est = 38.2°C), forming a similar operating condition set S_similar. The prediction deviation statistics of these 10 records were calculated. By analyzing the data from this similar operating condition set, we obtained: μ_V = -0.0183 (average voltage overestimation of 1.83%); σ_V = 0.0092 (relative standard deviation of voltage prediction of 0.92%); μ_P = -0.0265 (average power overestimation of 2.65%); and σ_P = 0.0153 (relative standard deviation of power prediction of 1.53%). Based on these statistical characteristics, an adaptive Bayesian correction function was applied to adjust the current prediction. Substituting these into numerical calculations, the following results were obtained: VMPP_corrected = 35.8•(1 + (-0.0183)) = 35.14V; PMPP_corrected = 9731.6•(1 + (-0.0265)) = 9473.7W. At the same time, the confidence interval of the corrected MPP is calculated: VMPP_CI = [35.14•(1-1.96•0.0092), 35.14•(1+1.96•0.0092)] = [34.5V, 35.78V]; PMPP_CI = [9473.7•(1-1.96•0.0153), 9473.7•(1+1.96•0.0153)] = [9183.5W, 9763.9W].

[0093] 2.3.3 MPPT Reference Trajectory and Power Forecast Generation. Based on the corrected maximum power point parameters (VMPP_corrected = 35.14V, PMPP_corrected = 9473.7W) and the current system operating point (V_current = 34.2V, P_current = 9126.5W), design the optimal trajectory from the current point to the target MPP. System operating constraints include: maximum voltage change rate (dV / dt)_max = 2V / s; maximum power change rate (dP / dt)_max = 500W / s; and voltage range [V_min, V_max] = [30V, 45V]. A trajectory planner based on model predictive control (MPC) is used to design the optimization problem. Solving the MPC problem yields an optimal trajectory from the current operating point to the target MPP, consisting of 10 points: t=0.0s: (34.20V, 9126.5W); t=0.1s: (34.39V, 9176.3W); t=0.2s: (34.58V, 9225.7W); ...; t=0.9s: (35.14V, 9473.7W). This trajectory forms the MPPT reference trajectory R_MPPT = {(V_0=34.20V, P_0=9126.5W), (V_1=34.39V, P_1=9176.3W), ..., (V_9=35.14V, P_9=9473.7W)}. Simultaneously, based on the system's dynamic characteristics and the PV curve, the power trend over the next 60 seconds is predicted to generate a short-term power forecast trajectory. Taking into account possible changes in sunlight (according to the weather forecast, there may be short-term cloud cover), the power prediction trajectory shows that there may be a power drop of about 5% during t = 30s-40s.

[0094] Step 3: Perform multi-timescale data coordination and processing to generate the final execution instruction set.

[0095] 3.1. Electrical Fast Response Data Processing (Millisecond Level). Collect high-frequency electrical signals from the photovoltaic system at a sampling rate of 10kHz, including: photovoltaic output voltage VPV; photovoltaic output current IPV; electrolyzer voltage VEL; and electrolyzer current IEL. Read the MPPT reference track R_MPPT and calculate the real-time tracking deviation. For example, at t = 0.5s: Reference value: (V_ref = 34.77V, P_ref = 9274.8W); Actual value: (V_actual = 34.65V, P_actual = 9245.3W); Tracking deviation: e_track = (V_ref - V_actual) / V_ref = 0.34%. Apply a wavelet transform denoising algorithm to the raw sampled data, using the db4 wavelet basis and a 4-level decomposition, to increase the signal-to-noise ratio from the original 15dB to 25dB, resulting in a filtered electrical signal. Combining converter control parameters and tracking error, the real-time power is calculated as follows: PPV = VPV•IPV = 34.65V•267A = 9255.6W; PEL = VEL•IEL = 62.3V•142A = 8846.6W; system efficiency = PEL / PPV = 95.6%. The power change rate dP / dt is monitored, and a rapid response mechanism is triggered when a sudden event (such as dP / dt exceeding the threshold of 200W / s) is detected. For example, a power change rate of -220W / s detected at t=0.621s triggers a rapid response, reducing the converter duty cycle D from 0.496 to 0.483 to stabilize the system. Generates electrical fast-response control instructions C_fast, including: converter duty cycle D = 0.483; PI controller parameters K_p = 0.095, K_i = 25.64; control period T_c = 0.1ms; response priority P_f = "high." Also records event characteristics: event type: "power drop"; event time: t = 0.621s; event amplitude: dP / dt = -220W / s; trigger source: "light fluctuation." C_fast is transmitted to the DC-DC converter control unit, and the event characteristics are uploaded to the medium-speed data processing unit.

[0096] 3.2. Medium-speed pressure regulation data processing (seconds). Collect medium-frequency data of the electrolysis system with a sampling rate of 10Hz, including: electrolytic cell pressure P_EL = 1.5MPa; hydrogen production H2_flow = 0.58Nm 3 / h; system power P_system = 8846.6W. Receive event characteristics and the short-term power prediction trajectory P_predicted to assess the impact of energy fluctuations. According to the power prediction, a 5% power drop is expected between t = 30s and 40s, requiring advance adjustment of the pressure setpoint. The electrolyzer pressure-efficiency characteristic curve set P_efficiency is read from the system parameter library, which contains power-efficiency functions at five pressure levels: P_1 = 0.5MPa: η_1(P) = 68 - 0.2•P / 1000 + 0.01•(P / 1000) 2 ;P_2 = 1.0MPa: η_2(P) = 70 -0.25•P / 1000 + 0.015•(P / 1000) 2 ;P_3 = 1.5MPa: η_3(P) = 72 - 0.3•P / 1000 + 0.02•(P / 1000) 2 ;P_4 = 2.0MPa: η_4(P) = 71 - 0.35•P / 1000 + 0.025•(P / 1000) 2 ;P_5 =2.5MPa: η_5(P) = 69 - 0.4•P / 1000 + 0.03•(P / 1000) 2Where P is power (W) and η is efficiency (%). A dynamic programming algorithm was applied to optimize the pressure setpoints: Under the current operating conditions (P_system = 8846.6W), the efficiencies at each pressure level are: P_1 = 66.3%; P_2 = 68.2%; P_3 = 70.1%; P_4 = 69.0%; and P_5 = 66.8%. Considering that the power forecast shows a 5% power drop between t = 30s and 40s, the power is expected to drop to approximately 8404W. At this power level, the highest efficiency, 69.0%, is achieved at pressure level P_2 = 1.0MPa. Balancing the efficiency improvement with the pressure conversion cost, the optimization algorithm determined the optimal pressure setpoint, P_EL_set = 1.0MPa, requiring a downward adjustment from the current pressure of 1.5MPa. Construct a pressure regulation trajectory planner, taking into account the dynamic characteristics of the electrolysis system: maximum pressure change rate (dP / dt)_max = 0.1MPa / s; maximum pressure overshoot ΔP_max = 0.05MPa; generate a smooth transition path from the current pressure P_EL = 1.5MPa to the target pressure P_EL_set = 1.0MPa, using an S-shaped curve to ensure a smooth transition: t=0s: P = 1.50MPa; t=1s: P = 1.45MPa; t=2s: P = 1.35MPa; t=3s: P = 1.25MPa; t=4s: P = 1.15MPa; t=5s: P = 1.05MPa; t=6s: P = 1.00MPa; generate a pressure regulation control command C_pressure, including: target pressure P_target = 1.0 MPa; pressure trajectory {(t=0s, P=1.50 MPa), (t=1s, P=1.45 MPa), ..., (t=6s, P=1.00 MPa)}; control parameters {K_p=2.5, K_i=0.8, K_d=0.3}; control period T_c = 100 ms. Simultaneously, the regulation status S_pressure is recorded: current pressure P_current = 1.5 MPa; target pressure P_target = 1.0 MPa; expected completion time t_expected = 6 s; regulation progress = 0%; expected efficiency improvement Δη = 69.0% - 70.1% = -1.1% (slight short-term efficiency decrease, but overall efficiency will increase considering the upcoming power fluctuation). The pressure regulation control command is sent to the electrolysis system pressure control unit, and the regulation status is uploaded to the slow data processing unit.

[0097] 3.3 Slow Temperature Optimization Data Processing (Minute-Level). Low-frequency data from the electrolysis system is collected at a sampling rate of 0.1 Hz, including: electrolytic cell temperature T_EL = 55°C; system efficiency η_system = 70.1%; and environmental condition data ENV = {T_ambient = 25°C, humidity = 45%, cloud_cover = 20%}. Regulation status and historical operating data are received to assess long-term operating trends. Weather forecast data is obtained through an external meteorological service interface, indicating that light intensity will gradually decrease by approximately 20% over the next four hours. Combined with historical light patterns, a long short-term memory (LSTM) model is applied to generate energy forecast data. The training dataset consists of 4,320 samples sampled every 10 minutes over the past 30 days. After training, the LSTM model achieves a root mean square error (RMSE) of 52 W on the test set, with a relative error of approximately 5.7%. The generated energy forecast shows that the system power will show a downward trend in the next 4 hours: t=0h: P = 8846.6W; t=1h: P = 8500W; t=2h: P = 8100W; t=3h: P = 7600W; t=4h: P = 7200W. The temperature-efficiency characteristic data set D_T_η is obtained from the system parameter library, which contains the efficiency functions at 7 temperature levels: T_1 = 40°C: η_1(P, P_EL) = 65 + 0.1•P / 1000 - 0.05•P_EL; T_2 = 45°C: η_2(P, P_EL) = 67 + 0.12•P / 1000 - 0.06•P_EL; T_3 = 50°C: η_3(P, P_EL) = 69 +0.15•P / 1000 - 0.07•P_EL; T_4 = 55°C: η_4(P, P_EL) = 70 + 0.17•P / 1000 - 0.08•P_EL; T_5 = 60°C: η_5(P, P_EL) = 71 + 0.2•P / 1000 - 0.1•P_EL; T_6 = 65°C: η_6(P, P_EL) = 70 + 0.18•P / 1000 - 0.12•P_EL; T_7 = 70°C: η_7(P, P_EL) = 68 + 0.15•P / 1000 - 0.15•P_EL, where P is power (W), P_EL is the electrolyzer pressure (MPa), and η is efficiency (%). A temperature dynamics model is established, taking into account the thermal inertia of the electrolyzer. A temperature decision-making algorithm based on a receding horizon optimization is designed to optimize the total energy efficiency over the next four hours. Given that power forecasts indicate a power decrease over the next few hours, and that lower temperatures are often more efficient at lower power, the optimization algorithm determines the long-term optimal temperature setpoint, T_EL_set = 50°C, which needs to be adjusted downward from the current temperature of 55°C.A multi-step forward temperature control strategy was constructed. Based on the thermodynamic model, the required heating power was calculated. To reduce the temperature from 55°C to 50°C, a slow cooling strategy (heating is stopped and natural heat dissipation is utilized) was selected. The required time was calculated: heat dissipation loss Q_loss = k•(T-T_ambient) = 10•(55-25) = 300W; temperature change rate dT / dt = -Q_loss / C_thermal = -300 / 12000 = -0.025°C / s = -1.5°C / min; and required time = (50-55) / (-1.5) = 3.33 minutes. Taking into account the system's thermal inertia, an adaptive model predictive controller (AMPC) was used to generate the temperature control sequence. Generate temperature optimization control instructions, including: target temperature T_target = 50°C; temperature trajectory {(t=0min, T=55.0°C), (t=1min, T=53.5°C), (t=2min, T=52.0°C), (t=3min, T=50.5°C), (t=4min, T=50.0°C)}; heating power sequence {(t=0min, P=0W), (t=1min, P=0W), (t=2min, P=0W), (t=3min, P=0W), (t=4min, P=250W)}; control parameters {K_p=5.0, K_i=0.2, K_d=1.0}; control period T_c = 10s. At the same time, the temperature control status data S_temp is recorded: current temperature T_current = 55°C; target temperature T_target = 50°C; expected completion time t_expected = 4 minutes; adjustment progress progress = 0%; and expected efficiency improvement Δη = 1.2% (under the expected power curve). The temperature optimization control command is sent to the electrolysis system temperature control unit, and the temperature control status data is saved for cross-timescale coordination.

[0098] 3.4. Data coordination and conflict resolution across time scales.

[0099] 3.4.1. Multi-timescale data buffering mechanism. Design a three-layer data buffer structure and implement a cross-buffer data synchronization mechanism to ensure that upper-level decisions are passed to lower-level execution. For example, when a pressure setpoint change is detected, the corresponding electrical control parameters are updated to adapt to the new pressure conditions. Implement a data time tag and synchronization mechanism to solve the time alignment problem of data at different time scales. For example, high-frequency data at t = 2.345s is associated with the most recent medium-frequency data (t = 2.3s) and low-frequency data (t = 0s).

[0100] 3.4.2. Dynamic Priority Allocation Algorithm. Define the system state space S and the set of objective functions, including: f_1: efficiency maximization objective, current value = 70.1%; f_2: safety assurance objective, current value = 95% (safety margin); f_3: equipment lifespan protection objective, current value = 98% (healthiness). Based on the current system state, evaluate the importance and conflict risk of each objective. In the current state, calculate the importance of each objective: I_1 = 0.6•0.701 + 0.4•0.15 = 0.4806 (efficiency maximization objective); I_2 = 0.8•0.95 + 0.4•0.05 = 0.78 (safety assurance objective); I_3 = 0.5•0.98 + 0.4•0.10 = 0.53 (equipment lifespan protection objective). Based on the importance assessment, implement an adaptive priority adjustment strategy. In the current state, the priority of each objective is calculated: P_1 = exp(0.4806 / 0.2) / Z = 0.414; P_2 = exp(0.78 / 0.2) / Z = 0.534; P_3 = exp(0.53 / 0.2) / Z = 0.053; where Z = exp(2.403) + exp(3.9) + exp(2.65) = 26.67. The system identifies the current operating condition as closest to "normal" and generates a priority assignment rule: {Safety > Efficiency > Stability}. Based on this rule, the system prioritizes safety constraints, followed by efficiency, and finally stability.

[0101] 3.4.3. Command Conflict Detection and Coordination Algorithm. This algorithm uses a multi-scale buffer dataset, priority assignment rules, and target importance assessment to detect potential conflicts between control commands at the three timescales. In the current state, the following conflict is detected: A pressure-temperature setpoint conflict: The pressure regulation command requires a pressure reduction of 1.0 MPa, while the temperature optimization command requires a temperature reduction of 50°C. Simultaneously reducing both pressure and temperature could cause system efficiency to temporarily drop beyond a safety threshold. For this conflict, a conflict resolution space A_pressure_temp is constructed, with the values ​​A_pressure_temp = {a1: "Pressure first, then temperature", a2: "Temperature first, then pressure", a3: "Adjust simultaneously but reduce the rate of change", and a4: "Maintain status quo"}. Evaluate the utility of each option: Calculate the utility of each option: U(a1) = 0.414•0.8 + 0.534•0.9 + 0.053•0.7 = 0.8489; U(a2) = 0.414•0.7 + 0.534•0.85 + 0.053•0.8 = 0.7861; U(a3) = 0.414•0.85 + 0.534•0.95 + 0.053•0.9 = 0.9069; U(a4) = 0.414•0.6 + 0.534•0.75 + 0.053•0.95 = 0.6993. Based on the utility evaluation, select option a3: "simultaneous adjustment but reduced rate of change." Implementation: Reduce the pressure change rate by 50%, from 0.1 MPa / s to 0.05 MPa / s; reduce the temperature change rate by 30%, from 1.5°C / min to 1.05°C / min. Electrical control-pressure regulation conflict: During pressure regulation, electrical control commands require corresponding adjustments to maintain system stability. Solution: Modify electrical control parameters and establish a feedforward compensation mechanism to pre-adjust electrical control parameters based on pressure changes. Specific adjustments: Adjust the converter duty cycle D from 0.483 to 0.491; increase the PI controller parameter K_p from 0.095 to 0.115; and reduce the PI controller parameter K_i from 25.64 to 23.12. After conflict resolution, a coordinated control instruction set is generated, which includes modified control instructions at all levels: modified electrical fast response control instruction C_fast_mod: converter duty cycle D = 0.491; PI controller parameters K_p = 0.115, K_i = 23.12; control period T_c = 0.1ms; response priority P_f = "high".Modified pressure regulation control command C_pressure_mod: target pressure P_target = 1.0 MPa; pressure change rate dP / dt = 0.05 MPa / s; pressure trajectory {(t=0s, P=1.50 MPa), (t=2s, P=1.45 MPa), ..., (t=12s, P=1.00 MPa)}; control parameters {K_p=2.2, K_i=0.7, K_d=0.25}; control period T_c = 100 ms. The modified temperature optimization control instruction C_temperature_mod specifies: target temperature T_target = 50°C; temperature change rate dT / dt = 1.05°C / min; temperature trajectory {(t=0 min, T=55.0°C), (t=1.43 min, T=53.5°C), ..., (t=5.71 min, T=50.0°C)}; the heating power sequence is appropriately adjusted; control parameters {K_p=4.8, K_i=0.18, K_d=0.9}; and control period T_c = 10s. A control translator is constructed to convert the coordinated abstract instructions into control signals for specific execution units. For example, the converter control signal is a PWM waveform with a frequency of 20 kHz and a duty cycle of 0.491; the pressure control signal is a 38% electric valve opening and a 65% compressor speed; and the temperature control signal is a 0W heater power (in the cooling phase) and a 20% cooling fan speed. Output the final execution instruction set C_final, which contains all modified instructions and control signals for system execution and reinforcement learning training.

[0102] Step 4: Based on the final execution instruction set, the constructed multi-domain coupling model and two-layer deep reinforcement learning network are used to generate system comprehensive control instructions.

[0103] 4.1. Construction of multi-domain coupled environment model.

[0104] 4.1.1. Physical Model Integration and Interface Design. Integrate the models constructed in the previous steps: the impedance-frequency characteristic mapping model; the combined thermoelectric model M_thermo_elec; and the electrolysis temperature-pressure characteristic model M_EL_PT. Build a modular model integration framework and define a unified data exchange interface: Photovoltaic model module PV_Model: Inputs: jω (angular frequency), S (light intensity), T (temperature); Outputs: ZPV (impedance), VOC (open-circuit voltage), ISC (short-circuit current), VMPP (maximum power point voltage), PMPP (maximum power point power); DC-DC converter model module DC_DC_Model: Inputs: ZPV (photovoltaic impedance), D (duty cycle), L (inductance), C (capacitance); Outputs: ZDC (output impedance), η_converter (conversion efficiency). Electrolyzer_Model: Inputs: VDC (DC voltage), IDC (DC current), P (pressure), T (temperature); Outputs: H2_flow (hydrogen production), η_EL (electrolysis efficiency), dT / dt (temperature change rate). This module implements inter-model data flow pipelines. For example, given the input conditions: ω = 2π•1000 rad / s, S = 800 W / m 2 , T = 38.2°C, the calculation process is: PV_out = PV_Model(2π•1000, 800, 38.2) = {ZPV=4.479∠41.9°Ω, VOC=44.43V, ISC=8.68A, VMPP=35.14V, PMPP=9473.7W}; DC_in = DataAdapter(PV_out) = {ZPV=4.479∠41.9°Ω, VPV=35.14V, IPV=269.6A}; DC_out = DC_DC_Model(DC_in, D=0.491, L=0.5mH, C=4700μF) = {ZDC=0.18+j0.05Ω, η_converter=0.969}; EL_in = DataAdapter(DC_out) = {VDC=68.9V, IDC=133.1A, P=1.5MPa, T=55°C}; EL_out = Electrolyzer_Model(EL_in) = {H2_flow=0.58Nm 3 / h, η_EL = 0.723, dT / dt = 0.002°C / s}. All sub-models are integrated to construct a complete system multi-physics model, SysModel, capturing the coupled relationships among the photovoltaic, power electronics, and electrolysis domains. The model solver employs a hybrid solution strategy: an explicit Runge-Kutta method with a 0.1ms step size for electrical fast dynamics; an implicit Backward Euler method with a 10ms step size for electrolysis chemical reactions; and an adaptive step size method with a step size range of 0.1s-10s for thermodynamic processes.

[0105] 4.1.2. Adaptive Parameter Identification Mechanism. A recursive extended Kalman filter (REKF) is constructed to identify model parameters online. Electrical parameters (L, C, R): High-frequency identification (1kHz), for example, the output inductance of the DC-DC converter is identified from the initial value L = 0.5mH to L = 0.485mH; electrolytic parameters (k0, α, Ea): Medium-frequency identification (1Hz), for example, the exchange current density of the electrolytic cell is identified from the initial value k0 = 0.1A / cm 2 Identified as k0=0.094A / cm 2 Thermal parameters (Cp, h): Low-frequency identification (0.01 Hz). For example, the electrolytic cell's heat dissipation coefficient is identified from its initial value of h = 10 W / °C to h = 10.38 W / °C. Update the parameters in the system multiphysics model SysModel to obtain the calibrated multidomain coupling model MCM_calibrated.

[0106] 4.1.3. Environment model encapsulation and reinforcement learning interface implementation. Based on a model credibility index (MCI) of 0.89 (maximum score: 1), adaptive adjustments are made to different submodules. For example, random perturbations are added to certain parameters in the electrolysis model with high uncertainty to enhance the robustness of the policy. An OpenAI Gym-compatible environment interface is designed to encapsulate the system model. Definitions: Observation space: Box(low=[-1, -1, ...], high=[1, 1, ...], shape=(20,)) containing 20 normalized state variables; Action space: Box(low=[-1, -1, ...], high=[1, 1, ...], shape=(8,)) containing 8 normalized control variables; Implement the step method to calculate state transitions and rewards; Build a differentiable reward function calculator to support multi-objective reinforcement learning. For example, for the state transition (state, action, next_state), the reward of each component is calculated as: system_efficiency(next_state) = 0.723 (electrolysis efficiency) × 0.969 (conversion efficiency) = 0.701; stability_measure(state, next_state) = exp(-10•||next_state - state|| 2 ) = 0.904; safety_margin(next_state) = min(1, distance_to_boundary(next_state) / safe_distance) = 0.85; total reward = 0.6 × 0.701 + 0.3 × 0.904 + 0.5 × 0.85 = 1.1168. Develop an environment randomization mechanism to enhance the generalization ability of the policy. For example, at the beginning of each episode, randomize parameters such as light intensity and ambient temperature within a range of ±5%. Output a reinforcement learning environment interface for subsequent policy network and execution network training.

[0107] 4.2. Policy network training and global energy optimization.

[0108] 4.2.1. Hierarchical Reward Function Design. Define a state space consisting of 11 global state variables, such as light intensity, ambient temperature, and system efficiency. Define an action space consisting of three high-level decision variables, such as power allocation ratio and operating mode. Design a hierarchical reward function structure, including a hierarchical reward function, an efficiency reward function, a stability reward function, and a safety reward function, and implement a dynamic weight adjustment mechanism. For example, w_safe(s) = 0.5 + 0.3•(1 - distance_to_boundary(s) / safe_distance) = 0.5 + 0.3•(1 - 0.85) = 0.545.

[0109] 4.2.2. Exploration Strategy Based on Posterior Sampling. Implement a policy network training process based on the PPO algorithm framework. The network structure is as follows: Input layer: 11 nodes, corresponding to global state variables; Hidden layer 1: 128 nodes, with ReLU activation function; Hidden layer 2: 64 nodes, with ReLU activation function; Output layer: Gaussian distribution with mean μ and standard deviation σ, representing each action dimension. Construct an exploration strategy based on posterior sampling to maintain the posterior distribution of the environment dynamics model. For example, sample five model instances from the posterior distribution, corresponding to different electrolyzer performance curves, and train the optimal sub-policy π_i for each model. Calculate the posterior probabilities for each model: M_1: w_1 = 0.35; M_2: w_2 = 0.25; M_3: w_3 = 0.20; M_4: w_4 = 0.15; M_5: w_5 = 0.05. For state s, the action distribution parameters output by each sub-strategy are: π_1(a|s): μ_1 = [0.3, 0.7, -0.2], σ_1 = [0.1, 0.1, 0.1]; π_2(a|s): μ_2 = [0.4, 0.6, -0.3], σ_2 = [0.1, 0.1, 0.1]; π_3(a|s): μ_3 = [0.35, 0.65, -0.25], σ_3 = [0.1, 0.1, 0.1]; π_4(a|s): μ_4 = [0.25, 0.75, -0.15], σ_4 = [0.1, 0.1, 0.1]; π_5(a|s): μ_5 = [0.45, 0.55, -0.35], σ_5 = [0.1, 0.1, 0.1]. The parameters of the exploration strategy distribution after mixing are: μ_explore = 0.35•[0.3, 0.7, -0.2]+ 0.25•[0.4, 0.6, -0.3] + 0.20•[0.35, 0.65, -0.25]+ 0.15•[0.25, 0.75, -0.15] + 0.05•[0.45, 0.55, -0.35]= [0.33, 0.665, -0.235]; σ_explore = sqrt(0.35•[0.1 2 + (0.3-0.33) 2 ] + ... + 0.05•[0.1 2 + (0.45-0.33) 2 ]) = [0.126, 0.131, 0.128]. The final hybrid exploration strategy π_explore(a|s) = N(μ_explore, σ_explore 2) = N([0.33, 0.665, -0.235], [0.126 2 , 0.131 2 , 0.128 2 ]). Implement progressive exploration - using a balancing mechanism to reduce the exploration ratio as training progresses.

[0110] 4.2.3. Implementation of a Neural-Symbolic Hybrid Architecture. Build a neural-symbolic hybrid architecture that combines connectionist and symbolic approaches. Examples of symbolic knowledge rules include: Rule 1: If light fluctuation is high, reduce the power allocation ratio; Rule 2: If the electrolyzer temperature is greater than 65°C, switch to low-temperature priority mode; Rule 3: If the system efficiency is less than 60%, recalculate the optimal operating point. Implement a neural network policy constraint mechanism and construct a differentiable logic rule encoder. For example, for rule 2: "If the electrolyzer temperature is greater than 65°C, switch to low temperature priority mode": f(x) = (T_EL - 65) / 10, which indicates the degree to which the temperature exceeds 65°C; when T_EL = 60°C, f(x) = (60-65) / 10 = -0.5, R(x) = sigmoid(10•(-0.5-0.5))≈ 0; when T_EL = 70°C, f(x) = (70-65) / 10 = 0.5, R(x) = sigmoid(10•(0.5-0.5)) = 0.5; when T_EL = 75°C, f(x) = (75-65) / 10 = 1.0, R(x) = sigmoid(10•(1.0-0.5))≈ 0.993. Implement the knowledge distillation process to inject expert knowledge into the neural network: For example, for the state s = {S=800W / m 2 , T_ambient=25°C, T_EL=70°C, ...}, the action distribution given by the expert strategy according to Rule 2 is π_expert(a|s) = N([0.2, 0.8, -0.4], [0.05 2 , 0.05 2 , 0.05 2 ]), indicating low temperature priority mode. The distribution of the current policy network output is π_θ(a|s) = N([0.3,0.7,-0.2],[0.1 2 , 0.1 2 , 0.1 2 By minimizing the KL divergence, the policy network is encouraged to learn expert knowledge. A dynamic node expansion mechanism is constructed to adaptively adjust the network structure based on task requirements. For example, when a large error is detected in a specific area, symbol nodes processing that area are added. After training, a knowledge-enhanced policy network (GPN+) is formed, combining data-driven learning with prior physical knowledge.

[0111] 4.3. Perform network training and local parameter optimization.

[0112] 4.3.1. Local Parameter Space Design. Read the global decision D_global = {ratio_power = 0.33, mode_selection = 0.665, stability_preference = -0.235} output from the optimization strategy network. Define the local state space S_local, which contains 16 device-level state variables: photovoltaic voltage VPV = 35.14V; photovoltaic current IPV = 269.6A; DC-DC output voltage VDC = 68.9V; DC-DC output current IDC = 133.1A; electrolyzer voltage VEL = 68.9V; electrolyzer current IEL = 133.1A; electrolyzer pressure PEL = 1.5MPa; electrolyzer temperature TEL = 55°C; and eight other state variables. Define a local action space A_local, containing eight specific control parameters: converter duty cycle D_converter = 0.491; switching frequency f_switching = 20 kHz; voltage controller proportional coefficient K_p_V = 0.115; voltage controller integral coefficient K_i_V = 23.12; current controller proportional coefficient K_p_I = 0.08; current controller integral coefficient K_i_I = 15.5; pressure valve opening valve_pressure = 0.38; and temperature control signal temp_control = 0.2. Establish a mapping between global decisions and local parameters. For example, in the global decision, ratio_power = 0.33 indicates that the ratio of electrolysis power to photovoltaic power should be 33%, which constrains the converter duty cycle D_converter to the range [0.47, 0.51]. Implement parameter normalization and scaling mechanisms. For example, for the converter duty cycle D_converter: Original range: [0.1, 0.9]; Current value: 0.491; Normalized value: 2•(0.491-0.1) / (0.9-0.1) - 1 = -0.022. Construct a parameter correlation analyzer to calculate the coupling relationship between parameters: C_ij = correlation(a_i, a_j) = E[(a_i - μ_i)(a_j - μ_j)] / (σ_i•σ_j); where C_ij is the correlation coefficient between parameters a_i and a_j; E[•] is the expected value; μ_i and μ_j are the means of parameters a_i and a_j, respectively; and σ_i and σ_j are the standard deviations of parameters a_i and a_j, respectively.By analyzing historical data, we derived the parameter correlation matrix C. For example, C_D_converter, K_p_V = 0.72 (strong correlation); C_D_converter, f_switching = 0.15 (weak correlation); and C_valve_pressure, temp_control = 0.68 (strong correlation). Based on this correlation matrix, we grouped and decoupled the eight parameters, dividing them into three groups: G_1 = {D_converter, K_p_V, K_i_V} (converter control parameter group); G_2 = {K_p_I, K_i_I, f_switching} (current control parameter group); and G_3 = {valve_pressure, temp_control} (electrolysis system parameter group). The output local parameter space, S_A_local, contains the state space, action space, and constraints, which are used for network training.

[0113] 4.3.2. Gradient-guided Exploration Strategy. Construct an execution network training framework based on the SAC algorithm. The network structure is as follows: policy network π_φ: a 3-layer MLP with hidden layer size [128, 64]; dual Q networks Q_θ1 and Q_θ2: 3-layer MLP with hidden layer size [128, 64]. Physical gradient information is extracted from the multi-domain coupling model MCM_optimal: for example, for the current state s = {VPV = 35.14V, IPV = 269.6A, ...} and action a = {D_converter = 0.491, K_p_V = 0.115, ...}, the model predicts the next state s' and reward R(s, a, s'), and then calculates the gradient of the reward with respect to each action component. For the G_1 group of parameters, the calculated physical gradients are: ΨR / ΨD_converter = 1.85; ΨR / ΨK_p_V = 0.73; ΨR / ΨK_i_V = -0.28, where Ψ is the partial derivative. Construct a gradient-guided exploration strategy: For example, for the policy network output a_mean = [0.491, 0.115, 23.12, ...], standard normal noise ε = [0.2, -0.3, 0.1, ...], physical gradient g_physics = [1.85, 0.73, -0.28, ...]: Normalize the physical gradient: normalize(g_physics) = [0.92, 0.36, -0.14, ...]; Calculate the exploration direction: ε + β•normalize(g_physics) = [0.2, -0.3, 0.1, ...] + 0.5•[0.92, 0.36, -0.14, ...] = [0.66, -0.12, 0.03, ...]; Generate exploration action: a_explore = [0.491, 0.115, 23.12, ...] + 0.2•[0.66, -0.12, 0.03, ...] = [0.623, 0.091, 23.126, ...]. Implements a gradient priority sampling mechanism: For example, for a state-action pair (s, a), if the physics gradient norm ||g_physics(s, a)|| = 2.0, its sampling probability will be exp(0.1•(2.0-1.0)) = exp(0.1) ≈ 1.105 times higher than that of a state-action pair with a gradient norm of 1.0. This allows the algorithm to explore high-gradient regions more frequently, accelerating the learning process.We built an adaptive exploration temperature adjustment mechanism. For example, if the current policy entropy H_t = -3.5 is higher than the target entropy H_target = -4, we update the temperature parameter: α_t+1 = α_t + 0.01•(-4 - (-3.5)) = α_t - 0.005, lowering the temperature parameter and reducing exploration. After training, we obtained the efficient execution network LEN+, which possesses gradient-guided exploration capabilities and achieves approximately 40% higher sample efficiency than the standard SAC algorithm.

[0114] 4.3.3. Continuous learning and adaptive adjustment mechanism. Build a sliding time window data buffer to maintain the latest 7 days of operating data, containing 8640 data points per day (sampling period 10 seconds). Implement an adaptive learning rate adjustment mechanism: For example, if the performance indicators of two consecutive evaluations are R_t-1 = 1.2 and R_t = 1.15, the adjusted learning rate is: η_t = 0.001•exp(-10•|1.15-1.2| / 1.2) = 0.000659. Build a multi-model integration architecture and maintain three execution network versions at the same time: LEN_1: Optimized for high light conditions (>700W / m 2 );LEN_2: Optimized for medium light conditions (400-700W / m 2 );LEN_3: Optimized for low light conditions (<400W / m 2 ). Implement the model selector: For example, the current light intensity is 800W / m 2 , calculate the similarity with each model condition: similarity(800W / m 2 , "high light") = 0.9; similarity(800W / m 2 , "medium light") = 0.4; similarity(800W / m 2 , "low light") = 0.1; select the LEN_1 model with the highest similarity. Build a knowledge transfer mechanism: For example, when the LEN_1 model learns new knowledge under high light conditions, it can partially transfer it to the LEN_2 model: θ_LEN_2 = θ_LEN_2 + 0.1•(θ_LEN_1 - θ_LEN_2). Implement an anomaly detection and recovery mechanism: For example, if the average performance over the past five moments is 1.3, and the current performance drops to 1.0, below the threshold of 0.8 × 1.3 = 1.04, network retraining is triggered. Build an incremental learning framework: Integrate all adaptive mechanisms and output an adaptive execution network.

[0115] 4.4. Double-layer network collaborative reasoning and system control execution. Read the current system state S_current = {S=800W / m 2, T_ambient=25°C, T_EL=55°C, P_EL=1.5MPa, ...}. The global decision is calculated using the optimization policy network GPN_optimal: D_global = forward_pass(GPN_optimal, S_current); where D_global is the global decision output; forward_pass is the forward propagation function of the neural network; GPN_optimal is the optimization policy network; and S_current is the current system state. In this example, the calculated global decision is D_global = {ratio_power=0.33, mode_selection=0.665, stability_preference=-0.235}. This indicates that 33% of the PV power should be allocated to the electrolysis system; the operating mode selection value of 0.665 corresponds to "efficiency priority mode"; and the stability preference value of -0.235 indicates a moderate tolerance for fluctuations in exchange for higher efficiency. Based on the global decision and the current local state S_local = {VPV = 35.14V, IPV = 269.6A, ...}, the adaptive execution network LEN_adaptive infers the specific control parameters: P_control = forward_pass(LEN_adaptive, S_local, D_global). P_control is the control parameter set; forward_pass is the forward propagation function of the neural network; LEN_adaptive is the adaptive execution network; S_local is the current local state; and D_global is the global decision. In this example, the calculated control parameter set P_control = {D_converter = 0.483, K_p_V = 0.098, K_i_V = 24.31, K_p_I = 0.075, K_i_I = 16.2, f_switching = 20kHz, valve_pressure = 0.42, temp_control = 0.15}. A feedforward-feedback hybrid control strategy is applied: u(t) =u_ff(t) + u_fb(t); u_ff(t) = f_ff(s(t), P_control); u_fb(t) = K•(r(t) - y(t)); where u(t) is the final control signal; u_ff(t) is the feedforward control component; u_fb(t) is the feedback control component; f_ff is the feedforward control function; s(t) is the system state; P_control is the control parameter set; K is the feedback gain; r(t) is the reference signal; y(t) is the system output.For the converter control signal, the feedforward component sets the duty cycle D = 0.483, and the feedback component is calculated and corrected based on the voltage error (the difference between the reference value of 35.14V and the actual value of 35.08V). The final control signal is D = 0.483 + 0.098 (35.14 - 35.08) / 35.14 = 0.48317. A real-time inference optimization accelerator was constructed, quantizing the two-layer network to INT8, reducing inference time from 15ms to 3.5ms, meeting the real-time control requirement of less than 5ms. The control instructions C_RL generated by the reinforcement learning network and the final execution instruction set C_final from S34 are weighted and fused: C_system = w_RL • C_RL + w_final • C_final, where C_system represents the system's integrated control instructions; w_RL represents the weights of the reinforcement learning control instructions; C_RL represents the weights of the reinforcement learning control instructions; w_final represents the weights of the final execution instruction set; and C_final represents the final execution instruction set. The weights are dynamically adjusted based on the system state. For example, the weight of C_final is increased during rapid transients, while the weight of C_RL is increased during steady-state conditions. In the current state, w_RL is set to 0.7 and w_final is set to 0.3, resulting in the fused system comprehensive control command. The final output is the comprehensive system control command C_system, which includes: converter control signal: PWM waveform, frequency 20kHz, duty cycle 0.483; electrolysis pressure control signal: electric valve opening 42%; electrolysis temperature control signal: cooling system power 15%. C_system is distributed to each execution unit to achieve coordinated control of the photovoltaic off-grid hydrogen production system. Execution feedback data D_feedback is also collected and stored in the system operation database DB_operation for model optimization and control adjustments in the next cycle.

[0116] By constructing a four-dimensional mapping relationship between light, temperature, frequency, and impedance, this method completely breaks free from the constraints of traditional equivalent circuit models and avoids the complex parameter identification process. It increases the data compression rate to 93.3% while maintaining a reconstruction error of less than 3%. It achieves accurate representation of complex nonlinear characteristics and real-time query response in milliseconds (average 0.5ms), which is approximately 10 times faster than traditional methods (typically 5-10ms). Under fluctuating light conditions, the system's response time to impedance characteristic changes is reduced from hundreds of milliseconds in traditional methods to less than 1 millisecond, effectively avoiding system instability caused by impedance mismatch and increasing the power output of the photovoltaic array by approximately 6-8%. Through indirect temperature estimation using electrical parameters, precise component temperature acquisition is achieved without the need for additional temperature sensors, with a measurement error within ±1.5°C, far superior to conventional irradiance models (typically ±5°C). By rotating sampling of each subarray, the system ensures continuous maximum power output. Dual estimation of the voltage and current paths effectively overcomes the limitations of single parameter estimation, improving temperature estimation accuracy by approximately 65% ​​and boosting MPPT tracking accuracy from 97% with traditional methods to 99.2%. Under rapidly changing partial shading conditions, the system recovery time is shortened from 11.5 seconds with traditional MPPT methods to 3.2 seconds, and power fluctuations are reduced from 38% to 14%, enhancing the system's adaptability and stability under complex lighting conditions. By constructing a three-tier control architecture consisting of millisecond-level electrical control, second-level pressure regulation, and minute-level temperature optimization, the system effectively addresses the significant differences in the dynamic characteristics of different physical processes, enabling the coordinated execution of control commands at different time scales and reducing the control command conflict rate from 35% with traditional methods to below 5%. Especially under rapidly changing system conditions, the system automatically adjusts the priority of each control objective based on the dynamic importance of safety, efficiency, and stability, ensuring the system maintains optimal operating conditions. Experiments show that this invention improves system efficiency by 9.2% and control stability by 25% during periods of fluctuating light intensity, while also reducing energy fluctuations by nearly 60%, enhancing the overall performance and reliability of the photovoltaic off-grid hydrogen production system. By incorporating electrolyzer pressure and temperature as optimization variables into the control system, dynamic matching of electrolyzer operating parameters with photovoltaic output is achieved. Based on short-term power forecast trajectories, electrolysis efficiency is improved by approximately 3% during power fluctuations. Combining long-term energy forecast data with rolling time-domain optimization, the system improves average efficiency by approximately 2.5% during varying daytime light intensity. Specifically, when an impending power drop (e.g., cloud cover) is detected, the system reduces pressure from 1.5 MPa to 1.0 MPa in advance, leaving sufficient margin for subsequent power recovery. When power is predicted to continue to decrease over the next few hours, the system gradually reduces temperature to the optimal low-power operating point (e.g., from 55°C to 50°C), avoiding energy waste. This results in an approximately 12% increase in hydrogen production during full-day operation, with particularly significant improvements on cloudy days when light intensity fluctuates frequently.By constructing a multi-domain coupled environment model and a two-layer deep reinforcement learning architecture, the limitations of traditional control methods have been overcome. The two-layer architecture achieves hierarchical optimization of global decision-making and local control, reducing learning complexity and increasing the convergence speed by about 3 times. Sample efficiency is improved, reducing the number of samples required for training by about 40%; the gradient information of the physical model is used to guide the exploration direction, further improving learning efficiency. Embedding physical knowledge as a priori constraints into the policy network improves the model's generalization ability by about 21%. Experimental results show that the present invention improves the system's efficiency by 4.3% under stable lighting conditions and 9.2% under lighting fluctuations. When shadows change rapidly, the system recovery time is shortened by 72%, the hydrogen production fluctuation is reduced by 60%, and the system's adaptability and robustness are enhanced.

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

Claims

1. A coordinated control method for a photovoltaic off-grid hydrogen production system, characterized in that: include: Collect the time domain signals of PV array voltage and current, build impedance frequency characteristic mapping model, and calculate converter control parameters; Collect ambient temperature data and light intensity data, perform temperature sensor-free self-correction photovoltaic characteristic prediction, generate MPPT reference trajectory and short-term power prediction trajectory, and combine them with converter control parameters to perform multi-time scale data coordination and processing to generate the final execution instruction set; Based on the final execution instruction set, the constructed multi-domain coupling model and two-layer deep reinforcement learning network are used to generate and issue system comprehensive control instructions.

2. The method according to claim 1, characterized in that Calculate converter control parameters, including: The time-domain signals of the PV array voltage and current under various light intensity and temperature conditions are collected and processed through FFT to obtain the impedance modulus and phase angle, forming the original impedance spectrum and performing tensor decomposition on it. Based on this, a radial basis function network is used to construct the four-dimensional mapping relationship between light, temperature, frequency and impedance, forming a global mapping model. The global mapping model is sparsely mapped and quantized, and an index structure is configured to obtain an indexed mapping model. The current operating parameters are used as input to obtain the current impedance parameters and calculate the converter control parameters accordingly.

3. The method according to claim 2, characterized in that Form a global mapping model, including: The original impedance spectrum is organized into a fifth-order tensor, and Tucker decomposition is applied to obtain the core tensor and factor matrix, generating a decomposed data tensor, from which the characteristic distribution map is extracted; Based on the feature distribution map, a three-layer radial basis function network structure is constructed and the initial RBF center point is determined to form an initial RBF network structure model and train it. The number of RBF nodes is dynamically adjusted during training, and a global mapping model is obtained through Bayesian hyperparameter optimization.

4. The method according to claim 2, characterized in that Obtain the current impedance parameters and calculate the converter control parameters based on them, including: Calculate the contribution of each RBF node in the global mapping model and arrange them in descending order. Filter the nodes whose contribution ranking exceeds the preset threshold as a sparse node set and quantize them at fixed points to obtain a quantization parameter set. Combine the sparse node set and the quantization parameter set to form a compressed mapping model. Based on the RBF center point and coverage radius in the compression mapping model, a tree index structure is constructed, the RBF node contribution is calculated for the leaf nodes to generate a node influence table, and a spatial hash index is constructed to form an indexed mapping model. The indexed mapping model is used to locate leaf nodes according to the current operating parameters, and the most relevant RBF nodes are extracted to calculate the current impedance parameters and generate the converter control parameters.

5. The method according to claim 1, wherein Generate MPPT reference trajectory and short-term power prediction trajectory, including: Collect ambient temperature data, light intensity data, and photovoltaic output power, establish a photovoltaic module thermal balance model and identify model parameters, and output a calibration parameter set and a combined thermal and power model; The equivalent open-circuit voltage and equivalent short-circuit current are obtained and combined with the calibration parameter set to calculate the fusion temperature estimate. This is combined with the current light intensity to predict the PV curve characteristics and perform error correction to generate the corrected maximum power point, forming the MPPT reference trajectory and short-term power prediction trajectory.

6. The method according to claim 5, characterized in that Calculates an estimate of the fusion temperature, including: Read the PV system topology data and sampling schedule, divide the PV array into predetermined sub-array areas, configure the bidirectional switch and sampling circuit, generate the sampling control sequence and sub-array configuration parameters, and execute the measurement program on the sub-arrays in sequence accordingly; During measurement, the original sampling data set is obtained and statistically filtered to form a filtered sampling data set. Taking into account the mismatch effect between sub-arrays, the equivalent open-circuit voltage and equivalent short-circuit current are calculated. The data are combined with the current light intensity and the calibration parameter set to calculate the temperature estimates of the voltage path and current path respectively. The data are then weighted, fused, and filtered to generate a fused temperature estimate.

7. The method according to claim 5, characterized in that Forming the MPPT reference trajectory and short-term power prediction trajectory, including: Calculate the operating parameters of the photovoltaic module based on the fusion temperature estimate and the current light intensity, solve the photovoltaic equation within the voltage range, and generate the PV curve characteristics; The maximum power point in the PV curve characteristics is searched, similar operating condition records are screened from the historical prediction error database, the statistical characteristics of the prediction deviation are calculated, and the corrected maximum power point is generated through correction function correction; it is combined with the pre-stored current system operating point and system operation constraints, and a trajectory planner based on model predictive control is used to solve the optimization problem considering the dynamic constraints of the system, generating the MPPT reference trajectory and short-term power prediction trajectory from the current operating point to the target maximum power point.

8. The method according to claim 1, characterized in that Generate the final execution instruction set, including: At the electrical control layer, the MPPT reference trajectory and converter control parameters are combined to monitor the power change rate and generate electrical fast-response control instructions and event characteristics; In the pressure control layer, based on event characteristics and short-term power prediction trajectory, the optimal pressure set point is calculated, and the pressure regulation control instructions and regulation status are generated; In the temperature control layer, the adjustment status is combined with the pre-stored historical operation data and weather forecast data to calculate the long-term optimal temperature set point and obtain the temperature optimization control instruction; The electrical fast response, pressure regulation and temperature optimization control instructions are stored in a multi-layer data buffer, and the priority allocation rules are used to generate the final execution instruction set.

9. The method according to claim 5, characterized in that Generate system comprehensive control instructions, including: Based on the impedance-frequency characteristic mapping model and the thermoelectric combined model, a photovoltaic-electrolysis multi-domain coupling model is constructed and encapsulated as a reinforcement learning environment interface. Based on this, a knowledge-enhanced policy network is trained to output an optimized policy network and global decision-making. Based on the global decision, define the local state space and action space, and train the adaptive execution network; Read the current system state, infer and generate a control parameter set by optimizing the strategy network and adaptive execution network, fuse it with the final execution instruction set, and output the system comprehensive control instructions.

10. The method according to claim 9, characterized in that By optimizing the strategy network and adaptive execution network, the system outputs comprehensive control instructions, including: The current system state is input into the optimization strategy network, and the global decision is calculated through forward propagation, including the power allocation ratio and working mode selection; The global decision and current local state are input into the adaptive execution network, which generates a specific control parameter set through forward calculation, including the converter duty cycle, controller parameters, and device operating parameters. Based on this, a feedforward-feedback hybrid control strategy is applied, combining model predictive control and feedback control, to generate reinforcement learning control instructions. The reinforcement learning control instructions are weighted and fused with the final execution instruction set. The weight of each instruction is dynamically adjusted according to the current system state, and the system comprehensive control instructions are output.

Images