Photovoltaic off-grid hydrogen production system operation mode switching instant DC bus voltage stability control method

Abstract

The invention discloses a photovoltaic off-grid hydrogen production system operation mode switching instant DC bus voltage stability control method, which comprises the steps of collecting a DC bus voltage signal, and reconstructing a high-dimensional phase space trajectory; calculating real geometric features such as phase space curvature and torsion of the trajectory; searching a historical similar trajectory segment in a historical trajectory database, and generating a final prediction trajectory at a future moment according to a known evolution path of the historical similar trajectory segment; and generating a feedforward control signal according to the deviation between the final prediction track and the reference voltage, and combining the feedforward control signal with traditional feedback control to form a final control instruction. The method does not depend on an accurate system mathematical model, and realizes rapid, smooth and stable control of voltage under mode switching by deeply excavating dynamic geometric information of voltage signals.

Description

Technical Field

[0001] This invention belongs to the field of power electronics and intelligent control, and in particular, it is a method for stabilizing the DC bus voltage during the instantaneous switching of the operation mode of a photovoltaic off-grid hydrogen production system. Background Technology

[0002] Off-grid photovoltaic hydrogen production systems, as an effective way to locally absorb unstable photovoltaic power and produce zero-carbon green hydrogen, have shown great application potential in remote areas, islands, and distributed energy scenarios. In such systems, the DC bus, as the energy collection and exchange hub connecting core components such as photovoltaic arrays, DC-DC converters, and electrolyzers, directly affects the efficiency of the entire system and the quality of hydrogen production. Especially during instantaneous switching of system operating modes (such as from standby to full-load startup, power step transitions, etc.), the severe power imbalance between the source and load ends can cause severe fluctuations or even instability in the DC bus voltage, potentially leading to device damage and system collapse. Therefore, developing a method for rapid, accurate, and stable control of the DC bus voltage during mode switching has crucial theoretical and engineering significance.

[0003] Currently, research on voltage stability control in DC microgrids or off-grid systems mainly focuses on several levels. At the basic control level, feedback control strategies based on proportional-integral-derivative (PID) are widely adopted to maintain basic system stability by adjusting real-time voltage deviations. To address the current sharing problem of multiple converters in parallel operation, droop control and its improved methods are commonly used to achieve autonomous power distribution by simulating the droop characteristics of synchronous generators. At a more advanced control level, researchers have introduced modern control theories such as model predictive control (MPC) and sliding mode variable structure control. These methods typically rely on establishing a relatively accurate mathematical model of the system and improving the dynamic response speed and robustness of the system through online optimization or preset switching logic. In addition, with the development of artificial intelligence technology, some studies have begun to explore using fuzzy logic or basic neural networks to tune the parameters of traditional PID controllers online, aiming to improve the adaptive capability of the control strategy to system changes. These methods have improved the control performance of DC bus voltage to a certain extent, laying the foundation for stable system operation.

[0004] However, facing the inherent strong nonlinearity, rapid time-varying nature, and high uncertainty of off-grid photovoltaic hydrogen production systems, existing technologies still have some limitations in terms of in-depth dynamic feature extraction, accurate system state characterization, and forward-looking control strategies. For example, existing methods, when attempting to extract dynamic features from voltage signals, use signal processing tools such as the Hilbert-Huang Transform (HHT) to analyze the DC bus voltage signal and extract the so-called instantaneous frequency as a feature. However, since the Hilbert Transform excels at analyzing narrowband signals, and the DC bus voltage is essentially a DC bias signal superimposed with ripple and disturbances, it is not entirely suitable. Therefore, the extracted features cannot truly and reliably reflect the dynamic behavior of the system, making it difficult to provide effective guidance for control, and may even be misleading. Summary of the Invention

[0005] The purpose of this invention is to provide a method for stabilizing the DC bus voltage during the switching of the operating mode of a photovoltaic off-grid hydrogen production system, so as to solve the above-mentioned problems existing in the prior art.

[0006] Technical solution: A method for stabilizing DC bus voltage during the instantaneous switching of operation modes in a photovoltaic off-grid hydrogen production system, including:

[0007] The preprocessed DC bus voltage signal is acquired, and the high-dimensional phase space trajectory of the voltage signal is reconstructed using an adaptive delay time.

[0008] Calculate the phase space geometric characteristics of high-dimensional phase space trajectories, including phase space curvature and phase space twist.

[0009] In the historical trajectory database, based on the current high-dimensional phase space trajectory and phase space geometric features, historical similar trajectory segments are searched, and based on the known evolution paths of historical similar trajectory segments, the final predicted trajectory for future moments is generated.

[0010] Based on the deviation between the final predicted trajectory and the reference voltage, a feedforward control signal is generated, and combined with a preset feedback control signal to form the final control command.

[0011] Beneficial effects: This invention achieves accurate causal prediction of voltage trajectories, solving the problem of control performance degradation caused by model mismatch; at the same time, it also realizes geometrically adaptive feedforward control, effectively suppressing DC bus voltage overshoot and oscillation, shortening settling time, and ensuring the speed, smoothness and robustness of control. Attached Figure Description

[0012] Figure 1 A flowchart illustrating the steps of a method for stabilizing DC bus voltage during the switching of operating modes in a photovoltaic off-grid hydrogen production system, as provided in this application embodiment.

[0013] Figure 2A flowchart illustrating the steps for reconstructing the high-dimensional phase space trajectory of a voltage signal according to an embodiment of this application.

[0014] Figure 3 A flowchart illustrating the steps for calculating the phase space geometric features of a high-dimensional phase space trajectory, as provided in an embodiment of this application.

[0015] Figure 4 A flowchart illustrating the steps for generating the final predicted trajectory at future moments, as provided in an embodiment of this application. Detailed Implementation

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

[0017] It should be noted that the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion, for example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0018] The study found that most existing control methods have a single dimension of perception for voltage dynamics, resulting in a lack of dynamic morphology awareness. Traditional methods treat voltage signals as one-dimensional scalar time series, and their control logic primarily responds to voltage amplitude deviations. This ignores the rich geometric information contained in the voltage evolution trajectory within phase space. When a system transitions from one stable point to another, the shape of its voltage trajectory (such as the severity of bending and the complexity of distortion) profoundly reflects the dynamic characteristics of the system. Existing methods cannot effectively distinguish between two fundamentally different dynamic processes: one is a benign transient process with a smooth trajectory, and the other is a malignant precursor to instability with a sharply bending or distorted trajectory. This lack of awareness of dynamic geometry leads to a lack of sufficient foresight in the controller, resulting in a passive response to deviations and hindering proactive, preventative, and precise regulation.

[0019] Furthermore, control methods relying on precise models face significant challenges in robustness and adaptability. Off-grid photovoltaic hydrogen production systems are complex, time-varying systems. Photovoltaic output fluctuates dramatically due to weather conditions, and the equivalent impedance of the electrolyzer varies with operating current, temperature, and aging. This makes establishing a globally accurate mathematical model that remains accurate under all operating conditions exceptionally difficult. Once the actual system parameters deviate from the model's preset values, the performance of methods such as model predictive control deteriorates significantly, especially during drastic dynamic processes like mode switching, where model mismatch problems are amplified dramatically, leading to poor control performance or even failure. Therefore, existing technologies generally lack a predictive and control capability that does not rely on precise mechanistic models and can autonomously learn from data and adapt to changes in system state.

[0020] like Figure 1 As shown, a method for stabilizing DC bus voltage during the switching of operation modes in a photovoltaic off-grid hydrogen production system is proposed, including the following steps:

[0021] The preprocessed DC bus voltage signal is acquired, and the high-dimensional phase space trajectory of the voltage signal is reconstructed using an adaptive delay time.

[0022] In this embodiment, the raw DC bus voltage signal of the off-grid photovoltaic hydrogen production system is acquired by a voltage sensor at a sampling frequency of, for example, 50 kHz. Since the acquired signal typically contains high-frequency noise and measurement interference, preprocessing is required. Preprocessing can employ filtering methods well-known to those skilled in the art, such as adaptive Kalman filters or Savitzky-Golay filters, to obtain a smoother preprocessed DC bus voltage signal that better reflects the true dynamics of the system, denoted as V. clean (t). In order to reveal the implications of the one-dimensional time series V clean The high-dimensional dynamic characteristics of the system underlying (t) are used to reconstruct its phase space trajectory using a time-delay embedding method. Specifically, the embedding dimension m and the adaptive delay time τ are selected. adaptive An m-dimensional state vector X(t) is constructed to represent the system state at time t. This state vector constitutes the high-dimensional phase space trajectory of the voltage signal. Its calculation formula is: X(t) = [V clean (t), V clean (t -τ adaptive V clean (t - 2τ adaptive ), ..., V clean (t - (m-1)τ adaptive [); where X(t) is the m-dimensional state vector at time t; V clean (t) represents the preprocessed voltage value at time t; τ adaptiveThe adaptive delay time, which can be adjusted according to changes in the system's dynamic characteristics, is crucial for ensuring the quality of phase space reconstruction. m represents the embedding dimension, typically chosen based on the system's complexity; for example, m=3 or m=4. This embodiment uses a one-dimensional time-series signal V... clean The t is expanded into a trajectory X(t) in a high-dimensional space, which can more comprehensively characterize the dynamic evolution process of the system.

[0023] Calculate the phase space geometric characteristics of high-dimensional phase space trajectories, including phase space curvature and phase space twist.

[0024] In this embodiment, after reconstructing the high-dimensional phase space trajectory X(t), its local geometric features are calculated to quantify its dynamic morphology. These geometric features can intuitively reflect the drasticness and complexity of voltage changes. Specifically, two core geometric features are mainly calculated: phase space curvature κ. phase (t): This feature quantifies the curvature of the trajectory at time t. A larger curvature indicates that the voltage trajectory is undergoing a sharp turn, corresponding to a rapid change in system state, such as a sudden load connection or disconnection. Phase space torsion τ phase (t): This feature quantifies the degree to which the trajectory deviates from its oscillating plane at time t, reflecting the three-dimensional twisting characteristics of the trajectory. The magnitude of the twist is related to the complexity and nonlinearity of the system dynamics. Optionally, the two geometric features can be calculated using Frenet-Serret frame theory based on a correct phase-space trajectory constructed from time-delay embeddings.

[0025] Based on the current high-dimensional phase space trajectory and phase space geometric features in the historical trajectory database, similar historical trajectory segments are searched, and based on the known evolution paths of the similar historical trajectory segments, the final predicted trajectory for future moments is generated.

[0026] In this embodiment, the system's past behavior is used to predict its future evolution. Specifically, a historical trajectory database needs to be constructed and maintained, which stores the system's phase space trajectory X(t) and its corresponding geometric features [κ] over a past period of time (e.g., the past hour). phase (t), τ phase When future predictions are needed, the following steps are performed: The current state X(t) and its geometric features are compared with all historical state points in the historical trajectory database. Instead of a simple Euclidean distance, a geometrically enhanced geodesic distance metric d is used for this comparison. enhancedThis metric considers both the similarity of the state vectors themselves and the similarity of their dynamic geometric shapes. Using this metric, the k most similar historical trajectory segments to the current state are retrieved from the historical database. For each retrieved historical similar trajectory segment, its subsequent evolutionary path is known. Utilizing this known evolutionary information, a weighted average is used to causally generate a predicted trajectory X for the future time t+Δt. pred (t+Δt), where Δt is the time increment. This causally emphasizes that the prediction is based entirely on historical information up to and including time t, without utilizing any future data. The final X is... pred (t+Δt) is a nonlinear combination of multiple historical possibilities, which constitutes a robust estimate of the future state.

[0027] Based on the deviation between the final predicted trajectory and the reference voltage, a feedforward control signal is generated, and combined with a preset feedback control signal to form the final control command.

[0028] Specifically, in obtaining the predicted trajectory X at future moments... pred After (t+Δt), it is used to generate control commands. From the predicted trajectory X pred Extract the core's predicted voltage component V from (t+Δt). pred (t+Δt). Calculate the predicted voltage and its relationship to the system's set reference voltage V. ref Prediction deviation e between (e.g., 800V) pred = V ref - V pred (t+Δt). Based on this prediction bias e pred And the geometric features of the predicted trajectory (such as the predicted curvature κ). pred ), generating a forward-looking feedforward control signal U feedforward This feedforward signal can anticipate and compensate for impending voltage fluctuations. Meanwhile, the traditional feedback control signal U... feedback It is also calculated in parallel, based on real-time voltage measurements V. clean (t) and reference voltage V ref real-time deviation e current This is typically generated by a proportional-integral (PI) controller. The feedforward control signal U... feedforward With feedback control signal U feedback To achieve collaborative integration and form the final control command U total Fusion typically employs a soft-switching mechanism to ensure a smooth, disturbance-free switching of the system when feedforward control is introduced. The final control command U... total It is converted into PWM duty cycle and other forms and sent to the DC-DC converter or related execution unit in the system to achieve stable control of the DC bus voltage.

[0029] This embodiment can be applied to photovoltaic off-grid hydrogen production systems to address the problem of severe fluctuations in DC bus voltage caused by sudden changes in illumination or start-up and shutdown of electrolyzers during operation mode switching. Through predictive control strategies, it achieves rapid and stable voltage control.

[0030] like Figure 2 As shown, according to one aspect of this application, reconstructing the high-dimensional phase space trajectory of a voltage signal includes:

[0031] Using an adaptive delay time as the basic interval, the voltage value at the current moment and the voltage values ​​at at least two different historical moments are sequentially extracted from the preprocessed DC bus voltage signal. The voltage value at the current moment and the voltage values ​​at at least two historical moments are combined in chronological order to form a multidimensional state vector. This multidimensional state vector constitutes a high-dimensional phase space trajectory that can reveal the inherent dynamic characteristics of the system.

[0032] In the phase space reconstruction of this embodiment, the choice of delay time τ is crucial: if τ is too small, the one-dimensional voltage signals V(t) and V(t-τ) are highly correlated, failing to provide new information, causing the reconstructed trajectory to cluster near the diagonal and unable to expand effectively; if τ is too large, the one-dimensional voltage signals V(t) and V(t-τ) may be completely uncorrelated, disrupting the dynamic continuity of the system. Specifically, the method for determining the adaptive delay time includes: performing statistical or spectral analysis on the time series characteristics of the preprocessed DC bus voltage signal to extract physical feature time scales that characterize the inherent time rhythm of the system's dynamic response; and calculating and obtaining the adaptive delay time for high-dimensional phase space trajectory reconstruction based on the extracted physical feature time scales through a preset ratio or functional relationship.

[0033] Optionally, physical feature timescales capable of characterizing the intrinsic time rhythm of the system's dynamic response are extracted, including:

[0034] The normalized autocorrelation function (NAC) of the preprocessed DC bus voltage signal is calculated, which describes the degree of self-similarity of the signal under different time delays. From the NAC, the minimum positive delay time at which the function value first crosses zero is identified as the characteristic decorrelation time of the system. From the NAC, the delay time when the function value decays from 1 to 1 / e is determined as the characteristic decay time of the system. By performing a geometric average of the characteristic decorrelation time and the characteristic decay time, a physical characteristic time scale for robustly evaluating the system's memory is obtained.

[0035] In this embodiment, the intrinsic time scale is determined by analyzing the self-similarity of the signal. Specifically, the preprocessed voltage signal V is calculated. cleanThe normalized autocorrelation function ACF(τ) of (t) is given by: ACF(τ) = E[(V clean (t) - μ v ) * (V clean (t+τ) -μ v )] / σ v 2 Where ACF(τ) is the autocorrelation function value with a delay time of τ; E[·] represents the mathematical expectation operation; μ v For signal V clean The mean of (t); σ v 2 For signal V clean The variance of ACF(t) is given by τ, where τ is the time delay. This function describes the self-similarity of the signal at different time delays, with ACF(0) = 1. ACF(τ) typically decays as τ increases. Two key time scales are identified from the calculated ACF(τ) function: the feature decorrelation time T. decorr The first positive delay is the time it takes for the ACF(τ) function value to first decrease from 1 and cross zero. This time point marks the point at which the signal becomes linearly independent of itself in time, and is an important time marker for the system to forget its past states. The characteristic decay time T1 / e is the delay time corresponding to the ACF(τ) function value decaying from 1 to 1 / e (approximately 0.368). This time constant is often used to describe the rate of exponential decay. To obtain a more robust assessment of the system's memory, the final physical characteristic timescale T can be obtained by geometrically averaging these two times. char :T char = sqrt(T decorr * T1 / e); where sqrt(·) is the square root operation. Based on this physical characteristic time scale, the adaptive delay time τ for phase space reconstruction is calculated through a preset proportional relationship. adaptive A preferred setting is: τ adaptive = T char / m; where m is the embedding dimension. For example, if the embedding dimension m=3, then τ adaptive = T char / 3. This ensures that the components of the embedding vector are both sufficiently decorrelated and retain the necessary dynamic relationships.

[0036] Optionally, the physical feature timescale that can characterize the intrinsic time rhythm of the system's dynamic response can also be:

[0037] Apply a fast Fourier transform to the preprocessed DC bus voltage signal to calculate the power spectral density of the DC bus voltage signal; identify the characteristic frequency in the power spectral density by peak detection or centroid calculation; and determine the reciprocal of the characteristic frequency as the physical characteristic time scale that can reflect the oscillation period of the system.

[0038] In this embodiment, the main oscillation period of the signal is determined by analyzing its energy distribution in the frequency domain. Specifically, the preprocessed voltage signal V... clean A Fast Fourier Transform (FFT) is applied to a sliding time window (e.g., containing 2048 sampling points) of signal (t) to calculate its power spectral density (PSD), denoted as PSD(f). PSD(f) reveals the distribution of signal energy at different frequencies f. From the calculated PSD(f), characteristic frequencies f that represent the system's main oscillation modes or average frequency behavior are identified using peak detection algorithms (e.g., finding local maxima exceeding three times the average power spectral density) or spectral centroid calculation. char The reciprocal of this characteristic frequency is determined as the physical characteristic time scale T that reflects the main oscillation period of the system. char :T char = 1 / f char Correspondingly, the adaptive delay time τ adaptive It can be set to a fraction of the period, for example: τ adaptive = T char / 4; For a periodic signal, when the sampling interval is 1 / 4 of its period, its sine and cosine components can be captured well, thereby reconstructing its dynamic characteristics.

[0039] Optionally, the physical feature timescale that can characterize the intrinsic time rhythm of the system's dynamic response can also be:

[0040] In a high-dimensional phase space trajectory, a pair of trajectory points that are adjacent to each other in the initial state are selected, and the distance between the evolution paths of the adjacent pair of trajectory points changes with time. By performing an exponential fitting on the distance change with time, the maximum Lyapunov exponent, which characterizes the average exponential divergence rate of the high-dimensional phase space trajectory, is estimated. The reciprocal of the absolute value of the maximum Lyapunov exponent is defined as the Lyapunov time of the system, and this Lyapunov time is used as the physical characteristic time scale for measuring the predictability time range of the system.

[0041] In this embodiment, from the perspective of nonlinear dynamics, the predictability time range of the system trajectory is determined by measuring the divergence rate of the trajectory. Specifically, in the high-dimensional phase space trajectory constructed by time delay embedding (where an initial τ value can be used, such as several times the sampling period), an initial state point X(t) is selected, and a neighboring reference point X is selected within its smallest neighborhood (e.g., at a distance less than 1% of the signal standard deviation). ref (t). Tracing the evolution paths of these two trajectories and calculating the distance d(t) between them. i ) = ||X(t i ) - X ref (t i )||With time t i The change in distance d(t). i The function d(t) is used to perform an exponential fit on the change over time. i ) ≈ d(t0) * exp(λ max * (t i - t0), where t0 is the initial time, can be used to estimate the maximum Lyapunov exponent λ, which characterizes the average exponential divergence rate of the phase space trajectory. max The reciprocal of the absolute value of the maximum Lyapunov exponent is defined as the Lyapunov time T of the system. lyap And use this as the physical characteristic time scale T char :T char = T lyap = 1 / |λ max The Lyapunov time represents the upper bound on the predictability of the chaotic behavior of a system. It is chosen as the characteristic time scale, and τ is set accordingly. adaptive = T char / m allows the reconstructed phase space to still have predictive significance when the system enters a state of chaos or complex dynamics.

[0042] This embodiment ensures that the selection of the adaptive delay time—a key parameter for phase space reconstruction—is always optimal and adaptable. Three methods are provided to calculate the time scale of physical characteristics from different dimensions and based on physical or statistical properties: the autocorrelation function method from the perspective of time memory, the spectral analysis method from the perspective of oscillation periodicity, and the Lyapunov time law from the perspective of the fundamental dynamic characteristics of the system's predictability boundary. In off-grid photovoltaic hydrogen production systems, the dynamic characteristics of the voltage signal change when the light intensity changes or the electrolyzer operating conditions change. The adaptive time scale calculation method can capture these changes in real time and dynamically adjust the adaptive delay time, ensuring that the reconstructed phase space trajectory always most clearly reflects the current dynamic structure of the system, improving the robustness and reliability of the entire analysis framework over a wide range of operating conditions.

[0043] like Figure 3 As shown, according to one aspect of this application, calculating the phase space geometric characteristics of a high-dimensional phase space trajectory includes:

[0044] The unit tangent vector of the trajectory is determined based on the displacement between adjacent data points on the high-dimensional phase space trajectory.

[0045] In this embodiment, the high-dimensional phase space trajectory is composed of a series of discrete state vectors {..., X(t-Δt), X(t), X(t+Δt),...} arranged in time order, where Δt is the time interval for data sampling. To determine the unit tangent vector T(t) at time t, the displacement vector ΔX(t) from time t to the next time t+Δt is calculated: ΔX(t) = X(t+Δt) - X(t); where X(t) is the m-dimensional state vector at time t; and X(t+Δt) is the m-dimensional state vector at time t+Δt. The direction of this displacement vector ΔX(t) approximates the tangent direction of the trajectory at time t. To eliminate the influence of its length, it is normalized to obtain the unit tangent vector T(t): T(t) = ΔX(t) / ||ΔX(t)||; where ||·|| represents the Euclidean norm (i.e., the magnitude) of the vector. T(t) is a vector with a magnitude of 1, which precisely points to the direction of the trajectory at time t.

[0046] Based on the rate of change of the unit tangent vector along the arc length of the trajectory, the principal normal vector is calculated, and the magnitude of this rate of change is determined as the phase space curvature.

[0047] In this embodiment, the phase space curvature κ phase The curve (t) measures the curvature of the trajectory. In the discrete case, directly calculating the rate of change of the unit tangent vector along the arc length is highly sensitive to noise. Therefore, a more robust geometric fitting method is preferred. Specifically, a three-point circle fitting method is used. For the current point X(t) on the trajectory and its two adjacent points X(t-Δt) and X(t+Δt), these three points uniquely define a plane and a circle in space. The radius R of this circle is then calculated. circle (t), thus obtaining the instantaneous curvature at time t: κ phase (t) = 1 / R circle (t); where κ phase (t) represents the phase space curvature at time t; R circle (t) is the radius of the fitted circle passing through the three points X(t-Δt), X(t), and X(t+Δt). The greater the curvature of the trajectory, the more easily it can be approximated by a circle with a smaller radius locally. To further improve the robustness of the computation, in some alternative implementations, the following optimization can be performed: when the three points are approximately collinear (i.e., R... circle(t) is maximized), and the results of three-point circle fitting are unstable. In this case, it is possible to automatically switch to five-point quadratic curve fitting or other higher-order fitting methods to calculate the curvature at that point, thereby improving accuracy. For the calculated original curvature sequence κ... phase The curve (t) can be smoothed using a recursive median filter or an exponentially weighted moving average filter to remove artifact pulses caused by signal noise. Simultaneously, the principal normal vector N(t) is determined. It points to the concave side of the trajectory, i.e., towards the center of the fitted circle, and is orthogonal to the unit tangent vector T(t).

[0048] The binormal vector is established by calculating the cross product of the unit tangent vector and the principal normal vector.

[0049] Specifically, the unit tangent vector T(t), the principal normal vector N(t), and the binormal vector B(t) together constitute a local right-handed orthogonal coordinate system, i.e., the Frenet-Serret frame, that accompanies the trajectory movement at time t. Once T(t) and N(t) are determined, the binormal vector B(t) can be directly calculated using the cross product: B(t) = T(t) × N(t); where × represents the cross product operation. The binormal vector B(t) is perpendicular to both T(t) and N(t), and its direction defines the direction in which the trajectory deviates from its oscling plane (the plane spanned by T(t) and N(t)).

[0050] The phase space torsion, which characterizes the degree of three-dimensional distortion of the trajectory, is quantified and obtained by projecting the rate of change of the secondary normal vector along the trajectory arc length onto the direction of the principal normal vector.

[0051] In this embodiment, the phase space twist τ phase (t) measures the tendency of the trajectory to deviate from a plane, i.e., the degree of its three-dimensional distortion. According to the Frenet-Serret formula, it is calculated as: τ phase (t) = -N(t) • (dB / ds); where • denotes the vector dot product operation; dB / ds represents the rate of change of the binormal vector B(t) along the arc length s of the trajectory. In discrete calculations, this formula is approximated as: τ phase N(t) = -N(t) • ((B(t+Δt) - B(t)) / ||ΔX(t)||); where B(t+Δt) is the binormal vector at the next moment; ||ΔX(t)|| is the magnitude of the displacement vector, approximately equal to the arc length element Δs. To handle singularity issues that may arise in practical calculations (e.g., when the trajectory is approximately a straight line and the curvature is close to zero, the definition of N(t) becomes unstable), the following measures are taken: When calculating the principal normal vector N(t), if its magnitude (i.e., ||ΔT(t)||) is less than a preset minimum threshold ε... N (e.g., 1e) -6If the value of τ is not found in the original value, then the point is considered a singular or near-singular point. In this case, the principal normal vector is not updated; instead, the value from the previous time step is retained, i.e., N(t) = N(t-Δt), to ensure the continuity of the calculation. The calculated original torsional τ is then... phase (t) is subject to a reasonable range limit based on the system's physical characteristics (e.g., [-5, 5] rad / ms) to prevent extreme values ​​caused by noise. Subsequently, it can be smoothed using filtering methods such as exponentially weighted moving averages to obtain the final, robust phase space torsion. This allows for the stable and reliable extraction of physically meaningful curvature and torsion from noisy, discrete phase space trajectories, providing crucial quantitative inputs for subsequent similarity measurements, multi-scale decomposition, and predictive control.

[0052] This embodiment reconstructs a single DC bus voltage time series to create a high-dimensional phase space trajectory that accurately reflects the system's intrinsic dynamic behavior. Based on this, geometric features such as phase space curvature and phase space torsion, calculated from this real trajectory using rigorous differential geometry tools such as the Frenet-Serret frame, are effective and have clear physical meaning. For example, the curvature directly quantifies the severity of the voltage change rate, i.e., the turning radius of the trajectory; while the torsion quantifies the tendency of the trajectory to deviate from its close plane, reflecting the complexity of the dynamic process and the degree of three-dimensional distortion. This enables the control system to, for the first time, possess the ability to gain deep geometric insight into the voltage dynamic process, providing high-quality, high-information-content input data for all subsequent algorithms (such as prediction and control), and solving the problem of dynamic morphological blindness.

[0053] According to one aspect of this application, the dynamic behavior of a photovoltaic off-grid hydrogen production system is a superposition of multiple physical processes at different time scales. For example, the charging and discharging of capacitors caused by the switching of power electronic devices is a fast-scale process on the order of milliseconds; the power regulation dynamics of the DC-DC converter are a medium-scale process on the order of tens of milliseconds; and the slow impedance change of the electrolyzer caused by temperature or operating condition changes is a slow-scale process on the order of seconds. Mixing these dynamics together and using the same model for prediction will inevitably lead to model mismatch and decreased accuracy. Therefore, before generating the final predicted trajectory for future moments, the process also includes: based on the phase space geometric characteristics, decomposing the feature space where the high-dimensional phase space trajectory is located into predetermined time-scale feature subspaces for describing different dynamic processes.

[0054] In this embodiment, the decomposition is based on phase space geometric features, as these features are closely related to the dynamic processes of the system. For example, curvature κ phase The rate of change over time dκ phase / dt can effectively characterize the intensity of dynamic changes, while the torsional τ phaseThe relative magnitude of the curvature reflects the complexity of the dynamics.

[0055] Furthermore, the decomposition is specifically as follows: data points whose temporal rate of change of phase space curvature is greater than a preset fast threshold are assigned to the fast-scale feature subspace; data points whose phase space torsion is numerically dominant and whose temporal rate of change of phase space curvature is lower than a preset fast threshold are assigned to the mesoscale feature subspace; and data points whose phase space curvature and phase space torsion are both less than a preset slow threshold are assigned to the slow-scale feature subspace.

[0056] Specifically, historical data on the geometric characteristics of the system over a period of time (e.g., the past 1000 ms) are collected, particularly the temporal rate of change of curvature |dκ. phase The distribution of / dt| is used. The probability density function of this distribution is nonparametrically fitted using the Kernel Density Estimation (KDE) method. It can be observed that this probability density function typically exhibits multiple peaks, corresponding to different dynamic modes. The Expectation-Maximization (EM) algorithm can be used to fit this mixed distribution as, for example, a superposition of three Gaussian components. These three Gaussian components correspond to fast, medium, and slow dynamic modes, respectively. The intersection point of the probability density function curves of two adjacent Gaussian components is adaptively determined as the threshold used for decomposition, for example, the fast-medium boundary threshold θ. fast_medium and the medium-slow boundary threshold θ medium_slow Based on the adaptively determined threshold, for each new data point X(t), according to its geometric features κ... phase (t) and τ phase (t), classify according to the following rules: if |dκ phase / dt| > θ fast_medium Then the point is assigned to the fast-scale feature subspace M. fast If |dκ phase / dt|≤θ fast_medium And |τ phase (t)| > |κ phase If (t)| (i.e., torsion is dominant), then this point belongs to the mesoscale eigenspace M. medium If κ phase (t) <θ medium_slow And |τ phase If (t)| is also extremely small, then the point is classified into the slow-scale eigenspace M. slow .

[0057] For each time-scale feature subspace decomposed, a dedicated trajectory prediction model is used in subsequent prediction steps.

[0058] In this embodiment, to achieve a smooth transition between different subspaces and avoid abrupt changes at the boundaries caused by rigid partitioning, a boundary fuzziness processing mechanism is introduced. Specifically, for data points located near the boundaries of two subspaces, fuzzy membership degrees can be defined. For example, for a |dκ... phase / dt| is slightly less than θ fast_medium A point can be calculated to simultaneously belong to M_ ast and M medium The membership degree (or weight). Its membership function can be a sigmoid function, such as the sigmoid function: μ fast (x) = 1 / (1 + exp(-c * (x - θ fast_medium )));where x is |dκ phase The value of / dt|; c is a coefficient controlling the steepness of the transition zone; μ fast (x) represents the membership degree of this point in the fast-scale subspace. Correspondingly, its membership degree in the mesoscale subspace is 1-μ. fast (x). Such a boundary point can be considered to belong to two subspaces simultaneously. In subsequent predictions, predictions can be made using dedicated models for each of the two subspaces, and then the prediction results can be weighted and fused based on membership. This allows for matching the most suitable prediction model to dynamic processes of different properties (e.g., using a simple geometric extrapolation model for fast-scale processes, a local linear model for mesoscale processes, and an optimal path model for slow-scale processes), thereby improving the overall prediction accuracy and robustness.

[0059] Furthermore, since we have already constructed dedicated prediction models for the three different dynamic processes (fast, medium, and slow) through multi-scale decomposition, we still need to Xize the prediction results of these three models. fast_pred X medium_pred and X slow_pred Intelligently combined. Therefore, this embodiment proposes a confidence-weighted fusion strategy based on prediction uncertainty, specifically: when combining data from M... fast M medium and M slow When predicting data points in the three subspaces, different models best suited to their dynamic characteristics are used. Because fast-scale processes are rapidly changing and highly nonlinear, but lack long-term memory, the fast-scale subspace M... fast The prediction can employ causal prediction methods based on historically similar trajectory segments. Mesoscale processes (such as converter power regulation) typically exhibit certain linear dynamic characteristics. Therefore, a mesoscale subspace M can be defined. mediumThe prediction identifies local autoregressive moving average (ARMA) models or state-space models. For example, by fitting historical data within the neighborhood of the current point using the least squares method, a local dynamic model X(t+1) = A·X(t) + B·U(t) is identified online, and multi-step predictions are performed based on this model; where X(t) is the current state vector, A is the state transition matrix, B is the control input matrix, and U(t) is the control input vector. Slow-scale processes (such as changes in electrolytic cell impedance) change gradually and have smooth trajectories. Therefore, a slow-scale subspace M can be defined. slow The prediction constructs a local quadratic surface model f(x) = x T A*x + b T The predicted trajectory is calculated as x + c, where x is the current state or position vector, A* is a symmetric matrix, b is the linear coefficient vector, and c is a constant term. T This is a transpose.

[0060] For each model's prediction, its reliability needs to be evaluated. This can be achieved by quantifying the uncertainty of the prediction. Specifically, during the prediction process in each subspace, a measure of the dispersion of the prediction results can be obtained. For example, in fast-scale prediction, the variance or covariance matrix of multiple historical evolution increments used for weighted averaging can be calculated, and its trace or largest eigenvalue can be used as the prediction uncertainty σ. 2 fast For mesoscale and slow-scale models, the prediction uncertainty σ can be estimated using the residuals of the model fit. 2 medium and σ 2 slow These uncertainty measures are converted into confidence scores C ranging from [0, 1]. i A preferred calculation method is to use the exponential decay function: C i = exp(-σ 2 i / σ 2 ref ); where C i Let σ be the confidence level for the i-th scale (where i can be fast, medium, or slow); 2 i The prediction uncertainty at this scale; σ 2 refThe reference variance is a hyperparameter determined based on the overall noise level of the system, used to adjust the confidence level's sensitivity to uncertainty. Higher uncertainty results in lower confidence levels. Besides the confidence level based on the model's own performance, the importance of predictions at different scales is also related to the system's macroscopic operational stage. Identifying the system's current operational stage can be based on the relative time t after the mode switching trigger signal. relative To determine this, different time weight vectors W are set for different stages. time = [W time_fast W time_medium W time_slow ]. Initial startup phase (e.g., 0 ≤ t) relative < 5ms): During this phase, the system response is drastic, requiring priority to be given to tracking fast-moving events. Therefore, W can be set... time = [2.0, 1.0, 0.5], giving the highest time weight to fast-scale predictions. Transition period (e.g., 5ms ≤ t) relative < 15ms): During this stage, the system is transitioning from one steady state to another, with dynamic interleaving across various scales. A smooth transition weighting function, such as a cosine function, can be used to smoothly transition the weights from the initial settings to the steady-state settings. The steady-state period (e.g., t...) relative ≥ 15ms): During this stage, the system gradually stabilizes, requiring greater attention to steady-state accuracy and suppression of slow drift. Therefore, W can be set... time = [0.5, 1.0, 2.0], giving the highest time weight to slow-scale predictions.

[0061] Confidence level C i With time weight W time_i Multiplying these results yields the combined influence at each scale. Normalization is then performed to calculate the final fusion weight α. i : α i = (C i * W time_i ) / Σ(C j * W time_j The prediction trajectory X is obtained by averaging the prediction results from the three scales: fast, medium, and slow. The fusion weights are then used to calculate the final predicted trajectory X. final_pred :X final_pred = α fast * X fast_pred + α medium * X medium_pred + α slow * X slow_pred ;where α fast α medium αslow These are the fusion weights for fast-scale, mesoscale, and slow-scale predictions, respectively; for the fused X... final_pred Perform a rationality check. For example, check whether its voltage component exceeds the system's physical operating range (e.g., [0.8*V]). nom 1.2*V nom ], V nom The system's rated voltage (or the jump between the current point and the previous trajectory point X(t)) is considered. If an anomaly is detected, soft limiting or exponential smoothing can be used for correction to ensure that the final predicted trajectory is physically reasonable and continuous. This embodiment can dynamically and intelligently combine information from different dedicated models, making full use of the advantages of each model, thereby obtaining more accurate and robust prediction results than any single model during complex operation mode switching processes.

[0062] like Figure 4 As shown, according to one aspect of this application, generating the final predicted trajectory for future moments includes:

[0063] In the historical trajectory database, based on a distance metric that integrates high-dimensional phase space trajectories and phase space geometric features, a predetermined number of historical trajectory segments most similar to the current state are retrieved.

[0064] Specifically, using the geometrically enhanced distance metric d enhanced The current state vector X(t) is compared with each historical state point X in the historical trajectory database. hist_i (t i (where t) i The comparison is performed using historical moments (strictly smaller than the current moment t to ensure causality). d is calculated. enhanced (X(t), X) hist_i (t i The values ​​of k are calculated and sorted in ascending order of distance. The k with the smallest distance is selected. hist A historical state point, for example, k hist =10. These points are considered to be the historical states most similar to the current state, and the trajectory segments they belong to are the retrieved historical similar trajectory segments.

[0065] For each retrieved historical trajectory segment, a corresponding similarity weight is assigned based on the degree of similarity between the historical trajectory segment and the current state.

[0066] In this embodiment, for each retrieved historical similarity point X hist_i (t i ), a similarity weight w needs to be assigned. i This weight should be positively correlated with the degree of similarity, i.e., with d. enhancedDistance is negatively correlated. A preferred method for weighting is to use a Gaussian kernel function: w i = exp( -d enhanced (X(t), X) hist_i (t i )) 2 / (2 *σ d 2 ) ); where w i σ is the weight of the i-th historical similarity point; exp(·) is the exponential function; σ d The distance-bandwidth parameter controls the rate at which the weights decay with distance; it can be set to k. hist The median of the distance values. Using this formula, historical points closer to the current state will be assigned a higher weight.

[0067] Extract the known evolution increments of each historical trajectory segment from the starting point to the predetermined time interval in the future. Perform a weighted average of the known evolution increments based on similarity weights and fuse them into the expected evolution vector. Superimpose the expected evolution vector onto the high-dimensional phase space trajectory at the current moment to construct the final predicted trajectory.

[0068] In this embodiment, for each selected historical similarity point X hist_i (t i Since its subsequent evolution is known, its true evolutionary increment ΔX within a future prediction step Δt can be extracted. hist_i ΔX hist_i =X hist_i (t i +Δt) - X hist_i (t i This represents the actual direction and magnitude of the system's evolution at a point in the past similar to the current state. This is achieved by analyzing all historical evolution increments ΔX. hist_i The similarity weight w is used to perform the calculation. i The weighted average of the determined values ​​is fused into the expected evolution vector ΔX. expected ΔX expected = Σ(w i * ΔX hist_i ) / Σ(w i Σ(·) represents the summation operation. The expected evolution vector can be understood as a comprehensive estimate of the most likely next evolution direction and magnitude of the current state based on all the most relevant historical experiences. It is not derived from a single historical precedent, but is the result of calculations from multiple historical possibilities, thus possessing better robustness. Superimposing the expected evolution vector onto the high-dimensional phase space trajectory X(t) at the current moment constructs the final predicted trajectory X at the future time t+Δt. pred (t+Δt): X pred(t+Δt) = X(t) + ΔX expected It can be seen that the superposition here is not a simple linear addition, but a non-linear, data-driven prediction process. It predicts the future by adding the evolutionary direction intelligently inferred from historical data to the current state. This ensures the causality of the prediction while improving its accuracy and resistance to noise by integrating multiple historical experiences. Furthermore, it can also be based on k... hist The degree of dispersion of the weighted evolution increments is used to calculate the uncertainty or confidence level of the prediction result, such as calculating its covariance matrix, to provide more dimensional information for subsequent control decisions.

[0069] According to one aspect of this application, a distance metric that integrates high-dimensional phase space trajectories and phase space geometric features is calculated as follows:

[0070] All state points in the historical trajectory database and the state point at the current moment are considered as a high-dimensional point set, and a nearest neighbor graph is constructed on the high-dimensional point set; each state point is connected only to its preset number of nearest neighbor state points.

[0071] In this embodiment, the system's phase space trajectory is not uniformly distributed throughout the high-dimensional space, but is constrained to one or more low-dimensional nonlinear manifolds. To approximate the intrinsic geometry of this manifold, a network graph reflecting its local topology needs to be constructed. Specifically, all state vectors {X} in the historical trajectory database are... hist_i The points are merged with the current state vector X(t) to form a large-scale high-dimensional point set. For each point in this set, the Euclidean distance between it and all other points is calculated, and the k nearest points are identified as its nearest neighbors. For example, k can be a fixed integer, such as k=15, or dynamically determined according to the total number of data points N, such as k=floor(log2(N)+6). An edge is established between each point and its k nearest neighbors, with the weight of the edge set to the Euclidean distance between them. This constructs a weighted undirected graph, namely the k-Nearest Neighbor Graph (k-NN Graph), which can be viewed as a discrete approximation of the system's dynamical manifold.

[0072] The distance between any two state points is determined as the shortest path length on the nearest neighbor graph connecting these two state points, i.e., the distance metric. The shortest path length aims to approximate the intrinsic distance between the two points on the system's dynamical manifold.

[0073] Furthermore, the determination of the shortest path length includes: defining the shortest path length as an approximation of the geodesic distance between two state points; and calculating the approximation of the geodesic distance iteratively by executing a shortest path search algorithm on the nearest neighbor graph.

[0074] In this embodiment, geodesic distance refers to the shortest path length between two points on a surface or manifold. Directly calculating the geodesic distance on a high-dimensional nonlinear manifold is difficult. Therefore, a method of calculating the shortest path length on a nearest neighbor graph is used to effectively approximate it. Specifically, two state points to be compared (e.g., the current point X(t) and a historical point X) are used. hist_i Using and as the starting and ending points, perform Dijkstra's algorithm on the nearest neighbor graph to calculate the shortest path length, denoted as d. geodesic (X(t), X) hist_i This is used as an approximation of the geodesic distance between two points. It is understood that although Dijkstra's algorithm itself is existing technology, this invention applies it to the manifold after voltage dynamic phase space reconstruction to address the technical problem that traditional Euclidean distance cannot accurately measure dynamic similarity.

[0075] Furthermore, it also includes optimizing the distance metric, specifically: extracting the phase space geometric features of the two state points to be compared; calculating the geometric adjustment factor based on the difference in phase space curvature and phase space torsion between the two state points; and obtaining the final distance metric by multiplying the approximate value of the geodesic distance with the geometric adjustment factor.

[0076] In this embodiment, using only geodesic distance is insufficient because it only considers the point's position on the manifold, without taking into account the dynamic shape of that position. Therefore, a geometric enhancement step is introduced to incorporate phase space geometric features into the distance metric. Specifically, for two state points X to be compared... i and X j Extract their corresponding phase space curvature κ respectively phase_i κ phase_j Phase space twist τ phase_i τ phase_j Based on these differences in geometric features, the geometric adjustment factor G is calculated. factor (i, j). The greater the difference in local geometry between two points, the larger the value of this factor, thus penalizing the distance. Its calculation formula can be: G factor (i, j) = 1 + α κ * |κ phase_i -κ phase_j | / κ max +α τ * |τphase_i -τ phase_j | / τ max ; where G factor (i, j) is the geometric adjustment factor; α κ and α τ The preset geometric weighting coefficients, such as α κ =0.3、α τ =0.2, used to balance the contributions of curvature and torsion; κ max and τ max The maximum value of each feature in the historical database is used for normalization. The geodesic distance approximation is multiplied by this geometric adjustment factor to obtain the final geometrically enhanced distance metric d. enhanced (i, j): d enhanced (i, j) = d geodesic (i, j) * G factor (i, j). The final distance metric d enhanced It measures not only the proximity of two points on the dynamic manifold, but also the similarity of their local dynamic behaviors. Only when two state points are close to each other on the manifold, and their voltage trajectories bend and twist in similar ways, is the d-value between them considered to be similar. enhanced Only then will the distance truly be small. This improves the accuracy of similarity retrieval.

[0077] This embodiment addresses the problem that Euclidean distance cannot measure true distances on nonlinear manifolds by constructing a nearest neighbor graph and calculating the shortest path length on the graph. Furthermore, it creates a composite metric by fusing geodesic distance with the local geometric features (curvature, torsion) of state points. This enables the system to find historical recurrences that not only have similar state values ​​but also highly consistent dynamic evolution behavior when searching for historical similarities, thereby improving prediction accuracy.

[0078] According to one aspect of this application, generating a feedforward control signal includes:

[0079] The final predicted trajectory is analyzed to obtain the predicted voltage at future time points, and the predicted voltage deviation between the predicted voltage and the reference voltage is calculated.

[0080] Differential geometric analysis is performed on the final predicted trajectory to obtain the predicted phase space curvature or phase space torsion.

[0081] Based on the predicted voltage deviation, the basic feedforward control quantity is calculated. Using the predicted phase space curvature or phase space torsion, one or more control gains of the basic feedforward control quantity are adjusted to generate a feedforward control signal that anticipates and adapts to the dynamic shape of the trajectory.

[0082] In this embodiment, the feedforward control signal Ufeedforward It consists of two control components with clearly defined objectives, which are fused together through a coordination mechanism. Specifically, these components include: a voltage compensation component, which is mainly calculated based on the magnitude of the predicted voltage deviation, aiming to quickly align the future voltage with the reference voltage; and a trajectory adjustment component, which is mainly calculated based on the difference between the predicted phase space curvature and the target curvature, aiming to actively shape the dynamic response process of the voltage to ensure a smooth transition. By weighting and superimposing or coordinating the voltage compensation component and the trajectory adjustment component, they together constitute the feedforward control signal.

[0083] Specifically, from the final predicted trajectory X pred Extract the predicted voltage component V from (t+Δt). pred (t+Δt). Calculate its relationship with the reference voltage V. ref The predicted voltage deviation e between (e.g., 800V) pred = V ref - V pred (t+Δt). Based on this prediction deviation, the voltage compensation component U is calculated. voltage_ff To provide stronger control under large deviations, a nonlinear gain design can be employed: U voltage_ff = K v * e pred * (1 + γ e * |e pred | / V nom ); where U voltage_ff For voltage compensation component; K v Based on the base voltage feedforward gain, for example, K v =1.2; γ e For nonlinear compensation coefficients, such as γ e =0.3;|e pred | represents the absolute value of the prediction bias; V nom The system's rated voltage is used for normalization. For the final predicted trajectory X... pred Differential geometric analysis is performed on (t+Δt) to obtain its predicted phase space curvature κ. pred Calculate the predicted curvature and the preset target curvature κ. target (e.g., κ) target =0.2, representing a smooth transition, the difference Δκ=κ between target -κ pred Based on this curvature difference, the trajectory adjustment component U is calculated. shape_ff :U shape_ff = K κ *Δκ; where U shape_ff Adjust the trajectory components; K κ For curvature control gain, such as K κ=0.8. When the predicted trajectory is more curved than the target (κ). pred > κ target When the voltage compensation component is at its maximum (0), it will have a suppressive effect; conversely, it will have an accelerating effect. The fundamental gain K of the voltage compensation component can be... v The design is such that K is a function that varies with the predicted geometric features, thus achieving adaptive adjustment. For example, K can be set to... v Compared with the predicted phase space curvature κ pred Related: K vadaptive = K v0 * (1 +γ κ *κ pred / κ nominal ); where K vadaptive K represents the adaptively adjusted voltage gain. v0 Based on the gain; γ κ κ is the curvature sensitivity coefficient. nominal This is the nominal curvature value. When the system predicts a sharp turn in the future trajectory, it automatically increases the feedforward gain to provide stronger control. The voltage compensation component, after adaptive gain adjustment, is weighted and superimposed or fused with the trajectory adjustment component to form the final feedforward control signal U. feedforward :U feedforward = w v * U voltage_ff + w s * U shape_ff ;where w v and w s Here are the weighting coefficients, and w v + w s = 1.

[0084] This embodiment, under the severe impact of mode switching, can not only anticipate future voltage drops or overshoots and compensate in advance, but also anticipate the drastic nature of voltage changes and actively apply damping to soften the process. While ensuring that the voltage quickly recovers to the rated value, it reduces the voltage overshoot and the number of oscillations, shortens the settling time, and improves the power quality of the entire transient process.

[0085] Furthermore, the final control commands are generated, including:

[0086] Based on the deviation between the real-time measured value of the DC bus voltage signal and the reference voltage, a feedback control signal is generated through a proportional-integral controller.

[0087] When a trigger signal indicating a change in system operating mode is detected, a soft switching function is enabled to generate dynamic fusion weights that change smoothly over time.

[0088] By using dynamic fusion weights, the feedforward control signal and the feedback control signal are weighted together to obtain the final control command.

[0089] In this embodiment, a standard proportional-integral (PI) feedback controller operates continuously. It calculates the voltage V measured in real time. clean (t) and reference voltage V ref Real-time deviation e between current = V ref - V clean (t), and based on this, a feedback control signal U is generated. feedback :U feedback = K p * e current + K i *∫e current dt; where K p For proportional gain; K i The integral gain is used. When a trigger signal indicating a change in system operating mode is detected (e.g., a step change in load current), the feedforward control process is initiated. To avoid system instability caused by the sudden addition of the feedforward signal, a soft-switching function can be used to generate a smoothly changing dynamic fusion weight α. switch (t). A preferred functional form is the hyperbolic tangent function (tanh): α switch (t) = 0.5 * (1 + tanh((t - t switch ) / T switch )); where α switch (t) represents the fusion weight at time t, whose value smoothly transitions from 0 to 1; t switch To detect the moment when the handover occurs; T switch The switching time constant, for example, T switch =1ms, which controls the speed of the transition process. Using this dynamic fusion weight, the feedforward control signal U is... feedforward With feedback control signal U feedback By performing coordinated weighting, the final control command U is obtained. total :U total = α switch (t) * U feedforward + U feedback In some preferred embodiments, to avoid integral saturation when introducing feedforward, the weights of the feedback portion in the fusion formula can be adjusted accordingly, for example: U total = α switch (t) *U feedforward + (1 - α switch (t) * w fb ) * U feedback;where w fb It is a coefficient less than 1. The final control command U total It can smoothly and without disturbance introduce the rapid response capability of feedforward control in the early stage of mode switching, while always maintaining the steady-state accuracy and robustness of feedback control.

[0090] This embodiment addresses the potential conflicts and disturbances that may arise at the moment of implementation between high-performance feedforward control and basic feedback control at the engineering level. When a mode switch is detected, directly superimposing a strong feedforward control signal onto the system would create a step in the control quantity, potentially triggering a secondary impact. By introducing a soft switching function, a smoothly changing dynamic weight is generated at the moment of switching. This weight smoothly and without impact transfers control dominance from the feedback controller to the feedforward + feedback co-controller within milliseconds. This ensures the continuity of control commands and achieves seamless and disturbance-free control mode switching. It improves the stability and reliability of the system at transient event entry points, avoids the instability risks caused by control switching itself, and is a crucial guarantee for the safe and effective application of the entire advanced control algorithm, comprehensively improving the system's overall performance and user experience.

[0091] In one specific embodiment, the scenario is set as follows: system rated voltage V nom = 800V, reference voltage V ref =800V. At the current time t, the collected and preprocessed voltage value is V. clean (t) = 790V. Current adaptive delay time τ adaptive = 2ms, embedding dimension m = 3. Prediction step size Δt = 1ms. The historical trajectory database contains some historical data. For simplicity, only the 3 historical points (k) most similar to the current state are listed. hist =3) and related information. Construct the current state vector X(t) based on the time-delay embedding formula. Assume the voltage values ​​at times t-2ms and t-4ms are 785V and 782V respectively. X(t) = [V clean (t), V clean (t - τ adaptive V clean (t - 2τ adaptive [790, 785, 782]. The three historical points most similar to X(t) were retrieved from the historical database, and their information is as follows: The state vector of historical point 1 is X. hist_i (t i The coordinates are [791, 786, 783], and the geodesic distance d is... geodesic The value is 1.8, the geometric feature [κ, τ] is [0.8, 0.3], and the evolution increment ΔX at the next time step is... hist_i[+5, +6, +5]; State vector X of history point 2 hist_i (t i The distance d between the geodesic lines is [789, 784, 781]. geodesic The value is 2.5, the geometric characteristic [κ, τ] is [0.9, 0.2], and the evolution increment ΔX at the next time step is... hist_i [+6, +5, +4]; State vector X of history point 3 hist_i (t i The coordinates are [792, 788, 785], and the geodesic distance d is... geodesic The value is 3.1, the geometric feature [κ, τ] is [0.7, 0.4], and the evolution increment ΔX at the next time step is... hist_i The geometric features of the current point are [+4, +4, +3]. t =0.85, τ t =0.25]. Taking historical point 1 as an example, calculate its geometrically enhanced distance d with the current point X(t). enhanced Assume α κ =0.3, α τ =0.2, κ max =2,τ max =1. Calculate the geometric adjustment factor G. factor (t, 1): G factor (t, 1) = 1 + 0.3 * |0.85 - 0.8| / 2 + 0.2 * |0.25 - 0.3| / 1 = 1 + 0.0075 + 0.01 = 1.0175. Calculate the geometric augmentation distance d. enhanced (t, 1): d enhanced (t, 1) = d geodesic (t, 1) * G factor (t, 1) = 1.8 * 1.0175 = 1.8315. Similarly, calculate d with respect to the other two historical points. enhanced The distances are 2.54 and 3.16, respectively. Calculate the similarity weight w. i (Let σ) d (2.54 meters from the median): w1 = exp(-1.8315) 2 / (2 * 2.54 2 )) = exp(-0.26) = 0.771; w2 = exp(-2.54 2 / (2 * 2.54 2 )) = exp(-0.5) = 0.607; w3 = exp(-3.16 2 / (2 * 2.54 2)) = exp(-0.77) = 0.463. Calculate the expected evolution vector ΔX. expected :Σ(w i ) =0.771 + 0.607 + 0.463 = 1.841; Σ(w i *ΔX hist_i ) = 0.771*[5, 6, 5] + 0.607*[6, 5, 4] + 0.463*[4, 4, 3] = [3.855, 4.626, 3.855] + [3.642, 3.035, 2.428] + [1.852, 1.852, 1.389] = [9.349, 9.513, 7.672]; ΔX expected = [9.349, 9.513, 7.672] / 1.841 = [5.078, 5.167, 4.167]. Construct the final predicted trajectory X. pred (t+Δt): X pred (t+1ms) = X(t) +ΔX expected = [790, 785, 782] + [5.078, 5.167, 4.167] = [795.078, 790.167, 786.167]. This embodiment predicts that after 1ms, the system's state vector will evolve to [795.078, 790.167, 786.167]. The most critical predicted voltage value is V. pred (t+1ms) = 795.078V. Based on this predicted voltage value, the prediction deviation e can be calculated. pred = 800 - 795.078 = 4.922V, and further calculate the corresponding feedforward control signal to achieve advance compensation for the voltage drop trend.

[0092] This invention employs time-delay embedding to reconstruct the phase space trajectory and combines it with differential geometry to calculate its curvature and torsion, thereby obtaining an analysis of the system's dynamic morphology. Simultaneously, it abandons the reliance on a precise system model and achieves accurate causal prediction of the voltage trajectory through historical similarity matching based on geometrically enhanced geodesic distance metrics, solving the problem of control performance degradation caused by model mismatch. It realizes forward-looking geometrically adaptive feedforward control, where the controller not only responds to predicted voltage deviations but also actively adjusts the control strategy according to the geometric morphology of the predicted trajectory to optimize the dynamic process. Under severe operating conditions such as mode switching in off-grid photovoltaic hydrogen production systems, it effectively suppresses DC bus voltage overshoot and oscillation, shortens settling time, and ensures the speed, smoothness, and robustness of control.

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

Claims

1. A method for stabilizing the DC bus voltage at the moment of switching the operation mode of a photovoltaic off-grid hydrogen production system, characterized in that, The method comprises the following steps: acquiring a preprocessed DC bus voltage signal, and reconstructing a high-dimensional phase space trajectory of the voltage signal by using an adaptive delay time; calculating phase space geometric features of the high-dimensional phase space trajectory, including phase space curvature and phase space torsion; in a historical trajectory database, searching for a historical similar trajectory segment according to the high-dimensional phase space trajectory and the phase space geometric features at the current time, and generating a final prediction trajectory at a future time based on a known evolution path of the historical similar trajectory segment; generating a feedforward control signal based on a deviation between the final prediction trajectory and a reference voltage, and combining a preset feedback control signal to form a final control instruction.

2. The method of claim 1, wherein, The method for reconstructing the high-dimensional phase space trajectory of the voltage signal comprises the following steps: sequentially extracting, from the preprocessed DC bus voltage signal, a voltage value at the current time and voltage values at at least two different historical times based on an adaptive delay time as a basis interval; combining the voltage value at the current time and the voltage values at the at least two historical times into a multi-dimensional state vector in chronological order, and the multi-dimensional state vector constitutes the high-dimensional phase space trajectory that reveals the inherent dynamic characteristics of the system.

3. The method of claim 1, wherein, The method for calculating the phase space geometric features of the high-dimensional phase space trajectory comprises the following steps: determining a unit tangent vector of the trajectory based on the displacement between adjacent data points on the high-dimensional phase space trajectory; solving a principal normal vector according to the rate of change of the unit tangent vector along the arc length of the trajectory, and determining the modulus of the rate of change as the phase space curvature; establishing a binormal vector by calculating the vector cross product of the unit tangent vector and the principal normal vector; quantifying and obtaining the phase space torsion representing the three-dimensional twisting degree of the trajectory by projecting the rate of change of the binormal vector along the arc length of the trajectory in the direction of the principal normal vector.

4. The method of claim 1, wherein, The method for generating the final prediction trajectory at the future time comprises the following steps: in the historical trajectory database, searching for a predetermined number of historical trajectory segments most similar to the current state based on a distance measurement that comprehensively considers the high-dimensional phase space trajectory and the phase space geometric features; assigning a similarity weight to each of the searched historical trajectory segments according to the similarity degree of the historical trajectory segment to the current state; extracting a known evolution increment of each historical trajectory segment from the starting point to the future predetermined time interval, and performing a weighted average on the known evolution increments based on the similarity weights to fuse into an expected evolution vector; stacking the expected evolution vector on the high-dimensional phase space trajectory at the current time to construct the final prediction trajectory.

5. The method of claim 1, wherein, The method for generating the feedforward control signal comprises the following steps: analyzing the final prediction trajectory to obtain a predicted voltage at the future time, and calculating a predicted voltage deviation between the predicted voltage and a reference voltage; performing differential geometry analysis on the final prediction trajectory to obtain a predicted phase space curvature or phase space torsion; calculating a basic feedforward control amount based on the predicted voltage deviation, adjusting one or more control gains of the basic feedforward control amount by using the predicted phase space curvature or phase space torsion, and generating a feedforward control signal that is predictive and adaptive to the dynamic form of the trajectory.

6. The method of claim 4, wherein, The distance measurement that comprehensively considers the high-dimensional phase space trajectory and the phase space geometric features is calculated in the following manner: All state points in the historical trajectory database and the state point at the current moment are considered as a high-dimensional point set, and a nearest neighbor graph is constructed on the high-dimensional point set; each state point is connected only to its preset number of nearest neighbor state points; The distance between any two state points is determined as the shortest path length in the nearest neighbor graph connecting the two state points, i.e., the distance metric.

7. The method of claim 6, wherein, Determining the shortest path length includes: The shortest path length is defined as an approximation of the geodesic distance between two state points; An approximate value for the geodesic distance is calculated iteratively by performing a shortest path search algorithm on the nearest neighbor graph.

8. The method of claim 7, wherein, This also includes optimizations to the distance metric, specifically: Extract the phase space geometric features of the two state points to be compared; The geometric adjustment factor is calculated based on the difference in phase space curvature and phase space torsion between two state points. The final distance metric is obtained by multiplying the approximate geodesic distance with a geometric adjustment factor.

9. The method of claim 1, wherein, Before generating the final predicted trajectory for future moments, the following steps are also included: Based on the geometric features of phase space, the feature space where the high-dimensional phase space trajectory is located is decomposed into a predetermined time-scale feature subspace for describing different dynamic processes; Specifically, the decomposition is as follows: data points whose temporal rate of change of phase space curvature is greater than a preset fast threshold are assigned to the fast-scale feature subspace; data points whose phase space torsion is numerically dominant and whose temporal rate of change of phase space curvature is lower than a preset fast threshold are assigned to the mesoscale feature subspace; and data points whose phase space curvature and phase space torsion are both less than a preset slow threshold are assigned to the slow-scale feature subspace. For each time-scale feature subspace decomposed, a dedicated trajectory prediction model is used in subsequent prediction steps.

10. The method of claim 1, wherein, The final control command is generated, including: Based on the deviation between the real-time measured value of the DC bus voltage signal and the reference voltage, a feedback control signal is generated through a proportional-integral controller. When a trigger signal indicating a change in system operating mode is detected, a soft switching function is enabled to generate dynamic fusion weights that change smoothly over time. By using dynamic fusion weights, the feedforward control signal and the feedback control signal are weighted together to obtain the final control command.

Images