Modular high-voltage SVG circulating current suppression method

Abstract

The invention discloses a modular high-voltage SVG circulating current suppression method. The method comprises the steps of obtaining and preprocessing module circulating current, voltage and state data; constructing an enhanced phase space fused with multiple physical quantities, extracting circulating current phase features and generating dynamic trend prediction data; calculating optimal feed-forward and feedback virtual impedance parameters based on phase characteristics and trend prediction; and generating a module control voltage correction value and a PWM control signal by applying a multi-time scale hybrid control strategy. Through hybrid control of phase space reconstruction of topology perception, multi-dimensional entropy analysis and stability guarantee, efficient suppression of the circulating current is realized, and the stability and robustness of the SVG system are improved.

Description

Technical Field

[0001] The present invention belongs to the field of power system control, and in particular to a method for suppressing circulating current of a modular high-voltage SVG. Background Art

[0002] With the expansion of the intelligentization of the power grid and the access scale of renewable energy, the high-voltage static var generator (SVG) has become the core equipment for reactive power compensation in modern power systems due to its advantages such as fast response, dynamic compensation, and harmonic suppression. In large-capacity application scenarios, SVG usually adopts a modular multilevel structure, and the capacity is expanded by connecting multiple power units in series and parallel to increase the voltage level. However, the circulating current problem caused by inconsistent parameters and switching characteristic differences between modules has become the key bottleneck restricting the improvement of the performance and reliability of high-voltage SVG, seriously affecting the safe operation of the system and the power quality. Therefore, the development of efficient circulating current suppression technology has important engineering value and research significance.

[0003] The current high-voltage SVG circulating current suppression technologies mainly include three categories: one is the method based on hardware improvement, such as adding a balancing reactor and a physical damping network, which restricts the circulating current path through an additional hardware circuit, but increases the system cost and loss; the second is the method based on improved PWM modulation, such as carrier phase-shifting technology and multi-carrier collaborative sorting strategy, which reduces the generation of circulating current by optimizing the switching timing, but has a limited application range; the third is the method based on virtual impedance, which simulates the impedance characteristics through a software algorithm and can suppress the circulating current without adding physical devices. Among them, the fixed virtual impedance method and the simple feedback regulation method are widely used, with the advantages of simple implementation and no increase in hardware cost.

[0004] However, the existing virtual impedance methods still have obvious deficiencies. Especially for high-voltage large-capacity SVG systems, with a large number of modules, complex topologies, and strong physical quantity coupling, it is difficult for traditional virtual impedance methods to achieve precise control. First, most traditional virtual impedance methods are designed based on simplified linear models, lacking in-depth analysis of the non-linear dynamic characteristics of SVG, resulting in limited suppression effects under complex working conditions; second, existing methods mostly adopt fixed parameters or simple adaptive regulation, and cannot effectively cope with the rapid changes in circulating current caused by power grid disturbances and load mutations. Summary of the Invention

[0005] The object of the invention is to provide a method for suppressing circulating current of a modular high-voltage SVG, in order to solve at least one technical problem existing in the prior art.

[0006] Technical solution: A method for suppressing circulating current of a modular high-voltage SVG, comprising: Obtain module data, perform preprocessing, and obtain preprocessed data, including circulating current data, voltage data, and status data; Construct an enhanced phase space based on the preprocessed data to generate feature data, including cyclic current phase features and dynamic trend prediction data; Based on the feature data and combined with the system physical constraints, calculate the virtual impedance parameters, including feedforward and feedback virtual impedance parameters; Use the virtual impedance parameters for gain modulation to generate a control signal and output control data, including the correction value of the module control voltage and the PWM control signal.

[0007] Advantageous effects: The present invention can capture the dynamic characteristics of the cyclic current more accurately, improving the cyclic current suppression efficiency; it can identify precursor signals before the cyclic current becomes unstable, initiate preventive control earlier, suppress the cyclic current peak below the safety threshold, and improve the reliability and robustness of the system. Description of the Drawings

[0008] Figure 1 It is a step flow chart of a modular high-voltage SVG cyclic current suppression method provided by an embodiment of the present application.

[0009] Figure 2 It is a step flow chart of generating feature data provided by an embodiment of the present application.

[0010] Figure 3 It is a step flow chart of generating reconstruction parameters provided by an embodiment of the present application.

[0011] Figure 4 It is a step flow chart of forming an enhanced state vector provided by an embodiment of the present application.

[0012] Figure 5 It is a step flow chart of generating cyclic current phase features provided by an embodiment of the present application. Detailed Embodiments

[0013] In order to enable those skilled in the art to better understand the solution 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 described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0014] It should be particularly noted that, for clearly showing the step flow of the present application, serial numbers are marked for each step in the specification. These serial numbers are only for the convenience of description and do not limit the execution order of the steps. In actual operation, according to the technical requirements of the specific implementation scenario, the steps can be executed in an order different from that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

[0015] It is found in the research that the existing control strategies are difficult to simultaneously consider the system dynamic response speed and the steady-state suppression accuracy, and it is easy to cause the instability of the control system when the module parameters vary greatly.

[0016] As Figure 1 shown, a modular high-voltage SVG circulating current suppression method is proposed, including the following steps: Obtain module data, preprocess it to obtain preprocessed data, including circulating current data, voltage data, and status data; Construct an enhanced phase space based on the preprocessed data to generate feature data, including circulating current phase features and dynamic trend prediction data; Based on the feature data, combined with the system physical constraints, calculate the virtual impedance parameters, including feedforward and feedback virtual impedance parameters; Use the virtual impedance parameters for gain modulation to generate control signals and output control data, including module control voltage correction values and PWM control signals.

[0017] Specifically, the module data includes module current data, module voltage data, switch status data, DC bus voltage data, and system temperature data. Noise filtering and data standardization are performed on the module data, and preprocessed circulating current data, preprocessed voltage data, and preprocessed status data are output. Using the preprocessed circulating current data, preprocessed voltage data, and preprocessed status data, an enhanced phase space based on the characteristics of the power electronic system is constructed, and circulating current phase features and dynamic trend prediction data are output. Based on the circulating current phase features and dynamic trend prediction data, combined with the system physical constraints, the optimal virtual impedance parameters are calculated, and the feedforward virtual impedance parameters and feedback virtual impedance parameters are output. Using the feedforward virtual impedance parameters and feedback virtual impedance parameters, adaptive gain modulation based on Lyapunov stability is realized to generate coordinated control signals for each module, and module control voltage correction values and PWM control signals are output.

[0018] According to one aspect of the present application, the preprocessed data is obtained, including: Collect the original current data, original voltage data, and switch status data from each module at a set sampling frequency (above 10 kHz), and simultaneously obtain the DC bus voltage data and temperature data.

[0019] Process the original current data using the switch state perception filtering technology, and obtain the filtered current data by identifying the switch transients and eliminating the corresponding noise; this technology adaptively removes the switch noise according to the switch time and state change characteristics, rather than using the traditional fixed parameter filtering; Analyze the current distribution among the modules, calculate the deviation of the current of each module from the average current, and extract the circulating current component; Compensate and correct the circulating current component according to the module characteristic differences (including gate voltage, junction temperature, and DC voltage) to obtain a more accurate compensated circulating current; consider the influence of actual hardware parameter differences on the circulating current; Apply a filtering algorithm adapted to the characteristics of the power electronic system to the original voltage data for denoising and normalization processing, and output the preprocessed voltage data; Sort and format the switch state data, construct a standard data structure including the relationship between the switch state and time, and form the preprocessed state data for subsequent analysis and use.

[0020] Specifically, collect the original current data Iraw, the original voltage data Vraw, and the switch state data Sraw from each module at a set sampling frequency fs (above 10 kHz), and simultaneously obtain the DC bus voltage data Vdc and the temperature data Tj. Apply switch state-aware filtering processing to the original current data Iraw to eliminate switch transient noise and obtain the filtered current data If: If(t) = Iraw(t) - Σ α_k·fsw(t-tk, Sraw(tk), Sraw(tk+1)); where tk is the switch time, fsw is the switch transient model function, α_k is the adaptive weight, and t is the current time. Calculate the current difference between modules and extract the circulating current component: Icirc,i(t) = Ii(t) - (1 / N)·Σ Ij(t); where N is the total number of parallel modules, and Ii(t) is the current of the i-th module. Compensate the circulating current component Icirc based on the module parameter differences to obtain the compensated circulating current: Icirc_comp,i(t)=Icirc,i(t)·(1+ηi(Vgs,i,Tj,i,Vdc)); where ηi is the compensation coefficient considering the gate voltage Vgs,i, the junction temperature Tj,i, and the DC voltage Vdc. Filter and normalize the original voltage data Vraw to obtain the preprocessed voltage data Vp. Convert the switch state data Sraw into a standard format to form the preprocessed state data Sp, including the specific switch state and time information of each module.

[0021] As Figure 2 shown, according to one aspect of the present application, generate feature data, including: Analyze the state data, split the circulating current data into a predetermined number of subsequences according to the switch state, and form the segmented circulating current data; Read the system topology and switch state, dynamically calculate the time delay and embedding dimension for the subsequences, and generate the reconstruction parameters; Extract key feature points from the circulating current data, and construct an enhanced state space in combination with the reconstruction parameters. Integrate the time-delay sequence in the circulating current data with physical quantities to form an enhanced state vector. Based on the enhanced state vector, calculate the geometric characteristics of the phase trajectory to generate the circulating current phase characteristics. Evaluate the coupling relationship of the circulating current between modules to form a module coupling matrix, and generate dynamic trend prediction data in combination with the circulating current phase characteristics.

[0022] Specifically, analyze the preprocessed state data, identify different switch combination modes in the system, and divide the compensated circulating current into multiple subsequences according to the switch states to form segmented circulating current data. This segmentation method is different from the traditional fixed-window segmentation. It is adaptively segmented based on the actual working state of the system. According to the current SVG system topology and switch states, dynamically calculate the optimal time-delay parameter and embedding dimension for each circulating current subsequence. The calculation process takes into account factors such as equivalent inductance, equivalent resistance, the number of modules, and the number of active modules to generate a set of reconstruction parameters. This dynamic parameter selection method overcomes the limitations of traditional fixed-parameter reconstruction. Implement the variable-length adaptive sliding window technique to determine the minimum necessary analysis time window according to the real-time characteristics of the circulating current, and generate adaptive window data. This method can dynamically adjust the size of the analysis window according to the periodicity and change rate of the circulating current, improving the calculation efficiency. Intelligently extract key turning points, extreme points, and rapidly changing points from the compensated circulating current instead of using all sampling points to form a set of key feature points. This selective sampling method reduces the computational amount while retaining the key dynamic characteristics of the circulating current. Construct an enhanced state space by combining the time-delay sequence of the compensated circulating current with key power electronics physical quantities (such as voltage, current change rate, and temperature difference) to form an enhanced state vector. This state space with multi-physical quantity fusion goes beyond the traditional single-variable reconstruction method and captures the system characteristics more comprehensively. Calculate the local geometric characteristics of the phase trajectory, including trajectory curvature, local stability index, and phase space velocity change, and generate the circulating current phase characteristics in combination with the energy flow characteristics and parameter sensitivity. These characteristic quantities directly reflect the dynamic essence of the circulating current rather than only focusing on time-domain or frequency-domain characteristics. Evaluate the mutual relationship between the circulating currents of different modules, calculate the coupling strength and phase difference between modules to form a module coupling matrix, which characterizes the dynamic interaction characteristics between modules. Based on the circulating current phase characteristics and the module coupling matrix, use recursive prediction technology to predict the trends and change characteristics of the circulating current at multiple future time points to generate dynamic trend prediction data. This prediction technology not only considers historical data but also integrates the physical characteristics of the system, improving the prediction accuracy.

[0023] In one embodiment of the present application, the compensation circulating current Icirc_comp is segmented into a set of subsequences Icirc_seg related to the switching state according to the preprocessing status data Sp: for time t, define the switching state matrix Sp(t); classify time points with the same or similar switching states: Tk = {t|Sp(t)≈Sk}; generate subsequences for each group of time points: Icirc_seg,k = {Icirc_comp(t)|t∈Tk}. Based on the SVG topology and the current switching state, calculate the optimal time delay τk and embedding dimension mk for each subsequence: τk = fτ(Leq(Sk),Req(Sk),freq_char(Icirc_seg,k)); mk = fm(N_module,N_active(Sk), complexity(Icirc_seg,k)); where Leq and Req are the equivalent inductance and resistance calculated based on the switching state, freq_char is the frequency characteristic function, N_module is the total number of system modules, N_active is the number of currently working modules, complexity is the signal complexity evaluation function, fτ and fm are functions for calculating the time delay and embedding dimension respectively, Sk is the switching state; k is the state index. Adopt the variable-length sequence adaptive sliding window technique to determine the minimum necessary time window Lwindow: Lwindow(t)=max{Lmin,β·fperiod(Icirc_comp(t),Sp(t))}; where Lmin is the minimum window length, fperiod estimates the main period of the current circulating current, and β is the safety factor. Extract the set of key feature points Kpoints from Icirc_comp instead of using all sampling points: Kpoints={t| there exists ε,|Icirc_comp(t)-Icirc_comp(t±ε)|>Δ(t) or |dIcirc_comp(t) / dt|>vth(t)}; where Δ(t) is the adaptive threshold, vth(t) is the change rate threshold, ε is the time offset, and dIcirc_comp(t) / dt represents the instantaneous change rate of the circulating current.

[0024] Construct the enhanced state space vector X(t) and add key power electronics physical quantities: X(t) = [Icirc_comp(t), Icirc_comp(t-τk),..., Icirc_comp(t-(mk-1)τk), Vp(t), dIcirc_comp(t) / dt, ΔTj(t)]; where Icirc_comp(t) is the compensated circulating current; τk is the time delay; Vp(t) is the preprocessed voltage, dIcirc_comp(t) / dt is the rate of change of the circulating current, ΔTj(t) is the junction temperature difference of the key power device, t is the time; k is the state index. Calculate the local geometric features of the circulating current phase trajectory to generate the circulating current phase feature CFPS: Trajectory curvature: κ(t) = |d 2 X(t) / dt 2 | / |dX(t) / dt| 3 ; Local Lyapunov exponent: λL(t) = (1 / Δt)·ln(|ΔX(t+Δt)| / |ΔX(t)|); Phase space velocity vector: v(t) = dX(t) / dt; Circulating current phase feature CFPS = [κ(t), λL(t), v(t), energy_flow(t), param_sensitivity(t)]; where energy_flow(t) is the energy flow feature, param_sensitivity(t) is the parameter sensitivity vector, Δt is the time interval, and ΔX is the small perturbation in the phase space. Calculate the inter-module coupling strength matrix C(t): C_{i,j}(t) = ∫Icirc_comp,i(τ)·Icirc_comp,j(τ+t)dτ / sqrt[∫Icirc_comp,i 2 (τ)dτ·∫Icirc_comp,j 2\(\int_{0}^{t}(\tau)d\tau\]; \(C(t)=[C_{i,j}(t)]_{N\times N}\). Combining the cyclic current phase feature CFPS and the coupling strength matrix \(C(t)\), applying recursive prediction technology to predict the future cyclic current trend, and generating dynamic trend prediction data DTP: Predicting the state at the next moment: \(X(t + \Delta t)=F(X(t), CFPS, C(t))\); Predicting the state at multiple steps: \(X(t + k\Delta t)=G(X(t+(k - 1)\Delta t), X(t+(k - 2)\Delta t),\cdots, X(t))\); Extracting the predicted cyclic current: \(I_{circ\_pred}(t + k\Delta t)=extract\_I(X(t + k\Delta t))\); where \(F()\) represents the state function for predicting the next moment, \(G()\) represents the state function for predicting multiple steps, and \(extract\_I()\) represents the function for extracting the cyclic current. The dynamic trend prediction data \(DTP = [I_{circ\_pred}(t+\Delta t), I_{circ\_pred}(t + 2\Delta t),\cdots, I_{circ\_pred}(t + n\Delta t), trend\_indicators(t)]\); where \(trend\_indicators(t)\) contains the key indicators of the cyclic current change trend.

[0025] As Figure 3 shown, according to one aspect of the present application, generating reconstruction parameters includes: Extracting the topological connection matrix at the current switching moment from the state data, and calculating the equivalent circuit parameters under the current topology based on the topological connection matrix; Analyzing the frequency characteristics of the segmented cyclic current data, determining the main frequency component and bandwidth, and calculating the optimal time delay parameter in combination with the equivalent circuit parameters; Evaluating the complexity of the segmented cyclic current data, and determining the optimal embedding dimension in combination with the total number of system modules and the number of currently active modules; Associating the optimal time delay parameter and the embedding dimension with the topological connection matrix to form reconstruction parameters.

[0026] Specifically, read the preprocessing status data Sp, extract the topological connection matrix Tconn at the current switching moment, and record the connection relationship and switching state of each module with the system. This matrix encodes the four switching states of each H-bridge module to form a topological description. Calculate the equivalent circuit parameters under the current topology based on Tconn, including the module equivalent inductance Leq, equivalent resistance Req, and equivalent capacitance Ceq. The calculation process considers the series and parallel relationships between modules. Different from the traditional fixed-parameter method, these parameters change in real time with the switching state. Analyze the frequency characteristics of the segmented circulating current data Icirc_seg, and calculate the main frequency component f_main, bandwidth BW, and harmonic distribution HarmDist. Different from the traditional FFT analysis, here a short-time adaptive spectrum estimation is used to adapt to the non-stationary characteristics of the circulating current. According to the information theory principle, calculate the optimal time delay parameter τk based on the main frequency component and bandwidth to ensure the maximum independence of data points in the phase space: For high-frequency dominant circulating current (f_main > f_threshold): Set a shorter delay to capture the rapid changes; For low-frequency dominant circulating current (f_main ≤ f_threshold): Set a longer delay to capture the complete cycle characteristics; Dynamically adjust the delay parameter: τk = α·(1 / f_main)·(Leq / Leq_ref) β , where α and β are parameters preset according to the system response characteristics, and f_threshold is the main frequency component threshold.

[0027] Evaluate the complexity of the segmented circulating current data \(I_{circ\_seg}\), calculate the complexity index \(CI\), including waveform complexity, irregularity, and entropy value. Different from the traditional signal complexity evaluation, here special attention is paid to the impact of switch jumps and state transitions in the power electronic system on the complexity. According to the total number of system modules \(N_{module}\), the current active modules \(N_{active}\), and the complexity index \(CI\), determine the optimal embedding dimension \(m_k\): Basic dimension determination: \(m_{base}=N_{active}+1\), which considers the degrees of freedom between modules; Complexity compensation: \(m_{comp}=\lceil CI\cdot\gamma\rceil\), where \(\gamma\) is the compensation coefficient, and \(\lceil\rceil\) represents rounding up; Final embedding dimension: \(m_k = \min(m_{base}+m_{comp},m_{max})\), ensuring that it does not exceed the maximum allowable dimension \(m_{max}\). Associate the calculated \(\tau_k\) and \(m_k\) with the current topological state \(T_{conn}\), and store them as the reconstruction parameter set \(RP\), including time delay, embedding dimension, and the corresponding topological state identifier. This parameter set will be used for subsequent phase space construction. Dynamically evaluate the performance of the reconstruction parameter set \(RP\), calculate the parameter sensitivity index \(PS\), and perform fine-tuning if necessary to optimize the reconstruction effect. This adaptive optimization mechanism is not available in traditional reconstruction methods. Feedback optimization based on the parameter sensitivity index \(PS\): Regularly evaluate the change trend of the parameter sensitivity index \(PS\). If \(PS\) continues to increase, it means that the current reconstruction parameters are gradually not suitable for the system changes, triggering parameter recalculation: If \(PS>\alpha\cdot PS_{previous}\) and lasts for \(N\) cycles, then: Re-execute the above steps to update the reconstruction parameter set \(RP\); Use the weighted average of the old and new parameters for smooth transition; Output the optimized reconstruction parameter set \(RP_{opt}\) for subsequent phase space construction. Where \(PS_{previous}\) is the parameter sensitivity index of the previous cycle.

[0028] As Figure 4 shown, according to one aspect of the present application, form an enhanced state vector, including: Construct a basic time-delay embedding vector according to the circulating current data and the reconstruction parameters; Extract key voltage features from the voltage data; Calculate the dynamic characteristics of the circulating current data; Based on the pre-stored system temperature data, calculate the temperature characteristics of the key power devices, including the junction temperature difference and the temperature change rate; Weightedly fuse the basic time-delay embedding vector with the key voltage features, dynamic characteristics, and temperature characteristics through a multi-physical quantity weight fusion matrix to form an enhanced state vector.

[0029] Specifically, read the preprocessed cyclic current data Icirc_comp and the reconstruction parameter set RP, and select the corresponding time delay τk and embedding dimension mk according to the current topological state. Construct the basic time-delay embedding vector Xbase, which contains the current moment and historical cyclic current data: extract the values at the current moment t and the past (mk - 1) delayed moments from Icirc_comp to form the vector Xbase(t) = [Icirc_comp(t), Icirc_comp(t - τk),..., Icirc_comp(t - (mk - 1)τk)]; this basic vector captures the time-evolution characteristics of the cyclic current. Read the preprocessed voltage data Vp and extract key voltage features, including the module output voltage Vout, the voltage change rate dV / dt, and the voltage unbalance Vunbal. Different from traditional methods, special attention is paid here to the correlation between voltage features and cyclic current. Calculate the dynamic characteristic indexes of the preprocessed cyclic current data Icirc_comp, including the change rate dI / dt, the acceleration d 2 I / dt 2 and the change trend Itrend. These derivative characteristics reflect the dynamic behavior of the cyclic current, which is often ignored in the traditional state space. Read the system temperature data and calculate the junction temperature difference ΔTj and the temperature change rate dTj / dt between key power devices. Incorporate the electro-thermal coupling effect into the state space and extract the dynamic characteristic parameters of the switching device, including the switching delay Tdelay, the rise / fall time Tr / Tf, and the on-state voltage Von. It reflects the actual switching behavior of the power device and affects the cyclic current distribution.

[0030] Construct the multi-physical quantity weight fusion matrix W_fusion, and adaptively adjust the weights of each physical quantity according to the current working condition: Under steady-state working conditions: increase the weights of temperature and voltage unbalance; Under dynamic working conditions: increase the weights of current change rate and acceleration; Under disturbance working conditions: increase the weights of switching characteristic parameters. This adaptive weight allocation mechanism is not available in traditional state space methods. Perform weighted fusion on the basic time-delay embedding vector Xbase, voltage features, current dynamic characteristics, temperature characteristics, and switching characteristics to construct the enhanced state space vector X(t): Read Xbase, Vout, dV / dt, Vunbal, dI / dt, d 2 I / dt 2 、Itrend、ΔTj、dTj / dt、Tdelay、Tr / Tf and Von; Apply W_fusion for multi-physical quantity weighted fusion and output the complete enhanced state vector X(t) as the basis for subsequent feature extraction and prediction. Normalize the enhanced state vector X(t) to ensure the dimensional consistency of different physical quantities and output the standardized enhanced state vector Xnorm(t). This normalization takes into account the typical ranges and dynamic characteristics of different physical quantities.

[0031] As Figure 5 shown, according to one aspect of the present application, generating a cyclic current phase feature includes: Reading and based on the enhanced state vector, in sequence: calculating the trajectory curvature and twist, forming the trajectory geometric feature; evaluating the local stability of the phase trajectory, forming the stability index; calculating the phase space velocity vector and acceleration vector, forming the dynamic evolution feature; Based on the cyclic current data and state data, forming the energy flow feature; Extracting the sensitivity of the cyclic current data to the module parameter differences, forming the parameter sensitivity feature; Determining the weights of the trajectory geometric feature, stability index, dynamic evolution feature, energy flow feature and parameter sensitivity feature, forming the cyclic current phase feature.

[0032] Specifically, reading a series of standardized enhanced state vectors Xnorm(t), constructing the phase trajectory Xtraj, which contains the state evolution information of consecutive time points. This trajectory captures the dynamic behavior of the system in the enhanced state space. Calculating the local geometric features of the phase trajectory, analyzing the shape and curvature of the trajectory: reading a sequence of consecutive points in the phase trajectory, calculating the trajectory curvature κ(t), characterizing the curvature of the trajectory and reflecting the non-linearity of the system state change; calculating the trajectory twist τ(t), characterizing the twist of the trajectory in the high-dimensional space; outputting the trajectory geometric feature GF, these features go beyond the traditional linear analysis methods and can capture the non-linear dynamics of the system. Evaluating the local stability of the phase trajectory, quantifying the sensitivity of the system to perturbations: selecting a reference point and its neighborhood points in the phase trajectory, tracking the evolution trajectories of these points, calculating the trajectory separation index; extracting the local Lyapunov exponent λL(t), characterizing the local stability of the system; outputting the stability index SI, these indexes can predict the behavior of the system under perturbations and provide forward-looking information for control.

[0033] Calculating the time evolution characteristics of the phase trajectory, capturing the rate of change of the system state: reading the states of consecutive time points in the phase trajectory, calculating the phase space velocity vector v(t)=dX(t) / dt, characterizing the rate of state change; calculating the phase space acceleration vector a(t)=d 2 X(t) / dt 2, characterizing the acceleration characteristics of state changes; outputting the dynamic evolution feature DF, which reflects the temporal evolution pattern of the system state. Analyze the periodicity and recurrence of the phase trajectory to identify regular patterns in the system: construct a delay coordinate map and analyze the self-similarity of the phase trajectory; calculate recurrence plot metrics to quantify the recurrence characteristics of the trajectory; extract the periodicity intensity PI and cyclic pattern CP; output the periodic feature PF, which helps to identify the regular components in the circulating current. Calculate the energy flow characteristics of the circulating current within the switching period: read the preprocessed circulating current data Icirc_comp and preprocessed state data Sp, divide the energy analysis window according to the switching period, and calculate the circulating current energy Ecirc(t) and power flow Pcirc(t); determine the energy flow direction index; output the energy flow feature EF, which establishes the connection between phase dynamics and physical energy. Extract the sensitivity of the circulating current to module parameter differences and analyze the impact of parameter changes on the system: construct a parameter change response model based on the enhanced state vector X(t) and calculate the sensitivities to resistance difference SR(t), inductance difference SL(t), and capacitance difference SC(t); comprehensively form the parameter sensitivity vector S(t) and output the parameter sensitivity feature SF, which reflects the response characteristics of the circulating current to system parameter changes.

[0034] Integrate all the above features to construct the complete circulating current phase feature CFPS: read the trajectory geometric feature GF, stability index SI, dynamic evolution feature DF, periodic feature PF, energy flow feature EF, and parameter sensitivity feature SF; determine the weights of each feature through feature importance evaluation; synthesize the final circulating current phase feature CFPS, which contains comprehensive information on the dynamics of the circulating current; this feature set will be used for subsequent trend prediction and control parameter calculation.

[0035] According to one aspect of the present application, generate dynamic trend prediction data, including: Based on the enhanced state vector and the circulating current phase feature, select a set of points similar to the current state from historical data and construct a local linear prediction model; Based on the local linear prediction model, predict the state at the next moment with the current state, and then use the prediction result as the new initial condition to predict the more distant future and generate a multi-step prediction sequence; Extract the circulating current prediction data from the multi-step prediction sequence and calculate the key dynamic indicators of the predicted circulating current, and integrate them into the dynamic trend prediction data.

[0036] Specifically, read the enhanced state vector X(t) of the historical time series and the cyclic current phase feature CFPS to construct the initial conditions for prediction. Ensure that the prediction is based on the complete state information of the system. Construct a local linear prediction model in the phase space, and build a linear approximation near the current state: select a set of points Xsim similar to the current state from the historical data; based on the evolution trajectories of these points, construct a local linear mapping Flocal; output the local prediction model LPM, which adapts to the local linearization of the nonlinear characteristics of the cyclic current. Implement a recursive rolling prediction strategy to gradually extend the prediction range: read the current state X(t) and the local prediction model LPM; predict the state at the next moment X(t+Δt) = Flocal(X(t)); use the prediction result as the new initial condition to predict the more distant future X(t+2Δt); repeat this process to generate a multi-step prediction sequence Xpred; this recursive prediction method is superior to the traditional single-step prediction and can capture long-term trends. Introduce an error correction mechanism for the phase feature to improve the recursive prediction accuracy: analyze the systematic error patterns in the historical predictions; adjust the prediction correction factor CF based on the cyclic current phase feature CFPS; apply the correction factor to optimize the prediction result: Xpred_corr = Xpred + CF; output the corrected prediction sequence Xpred_corr, which reduces the error accumulation in the recursive prediction.

[0037] Implement multi-model prediction fusion to integrate the advantages of different prediction strategies: run multiple prediction models in parallel, including local linear models, nonlinear dynamics models, and pattern matching models; dynamically adjust the weights Wmodel of each model according to the historical accuracy and the current cyclic current phase feature CFPS; fuse the prediction results of each model with weights to form an integrated prediction Xpred_ensemble; this multi-model fusion method avoids the limitations of a single model and improves the prediction robustness. Extract the cyclic current trend information from the predicted state sequence: extract the cyclic current prediction component Icirc_pred from Xpred_ensemble; analyze the change trend of the predicted cyclic current, calculate the peak value Ipeak, the root mean square value Irms, and the frequency characteristics Ifreq; evaluate the reliability of the prediction and generate a credibility index Conf; output the cyclic current prediction data IP, which includes the amplitude and characteristic predictions of the future cyclic current. Calculate the key dynamic indicators of the predicted cyclic current, including the maximum change rate max(dI / dt), the energy change trend ΔE, and the stability expectation Stability, to form the trend indicator data TI. Integrate the cyclic current prediction data IP and the trend indicator data TI to form the complete dynamic trend prediction data DTP, providing forward-looking information for the subsequent virtual impedance parameter calculation.

[0038] According to one aspect of the present application, calculate the virtual impedance parameters, including: Analyze the energy distribution and flow characteristics of the circulating current data within the switching period to generate energy flow characteristic data; Evaluate the sensitivity of the circulating current data to the module parameter differences to form parameter sensitivity data; Calculate the sample entropy and entropy change rate of the enhanced phase space to analyze the complexity characteristics of the circulating current; Combine the complexity characteristics of the circulating current with the dynamic trend prediction data, energy flow characteristic data, and parameter sensitivity data to calculate the optimal virtual impedance value at the future moment and form the feedforward virtual impedance parameter; Construct an adaptive feedback control strategy based on the circulating current phase characteristics, calculate the virtual impedance value of the real-time response, and generate the feedback virtual impedance parameter.

[0039] Specifically, comprehensively analyze the preprocessed circulating current data, preprocessed voltage data, and module coupling matrix, calculate the operating condition index vector characterizing the current operating state of the system, including load change rate, grid fluctuation degree, active power, reactive power, harmonic distortion degree, and module coupling degree, to form the operating condition index data. Monitor the change trend of the operating condition index data. When a significant change is detected (such as load mutation or grid disturbance), mark it as an operating condition conversion point, trigger the corresponding adjustment of the control strategy, and output the operating condition conversion flag. This active operating condition identification method can respond to system state changes in advance. Analyze the energy distribution and flow characteristics of the circulating current within the switching period, calculate the circulating current energy and its change rate, determine the energy flow direction and intensity, and generate the energy flow characteristic data. These data reflect the energy dynamic characteristics of the circulating current and are an important basis for optimizing the virtual impedance. Evaluate the sensitivity of the circulating current to the parameter differences between modules, calculate the sensitivities to resistance differences, inductance differences, and capacitance differences respectively, to form the parameter sensitivity data. This sensitivity analysis helps to identify the main influencing factors of the circulating current and optimize the control parameters targeted. Calculate the complexity index of the circulating current phase space, including the sample entropy value and entropy change rate, and extract the signal complexity characteristics. Reflects the non-stationarity and uncertainty of the circulating current and guides the dynamic adjustment of the control strategy. According to the dynamic trend prediction data, energy flow characteristic data, parameter sensitivity data, and signal complexity characteristics, comprehensively calculate the optimal virtual impedance value at the future moment to form the feedforward virtual impedance parameter. This parameter is optimized for the predicted circulating current trend and has a preventive control effect. Design an adaptive feedback control strategy based on the current circulating current phase characteristics, calculate the virtual impedance value of the real-time response, and generate the feedback virtual impedance parameter. This parameter is responsible for the immediate adjustment of the current circulating current state to ensure the accuracy of the control. Considering the physical limitations of power electronic devices, impose safety constraints on the virtual impedance parameters to ensure that the parameter values and change rates are within the tolerable range of the system, and output the constrained virtual impedance parameters. Ensures the realizability of the control parameters and the safe operation of the system.

[0040] In one embodiment of the present application, based on the preprocessed circulating current data, preprocessed voltage data, and the inter-module coupling strength matrix C(t), the system operating condition index vector W(t) is calculated: W(t) = [|dIload / dt|, |dVgrid / dt|, Pactive, Qreactive, THD, C_norm(t)]; where C_norm(t) is the norm of the coupling matrix, characterizing the degree of coupling between modules. Detect significant changes in the operating conditions and trigger adjustments to the control strategy if necessary: if ||W(t) - W(t-Δt)|| > θW·||W(t-Δt)||: Mark it as a point of operating condition transition and adjust the calculation strategy of the control parameters; where θW is the operating condition change threshold. Calculate the circulating current energy feature Ecirc(t) and the power flow feature Pcirc(t) within the switching period: Ecirc(t) = ∫ t t+Tsw L(τ)·Icirc_comp 2 (τ)dτ / 2; Pcirc(t) = dEcirc(t) / dt; Energy flow direction index: Denergy(t) = sign(Pcirc(t))·|Pcirc(t)| γ; where Tsw is the switching period, L(τ) is the time-varying equivalent inductance, and γ is the energy exponential parameter. Extract the sensitivity vector S(t) of the circulating current to the module parameter differences: SR(t) = ΨIcirc_comp(t) / ΨΔR; SL(t)=ΨIcirc_comp(t) / ΨΔL; SC(t)=ΨIcirc_comp(t) / ΨΔC; S(t)=[SR(t),SL(t),SC(t)]; where Ψ is the partial derivative. Calculate the sample entropy SampEn(t) and the entropy change rate dSampEn(t) of the circulating current phase space: SampEn(t)=-ln[Am(r) / Bm(r)]; dSampEn(t)= [SampEn(t)-SampEn(t-Δt)] / Δt; where Am(r) and Bm(r) are the ratios of the number of point pairs with distances less than r in the m+1 and m-dimensional phase spaces, respectively. Based on the dynamic trend prediction data DTP, the energy feature Ecirc(t), the sensitivity vector S(t), and the entropy change rate dSampEn(t), calculate the feedforward virtual impedance parameter Zv_ff: At the prediction time t+kΔt: Zv_ff(t+kΔt) = h(Icirc_pred(t+kΔt),Ecirc(t),S(t),dSampEn(t),W(t)); where Icirc_pred is the predicted circulating current; Ecirc is the circulating current energy; S is the parameter sensitivity vector; dSampEn is the entropy change rate; W is the operating condition index vector; h is the mapping function; t is the time; k is the prediction step; Δt is the time interval; h is the mapping function that maps the predicted circulating current characteristics to the optimal virtual impedance value. Based on the current circulating current phase feature CFPS and the system state, calculate the feedback virtual impedance parameter Zv_fb: Zv_fb(t)=KP·Icirc_comp(t)+KI·∫Icirc_comp(τ)dτ + KD·dIcirc_comp(t) / dt; where KP is the proportionality coefficient; KI is the integral coefficient; KD is the differential coefficient; dynamically adjusted with CFPS and the operating condition index W(t); Icirc_comp is the compensated circulating current; t is the time; τ is the integration variable. Apply the physical limitations of the power electronic devices to ensure that the virtual impedance parameters satisfy the realizability constraints: |Zv_ff(t)|≤Zmax(t); |dZv_ff(t) / dt|≤η(t)·Vdc / (Irated·Tsw); where Zmax(t) is the maximum allowable impedance value, η(t) is the dynamic safety factor, Vdc is the DC bus voltage, and Irated is the rated current.

[0041] According to one aspect of the present application, the process of forming the operating condition index data is specifically as follows: Read the preprocessed circulating current data, calculate the circulating current amplitude, phase, frequency, and form factor to form the basic characteristics of the circulating current. Analyze the preprocessed voltage data, extract the change characteristics of the grid voltage, including the voltage change rate dVgrid / dt, fluctuation index, and unbalance degree, and output the grid state characteristics. Read the preprocessed state data and the module coupling matrix, and evaluate the coupling state between modules: Calculate the norm of the coupling matrix to characterize the overall coupling strength; extract the maximum coupling value and the corresponding module pair; analyze the uniformity of the coupling distribution and calculate the coupling imbalance index; output the module coupling characteristics, which reflect the interaction between the internal modules of the system. Calculate the system power characteristics, including active power Pactive, reactive power Qreactive, and fluctuating power Sosc: Read the preprocessed voltage data Vp and the preprocessed circulating current data Icirc_comp; calculate the instantaneous power S(t) = Vp(t)·Icirc_comp(t); decompose it into active, reactive, and fluctuating components; output the power characteristic data PD, which characterizes the energy transfer state of the system. Analyze the system harmonic state, calculate the total harmonic distortion THD and the main harmonic components Harm_n: Perform harmonic analysis on the preprocessed circulating current data; calculate the ratio of the amplitude of each harmonic to the fundamental wave; determine the dominant harmonics and their phase relationships; output the harmonic characteristic data HD, which characterizes the nonlinear distortion state of the system. Evaluate the system load characteristics, including the load change rate dIload / dt, load type index, and load stability: Analyze the load-related patterns in the preprocessed circulating current data; judge the load type through the power factor and harmonic characteristics; evaluate the regularity and predictability of the load change; output the load characteristic data LD, which characterizes the external load condition of the system. Construct an operating condition evaluation weight matrix W_eval, and dynamically adjust the weights of each characteristic according to the key points of system operation: In the high-precision circulating current control mode: Increase the weights of the circulating current characteristics and the module coupling characteristics; In the fast dynamic response mode: Increase the weights of the load characteristics and the grid state characteristics; In the harmonic suppression mode: Increase the weights of the harmonic characteristic data. This adaptive weight adjustment mechanism is not available in traditional operating condition evaluation methods. Integrate all characteristics to construct an operating condition index vector W(t): Read the basic characteristics of the circulating current IB, the grid state characteristics GS, the module coupling characteristics MC, the power characteristic data PD, the harmonic characteristic data HD, and the load characteristic data LD; perform weighted fusion using W_eval; output the final operating condition index data W(t), which comprehensively characterizes the current operating condition of the system; this index will be used for subsequent control parameter calculation and strategy adjustment.

[0042] According to one aspect of the present application, analyze the complexity characteristics of the circulating current, including: Read the circulating current data and perform the following operations: Segment the cyclic current data into subsequences of lengths m and m + 1, calculate the proportion of similar subsequences, determine the conditional probability, and obtain the sample entropy of the original signal; where m is a natural number greater than 0. Perform coarse-graining processing on the cyclic current data, calculate the sample entropy for each coarse-grained sequence, form an entropy-scale curve, and generate multi-scale entropy features. Based on the cyclic current data, analyze the sequence pattern and fuzzy boundary, and calculate the permutation entropy and fuzzy entropy. Calculate the change rates of the sample entropy of the original signal, multi-scale entropy features, permutation entropy, and fuzzy entropy, and perform weighted fusion on each entropy based on the change rates to generate the cyclic current complexity features.

[0043] Specifically, read the preprocessed cyclic current data Icirc_comp within a time window, ensuring that it contains sufficient dynamic information. The time window length is adaptively adjusted according to the system characteristics instead of the traditional fixed window length. Calculate the sample entropy to quantify the complexity and irregularity of the time series: Segment Icirc_comp into subsequences of lengths m and m + 1; calculate the proportion of similar subsequences and determine the conditional probability; calculate the sample entropy value SampEn = -ln(A / B), where A is the m + 1-dimensional similarity proportion and B is the m-dimensional similarity proportion; output the sample entropy value, which characterizes the randomness and unpredictability of the cyclic current. Implement multi-scale entropy analysis to evaluate the complexity at different time scales: Apply coarse-graining processing with different scale factors τ to Icirc_comp; calculate the sample entropy for each coarse-grained sequence; generate a multi-scale entropy curve; output the multi-scale entropy feature MSE, which characterizes the complex dynamics of the cyclic current at different time scales. Calculate the permutation entropy to capture the sequence pattern of the time series: Divide Icirc_comp into overlapping subsequences; analyze the relative magnitude order of the data points in each subsequence; count the occurrence frequencies of different permutation patterns; calculate the permutation entropy PermEn = -Σ p(π)·log(p(π)), where p(π) is the probability of the permutation pattern π; output the permutation entropy value PermEn, which characterizes the sequence complexity of the cyclic current.

[0044] Perform fuzzy entropy analysis to process the noise and fuzzy boundaries in the circulating current: Apply the fuzzy membership function to Icirc_comp to calculate the sequence similarity; Calculate the fuzzy entropy value based on the fuzzy similarity; Output the fuzzy entropy value FuzzyEn to provide a more robust complexity assessment for the noise. Calculate the entropy change rate dEntropy to characterize the dynamic characteristics of the system state transition: Track the time variations of SampEn, MSE, PermEn, and FuzzyEn; Calculate the change rates of each entropy value: dSampEn, dMSE, dPermEn, and dFuzzyEn; Synthesize to form the entropy change rate vector dEntropy; Output the entropy change rate data dEntropy to characterize the dynamic change of the system complexity. Implement entropy gradient analysis to identify the system state transition points: Calculate the spatial gradient and time gradient of dEntropy; Identify the points with significant gradient changes and mark them as transition points; Conduct key analysis on the intervals near the transition points; Output the state transition data ST, including the positions and characteristics of the transition points. Integrate each entropy index to construct a comprehensive complexity assessment model: Design the entropy feature fusion weights and adjust the weights of each entropy index according to the characteristics of the current circulating current; Fuse SampEn, MSE, PermEn, FuzzyEn, and dEntropy; Generate the comprehensive complexity index CI; Output the signal complexity feature CI to comprehensively characterize the complex dynamic characteristics of the circulating current; for subsequent control parameter optimization.

[0045] According to one aspect of the present application, forming the feedforward virtual impedance parameter includes: Predict the dynamic characteristics of the circulating current based on the dynamic trend prediction data analysis to determine the key control moments; Based on the energy flow feature data, parameter sensitivity data, and circulating current complexity features, predict the energy flow characteristics, parameter sensitivity, and the change amount of the circulating current complexity at future moments to form the prediction feature data; For the determined key control moments, apply the prediction feature data and the pre-constructed impedance optimization model to solve for the optimal virtual impedance values within the prediction time window, generate the virtual impedance time series, and form the feedforward virtual impedance parameter.

[0046] Specifically, read the dynamic trend prediction data DTP and extract the predicted circulating current and trend indicators at multiple future moments. These forward-looking data are the basis for implementing preventive control. Analyze the dynamic characteristics of the predicted circulating current to determine the key control moments: identify the peak moment t_peak and the zero-crossing moment t_zero of the predicted circulating current; determine the moment with the maximum rate of change and the moment with the maximum energy; calculate the importance weight W_time of the critical moment; output the critical moment data KT to guide the timing optimization of the control parameters. Read the energy flow characteristic data EF and analyze the energy dynamic characteristics of the circulating current: evaluate the energy flow direction Denergy_pred at future moments; predict the energy peak Epeak_pred and the energy fluctuation characteristics Efluc_pred; calculate the sensitivity E_sensitivity of the energy change to the system impact; output the predicted energy characteristics PE to provide an energy-level basis for the calculation of control parameters. Read the parameter sensitivity data S(t) and predict the impact of parameter changes on the future circulating current: evaluate the sensitivity SR_pred of the future circulating current to the resistance difference; evaluate the sensitivity SL_pred of the future circulating current to the inductance difference; evaluate the sensitivity SC_pred of the future circulating current to the capacitance difference; output the predicted sensitivity data PS to guide the precise adjustment of the virtual impedance parameters.

[0047] Read the signal complexity characteristic CI and analyze the non-stationary characteristics of the circulating current: predict the entropy value change CI_pred at future moments; identify possible state transition points; design control strategies for different complexity intervals; output the predicted complexity characteristic PC to provide a basis for dealing with non-stationary circulating current. Construct a mapping model of the impact of virtual impedance parameters on the circulating current: based on historical data, establish a relationship model Z_I_model between the virtual impedance and the circulating current response; this model considers the non-linear mapping relationship under different working conditions and goes beyond the traditional linear impedance model; output the impedance mapping model ZIM to provide a basis for solving the optimal parameters. Construct a multi-objective optimization function to balance the circulating current suppression and dynamic response: define the circulating current suppression objective function J_circ = f(Icirc_pred); define the dynamic response objective function J_resp = g(dIcirc_pred / dt); construct the comprehensive objective function J = α·J_circ + (1-α)·J_resp, where α is a weight dynamically adjusted according to the working conditions; output the optimization objective function OF to guide the solution of the optimal parameters. Apply the prediction optimization algorithm to solve the optimal virtual impedance parameters at future moments: read the impedance mapping model ZIM and the optimization objective function OF; for each critical moment within the prediction time window, solve the virtual impedance value corresponding to the minimum value of the objective function; generate the virtual impedance time series Z_seq; apply a smoothing constraint to ensure the continuity of the impedance change; output the feed-forward virtual impedance parameter Zv_ff, which contains the optimal impedance values at multiple future moments.

[0048] According to one aspect of the present application, generating feedback virtual impedance parameters includes: Reading the phase characteristics of the circulating current and: Extracting trajectory geometric characteristics, stability indicators, and dynamic evolution characteristics from the phase characteristics of the circulating current to construct PID control parameters; Extracting parameter sensitivity characteristics and period characteristics from the phase characteristics of the circulating current, adjusting the PID control parameters and evaluating the stability indicators to form characteristic matching parameters, resonance suppression components, and stability enhancement components; generating a basic impedance value through the PID control parameters and the characteristic matching parameters; Dynamically adjusting a preconfigured state feedback gain matrix according to the phase characteristics of the circulating current to form a state feedback component; Performing weighted fusion on the basic impedance value, the state feedback component, the resonance suppression component, and the stability enhancement component to form feedback virtual impedance parameters.

[0049] Specifically, read the current loop current phase characteristic CFPS, extract the key characteristics required for real-time control, including the trajectory geometric characteristic GF, the stability index SI, and the dynamic evolution characteristic DF. Design an adaptive PID controller based on the phase characteristic to dynamically adjust the control parameters: adjust the proportional coefficient KP according to the trajectory curvature κ(t): increase KP when the curvature is large to improve the response speed, and decrease KP when the curvature is small to avoid over-regulation; adjust the integral coefficient KI according to the local Lyapunov exponent λL(t): increase KI when the stability is high to eliminate the steady-state error, and decrease KI when the stability is low to prevent integral saturation; adjust the differential coefficient KD according to the phase space velocity v(t): increase KD when the velocity is large to provide damping, and decrease KD when the velocity is small to avoid noise amplification; output the adaptive PID parameters KPID, and this parameter adjustment mechanism based on the phase characteristic is different from the traditional error-driven parameter tuning. Read the parameter sensitivity characteristic SF and the period characteristic PF in the loop current phase characteristic CFPS, and fine-tune the control parameters for different loop current characteristics: focus on optimizing the high-sensitivity parameters: when the inductance sensitivity SL(t) is greater than the threshold, adjust the ratio of KP and KD to improve the inductance control accuracy; adjust the control frequency according to the periodicity intensity PI: synchronize the control update frequency with the loop current period when the periodicity is strong; preset the control parameters based on the loop pattern CP: apply the pre-optimized parameter combination when the typical pattern is recognized; output the feature matching parameters FMP to achieve more refined parameter adjustment. Construct the feedback virtual impedance base value to achieve an instant response to the current loop current: basic PID control law: Zv_base = KPID(0)·Icirc_comp(t) + KPID(1)·∫Icirc_comp(τ)dτ + KPID(2)·dIcirc_comp(t) / dt; fine-tune using the feature matching parameters FMP: Zv_base = Zv_base·FMP; output the base impedance value Zv_base as the basic component of the feedback control.

[0050] Read the enhanced state vector X(t), design a state feedback enhancer to supplement the basic PID control: extract a subset of state variables Xsub that is directly related to the circulating current; calculate the state feedback control quantity: Zv_state = K_state · Xsub, where K_state is the state feedback gain matrix; dynamically adjust K_state according to the circulating current phase characteristic CFPS to achieve adaptive state feedback; output the state feedback component Zv_state to enhance the direct response ability to the system state. Introduce a resonance suppression compensator to provide additional damping for the resonant components in the circulating current: analyze the periodic feature PF in the circulating current phase characteristic CFPS to identify the potential resonant frequency f_res; design a targeted resonance suppressor: Zv_res = A_res sin(2πf_res·t + φ_res); dynamically adjust the amplitude A_res and phase φ_res of the suppressor to form a 180° phase difference with the resonant component of the circulating current; output the resonance suppression component Zv_res to effectively suppress the oscillation of the circulating current at a specific frequency. Based on the stability index SI in the circulating current phase characteristic CFPS, design a stability enhancer: when the local Lyapunov exponent λL(t) approaches the critical value, add an additional stabilization component; calculate the stabilization control quantity: Zv_stab = K_stab · (λL_threshold - λL(t)) when λL(t)>λL_threshold; output the stability enhancement component Zv_stab to prevent the system from becoming unstable in the critical state. Integrate each control component to generate the final feedback virtual impedance parameter: read the basic impedance value Zv_base, the state feedback component Zv_state, the resonance suppression component Zv_res, and the stability enhancement component Zv_stab; construct a component weighted fusion strategy: dynamically adjust the weights W_fb of each component based on the circulating current phase characteristic CFPS; calculate the final feedback virtual impedance: Zv_fb = W_fb(0)·Zv_base + W_fb(1)·Zv_state + W_fb(2)·Zv_res + W_fb(3)·Zv_stab; output the feedback virtual impedance parameter Zv_fb to achieve precise response control of the current circulating current state.

[0051] According to one aspect of the present application, output control data, including: Construct a Lyapunov energy function based on the circulating current and control weights to evaluate the system stability and control effect, and generate stability evaluation data; this energy-based stability evaluation method can theoretically ensure the stable operation of the system; Construct a stability-oriented gain modulation mechanism to dynamically adjust the control weights according to the circulating current state and virtual impedance difference, ensure that the system evolves in a stable direction, and output the basic hybrid weights; Implement a multi-time-scale control strategy, design control layers with different response speeds for high-frequency disturbances, medium-frequency load changes, and low-frequency parameter drifts respectively, calculate the weights of each layer and synthesize them into the final hybrid weight to form a multi-scale hybrid weight; this hierarchical control structure can cope with system changes at different frequencies simultaneously; Based on the multi-scale hybrid weight, weight and fuse the feedforward virtual impedance parameter and the feedback virtual impedance parameter, and calculate the final comprehensive virtual impedance value; this feedforward-feedback hybrid control structure combines the advantages of predictive control and real-time feedback; According to the comprehensive virtual impedance value, calculate the voltage compensation amount that each module needs to apply to form a preliminary module voltage correction value; Considering the power balance and thermal stress distribution among modules, coordinate and optimize the module voltage correction value to ensure that each module maintains a reasonable power distribution while suppressing the circulating current and outputs a coordinated voltage correction value; Combine the coordinated voltage correction value with the reference voltage to generate the final control voltage command for each module and form a module control voltage correction value; Based on the module control voltage correction value, apply an appropriate PWM modulation strategy, considering factors such as switching frequency and dead time, to generate a PWM control signal for driving the power device.

[0052] Specifically, construct a Lyapunov energy function based on the circulating current and control weights to evaluate the system stability and control effect: Define the system state vector X(t)=[Icirc_comp,ω], which includes the compensated circulating current and the feedforward-feedback hybrid weight; construct the Lyapunov candidate function: V(X) = (1 / 2)L·Icirc_comp T ·Icirc_comp + (1 / 2)kω·(ω-ω*) T ·Q·(ω-ω*); where V is the Lyapunov energy function, L is the system equivalent inductance matrix, reflecting the electromagnetic energy of the circulating current among modules; ω is the feedforward-feedback hybrid weight vector, ω* is the ideal weight vector; kω is the weight adjustment coefficient matrix, Q is a positive definite weight matrix, reflecting the energy cost of the weight deviation; T represents the transpose of a vector or matrix. Calculate the time derivative of the Lyapunov function: dV / dt = Icirc_comp T ·L·(dIcirc_comp / dt)+kω·(ω-ω*) T·Q·(dω / dt). Construct a stability criterion to ensure dV / dt < 0: Construct a stability index SI = -dV / dt. When SI is positive and greater than the threshold, the system evolves towards stability; output stability evaluation data SE = {V(t), dV / dt, SI} as the direct basis for weight modulation. Construct a gain modulation mechanism based on Lyapunov stability to ensure the stable convergence of the system: Read the stability evaluation data. According to the Lyapunov stability theory, design a weight dynamic adjustment law to ensure dV / dt < 0: dω / dt = -γ·▽ω(dV / dt); where γ is a positive learning rate parameter that is dynamically adjusted to accelerate convergence; ▽ω(dV / dt) is the gradient of dV / dt with respect to ω, pointing in the direction where dV / dt increases fastest. Expand the stability-oriented weight adjustment equation: dω / dt = -γ·[kω·Q·(ω - ω*) + Icirc_comp T ·L·(Ψ(dIcirc_comp / dt) / Ψω)]; By analyzing the effects of Zv_fb and Zv_ff on the circulating current, it is further simplified to: dω / dt = -γ·sign(Icirc_comp·d(Zv_fb - Zv_ff) / dt)·|Icirc_comp|·|Zv_fb - Zv_ff| α ; where α is an exponential parameter that is dynamically adjusted according to the circulating current state. Apply the stability constraint to increase the damping term when the stability index SI is below the threshold: dω / dt = dω / dt - β·(ω - ω*), where dω / dt is the weight change rate; γ is the learning rate parameter; Icirc_comp is the compensated circulating current; α is the exponential parameter; β is the damping coefficient, and output the basic hybrid weight ω(t) as the initial weight for multi-time scale control.

[0053] Implement a multi-time scale control strategy to perform hierarchical weight optimization based on the stability evaluation data SE: Read the basic hybrid weight and the stability evaluation data; Decompose the circulating current into high-frequency, medium-frequency, and low-frequency components: Icirc_comp(t) = Icirc_hf(t) + Icirc_mf(t) + Icirc_lf(t); Correspondingly construct hierarchical Lyapunov functions: V_hf(t) = (1 / 2)L·Icirc_hf T ·Icirc_hf + (1 / 2)kω_hf·(ω_hf - ω_hf*) T ·Q_hf·(ω_hf - ω_hf*); V_mf(t) = (1 / 2)L·Icirc_mf T ·Icirc_mf + (1 / 2)kω_mf·(ω_mf - ω_mf*) T·Q_mf·(ω_mf - ω_mf*); V_lf(t) = (1 / 2)L·Icirc_lf T ·Icirc_lf + (1 / 2)kω_lf·(ω_lf - ω_lf*) T ·Q_lf·(ω_lf - ω_lf*); Calculate the optimal weights for high - frequency perturbations, medium - frequency load changes, and low - frequency parameter drifts respectively: dω_hf / dt = -γ_hf·▽ω_hf(dV_hf / dt); dω_mf / dt = -γ_mf·▽ω_mf(dV_mf / dt); dω_lf / dt = -γ_lf·▽ω_lf(dV_lf / dt); where γ_hf > γ_mf > γ_lf to ensure the fastest high - frequency response speed. Combine the weights of each frequency band to form a multi - scale hybrid weight: ω_ms(t) = λ_hf(t)·ω_hf(t) + λ_mf(t)·ω_mf(t) + λ_lf(t)·ω_lf(t); where λ_hf(t), λ_mf(t), and λ_lf(t) are dynamic frequency - band weights satisfying λ_hf(t)+λ_mf(t)+λ_lf(t)=1; ω_hf, ω_mf, ω_lf are the optimal weights of each frequency band. Verify the stability of the global Lyapunov function under the combined weight: dV_total / dt = dV_hf / dt + dV_mf / dt + dV_lf / dt < 0; Output the final multi - scale hybrid weight ω_ms(t) for feed - forward and feedback virtual impedance fusion. Based on the multi - scale hybrid weight, fuse the feed - forward and feedback virtual impedances to ensure global stability: Read the multi - scale hybrid weight ω_ms(t), the feed - forward virtual impedance parameter Zv_ff, and the feedback virtual impedance parameter Zv_fb; Apply the weights for weighted fusion: Zv(t) = ω_ms(t)·Zv_fb(t)+(1 - ω_ms(t))·Zv_ff(t); where Zv is the combined virtual impedance value; ω_ms is the multi - scale hybrid weight. Calculate the change rate of the Lyapunov function under the fused impedance: dV / dt|Zv = Icirc_comp T ·L·(dIcirc_comp / dt)|Zv; Verify the stability condition dV / dt|Zv < 0, if not satisfied, perform weight fine - tuning: ω_ms_adj(t) = ω_ms(t)-η·sign(dV / dt|Zv)·|dV / dt|Zv|, where η is the fine - tuning step size; Output the final combined virtual impedance value Zv(t) to ensure the stability and dynamic performance of the system.

[0054] According to another aspect of the present application, calculating the optimal weights for high-frequency perturbations, mid-frequency load variations, and low-frequency parameter drifts can also be: using the state transition data ST to optimize the multi-time scale weight allocation: reading the state transition data ST and analyzing the transition characteristics of different frequency bands; dynamically adjusting the frequency band weights according to the type of transition point: if ST shows that a high-frequency transition is about to occur, increase the weight of the high-frequency layer: λ_hf(t) = λ_hf(t) + Δλ_hf; if ST shows that a mid-frequency transition is about to occur, increase the weight of the mid-frequency layer: λ_mf(t) = λ_mf(t) + Δλ_mf; if ST shows that a low-frequency transition is about to occur, increase the weight of the low-frequency layer: λ_lf(t) = λ_lf(t) + Δλ_lf; renormalize the frequency band weights: λ_hf(t)+λ_mf(t)+λ_lf(t) = 1; output the optimized multi-scale hybrid weights considering state transitions.

[0055] In another embodiment of the present application, the output control data can also be: constructing a Lyapunov energy function for stability assurance: V(Icirc_comp,ω) = (1 / 2)L·Icirc_comp 2 + (1 / 2)kω·(ω-ω*) 2; where ω is the feedforward-feedback hybrid weight, ω* is the ideal weight, and kω is the weight adjustment coefficient. A gain modulation law based on Lyapunov stability is constructed to calculate the hybrid weight ω(t): dω(t) / dt = -γ·sign(Icirc_comp·d(Zv_fb - Zv_ff) / dt) - β·ΨV / Ψω; where γ is the adjustment rate and β is the damping coefficient. The multi-time scale gain modulation strategy is applied to calculate the fast-layer weight ωf(t), the medium-speed layer weight ωm(t), and the slow-layer weight ωs(t) respectively: dωf(t) / dt = -γf·ff(Icirc_comp, Zv); dωm(t) / dt = -γm·fm(Icirc_comp, Zv); dωs(t) / dt = -γs·fs(Icirc_comp, Zv); The combined weight: ω(t) = αf·ωf(t) + αm·ωm(t) + αs·ωs(t); where γf, γm, and γs are the adjustment rates of each layer, and αf, αm, and αs are the weight coefficients of each layer. Calculate the final combined virtual impedance value Zv, and calculate the voltage correction value ΔVi of each module according to the combined virtual impedance value Zv: ΔVi(t) = -Zv(t)·Icirc_comp,i(t); Considering the power balance between modules, the voltage correction value is coordinated and optimized to obtain the coordinated voltage correction value ΔVi_coord: ΔVi_coord(t) = ΔVi(t) + fi_balance(Pi(t), P*(t), Tj,i(t), T*j(t)); where fi_balance is the power balance function, Pi(t) is the power of module i, P*(t) is the average power, Tj,i(t) is the junction temperature of module i, and T*j(t) is the average junction temperature. Generate the final module control voltage command Vi_cmd: Vi_cmd(t) = Vref(t) + ΔVi_coord(t), where Vref(t) is the reference voltage. According to the module control voltage command Vi_cmd, the final PWM control signal PWMi is generated through PWM modulation: PWMi(t) = modulatePWM(Vi_cmd(t), Vdc, Tsw, modulation_strategy); where modulatePWM is the PWM modulation function and modulation_strategy is the modulation strategy.

[0056] In a specific embodiment of the present application, based on a three-phase four-module high-voltage SVG system with a total capacity of 10 MVA, a DC-side voltage of 3 kV, and an AC-side phase voltage of 10 kV / 50 Hz. The system includes 4 series H-bridge modules per phase, and parameter differences between the modules lead to circulating current problems. The main parameters of the experimental system are as follows: Rated capacity: Sn = 10 MVA; DC bus voltage: Vdc = 3000 V ± 5%; AC phase voltage: Vac = 10 kV / 50 Hz; Number of modules: 4 series modules per phase, a total of 12 modules; Single-module inductor: Lm = 2.0 mH ± 10%; Single-module capacitor: Cm = 5000 μF ± 7%; Sampling frequency: fs = 20 kHz; Switching frequency: fsw = 1000 Hz. The specific steps are as follows: Step 1: Data acquisition and preprocessing.

[0057] 1.1. Collect the original data of the four modules in phase A at a sampling frequency of 20 kHz for 100 ms (a total of 2000 data points): Original current data Iraw: [23.5 A, 24.2 A, 23.9 A, 24.7 A] (instantaneous values of the four modules at t = 10 ms); Original voltage data Vraw: [753 V, 745 V, 750 V, 748 V] (instantaneous values of the four modules at t = 10 ms); Switching state data Sraw: [1, 0, 0, 1; 1, 0, 0, 1; 1, 0, 0, 1; 1, 0, 0, 1] (switching states of the four modules at t = 10 ms); DC bus voltage data Vdc: [3050 V, 3030 V, 3070 V, 3020 V] (DC-side voltages of the four modules); Temperature data Tj: [45.2 °C, 47.5 °C, 44.8 °C, 48.3 °C] (IGBT junction temperatures of the four modules).

[0058] 1.2. Switching state sensing filtering calculation If(t)=Iraw(t)-Σ[α_k·fsw(t-tk,Sraw(tk), Sraw(tk+1))]; where If(t) is the filtered current data; Iraw(t) is the original current data; tk is the switching time; fsw is the switching transient model function; α_k is the adaptive weight; t is the time. Taking the switching transient at t = 10.005 ms as an example, it is recognized that the switching state of module 1 changes from [1, 0, 0, 1] to [1, 0, 1, 0]: The calculated value of the switching transient model fsw is 2.7 A; The adaptive weight α_k is 0.85; The original current Iraw(10.005 ms) is the peak value of 30.2 A; The filtered current If(10.005 ms)=30.2 A - 0.85×2.7 A = 27.905 A.

[0059] 1.3. Calculation of circulating current: Icirc,i(t) = Ii(t) - (1 / N)·Σ[Ij(t)]; where Icirc,i(t) is the circulating current component of the i-th module; Ii(t) is the current of the i-th module; N is the total number of parallel modules; j is the module index; t is time. Calculated with data at t = 10 ms: Average current: (23.5 A + 24.2 A + 23.9 A + 24.7 A) / 4 = 24.075 A; Circulating current of module 1: Icirc,1 = 23.5 A - 24.075 A = -0.575 A; Circulating current of module 2: Icirc,2 = 24.2 A - 24.075 A = 0.125 A; Circulating current of module 3: Icirc,3 = 23.9 A - 24.075 A = -0.175 A; Circulating current of module 4: Icirc,4 = 24.7 A - 24.075 A = 0.625 A.

[0060] 1.4. Compensation loop current calculation: Icirc_comp,i(t) = Icirc,i(t)·(1 + ηi(Vgs,i, Tj,i, Vdc)); where Icirc_comp,i(t) is the compensated loop current; Icirc,i(t) is the original loop current component; ηi is the compensation coefficient; Vgs,i is the gate voltage; Tj,i is the junction temperature; Vdc is the DC voltage; t is the time; and i is the module index. Calculation of the parameter difference compensation coefficient: η1 = 0.03×(Vgs,1 - 15) / 15 + 0.05×(Tj,1 - 45) / 45 + 0.02×(Vdc,1 - 3000) / 3000 = 0.00095; η2 = 0.03×(15.1 - 15) / 15 + 0.05×(47.5 - 45) / 45 + 0.02×(3030 - 3000) / 3000 = 0.00518; η3 = 0.03×(14.9 - 15) / 15 + 0.05×(44.8 - 45) / 45 + 0.02×(3070 - 3000) / 3000 = 0.00005; η4 = 0.03×(15.0 - 15) / 15 + 0.05×(48.3 - 45) / 45 + 0.02×(3020 - 3000) / 3000 = 0.0038; Compensated loop current: Icirc_comp,1 = -0.575A×(1 + 0.00095) = -0.5756A; Icirc_comp,2 = 0.125A×(1 + 0.00518) = 0.1256A; Icirc_comp,3 = -0.175A×(1 + 0.00005) = -0.1750A; Icirc_comp,4 = 0.625A×(1 + 0.0038) = 0.6274A.

[0061] 1.5. The preprocessed voltage data Vp at t = 10 ms: [750.2V, 742.8V, 748.1V, 745.3V]; The preprocessed status data Sp at t = 10 ms: {Module 1: [1, 0, 0, 1], Module 2: [1, 0, 0, 1], Module 3: [1, 0, 0, 1], Module 4: [1, 0, 0, 1], Timestamp: 10 ms}.

[0062] Step 2. Topology-aware phase space reconstruction.

[0063] 2.1. Divide the compensated circulating current Icirc_comp within the 20ms window according to the switch states. The following typical switch states are identified: State 1 [1, 0, 0, 1]: t ∈ [0 - 5ms, 10 - 15ms]; State 2 [1, 0, 1, 0]: t ∈ [5 - 10ms, 15 - 20ms]. Example values of the segmented circulating current data Icirc_seg: State 1 (t = 10ms): [-0.5756A, 0.1256A, -0.1750A, 0.6274A]; State 2 (t = 8ms): [-0.4923A, 0.1302A, -0.1829A, 0.5450A].

[0064] 2.2. Extract the topological connection matrix Tconn according to the switch states: Topological matrix for State 1 [1, 0, 0, 1]: [[1, 0], [0, 1], [0, 1], [1, 0]]; Topological matrix for State 2 [1, 0, 1, 0]: [[1, 0], [0, 1], [1, 0], [0, 1]]. Calculate the equivalent circuit parameters: Equivalent inductance for State 1: Leq(S1) = 1.95mH; Equivalent resistance for State 1: Req(S1) = 0.12Ω; Equivalent inductance for State 2: Leq(S2) = 2.05mH; Equivalent resistance for State 2: Req(S2) = 0.13Ω. Analyze the frequency characteristics of the circulating current: Main frequency for State 1: f_main(S1) = 270Hz, bandwidth BW(S1) = 150Hz; Main frequency for State 2: f_main(S2) = 250Hz, bandwidth BW(S2) = 140Hz. Calculate the optimal time delay: Delay factors: α = 0.25, β = 0.5, reference inductance Leq_ref = 2.0mH; Time delay for State 1: τ1 = 0.25×(1 / 270)×(1.95 / 2.0) 0.5 ≈ 0.9ms; Time delay for State 2: τ2 = 0.25×(1 / 250)×(2.05 / 2.0) 0.5≈ 1.0 ms. Evaluate the signal complexity: Complexity index for state 1: CI(S1) = 0.65; Complexity index for state 2: CI(S2) = 0.71. Calculate the optimal embedding dimension: Basic dimension: m_base = N_active + 1 = 4 + 1 = 5; Complexity compensation (γ = 0.5): m_comp(S1) = ⌈0.65 × 0.5⌉ = 1; Complexity compensation (γ = 0.5): m_comp(S2) = ⌈0.71 × 0.5⌉ = 1; Embedding dimension for state 1: m1 = min(5 + 1, 8) = 6; Embedding dimension for state 2: m2 = min(5 + 1, 8) = 6. Form the reconstruction parameter set RP: State 1: {τ1 = 0.9 ms, m1 = 6, Tconn1 = [[1, 0], [0, 1], [0, 1], [1, 0]]}; State 2: {τ2 = 1.0 ms, m2 = 6, Tconn2 = [[1, 0], [0, 1], [1, 0], [0, 1]]}. Calculate the parameter sensitivity index PS: State 1: PS(S1) = 0.72; State 2: PS(S2) = 0.85.

[0065] 2.3 Variable - length adaptive sliding window: Main period of the circulating current: T_circ = 1 / 260 Hz ≈ 3.85 ms; Safety factor: β = 2.5; Minimum window length: Lmin = 5 ms; Window length for state 1: Lwindow(S1) = max(5 ms, 2.5 × 3.85 ms) = 9.63 ms; Window length for state 2: Lwindow(S2) = max(5 ms, 2.5 × 4.0 ms) = 10.0 ms; Extraction of key feature points (taking state 1 as an example): Adaptive threshold: Δ(t) = 0.1 × |Icirc_comp|_max = 0.1 × 0.63 A = 0.063 A; Rate - of - change threshold: vth(t) = 0.15 × |dIcirc_comp / dt|_max = 0.15 × 250 A / s = 37.5 A / s; 23 key feature points were identified in a 10 - ms window, accounting for approximately 11.5% of the total number of points.

[0066] 2.4. Base Delay Embedding Vector: Taking t = 10 ms, State 1, and Module 1 as examples: τ1 = 0.9 ms, m1 = 6; Xbase(10 ms)=[Icirc_comp,1(10 ms),Icirc_comp,1(9.1 ms),Icirc_comp,1 (8.2 ms),Icirc_comp,1(7.3 ms),Icirc_comp,1(6.4 ms),Icirc_comp,1(5.5 ms)]; Xbase(10 ms)=[-0.5756 A, -0.5623 A, -0.5342 A, -0.4981 A, -0.4625 A, -0.4283 A]. Extract voltage features: Output voltage of Module 1: Vout,1(10 ms)=750.2 V; Voltage change rate: dV / dt(10 ms)=(750.2 V - 745.8 V) / 0.1 ms = 44 kV / s; Voltage unbalance: Vunbal(10 ms)=0.0045. Calculate the dynamic characteristics of the circulating current: Change rate: dI / dt(10 ms)=(-0.5756 A + 0.5623 A) / 0.9 ms = -14.8 A / s; Acceleration: d 2 I / dt 2 (10 ms)= -252 A / s 2 ; Change trend: Itrend(10 ms)=-1 (decreasing trend). Calculate temperature characteristics: Junction temperature difference: ΔTj(10 ms)=48.3°C - 44.8°C = 3.5°C; Temperature change rate: dTj / dt(10 ms)=0.042°C / s. Construct the weight fusion matrix (steady-state operating condition): W_fusion = [0.35, 0.15, 0.15, 0.15, 0.10, 0.10] (corresponding to the delay sequence, voltage, current dynamics, temperature difference, temperature change rate, switch characteristics); Generate the enhanced state space vector (Module 1, t = 10 ms): Delay sequence component: 0.35×[-0.5756 A, -0.5623 A, -0.5342 A, -0.4981 A, -0.4625 A, -0.4283 A]; Voltage feature component: 0.15×[750.2 V, 44 kV / s, 0.0045]; Current dynamic component: 0.15×[-14.8 A / s, -252 A / s 2 , -1]; Temperature characteristic component: 0.10×[3.5°C, 0.042°C / s]; Switch characteristic component: 0.10×[0.12 μs, 0.25 μs, 2.1 V]; The complete enhanced state vector X(10 ms) contains all the above components. Normalization processing gives the standardized enhanced state vector Xnorm(10 ms).

[0067] 2.5. Calculation of geometric features of the phase trajectory: Curvature calculation (t = 10 ms): κ(10 ms) = 25.7; Trajectory twist: τ(10 ms) = 158.3; Geometric feature vector: GF(10 ms)=[25.7, 158.3]. Local stability assessment: Trajectory separation of the tracking points (10 ms, 11 ms, 12 ms); Calculation of the local Lyapunov exponent: λL(10 ms) = -0.032; Stability index: SI(10 ms)=0.87. Calculation of the phase space velocity and acceleration: Velocity vector: v(10 ms)=[-14.8 A / s, 178.5 V / s, -252 A / s 2 , 0.042 °C / s]; Acceleration vector: a(10 ms)= [-128 A / s 2 , 2350 V / s 2 , -3680 A / s 3 , 0.009 °C / s 2 ; Dynamic evolution feature: DF(10 ms)=[||v(10 ms)||, ||a(10 ms)||, v(10 ms)·a(10 ms) / ||v(10 ms)||]=[294.8, 4327.6, -1215.4]. Periodicity and recurrence analysis: Recurrence plot index: RQA = 0.73; Periodicity intensity: PI(10 ms)=0.68; Cycle pattern: CP(10 ms)=2 (corresponding to pattern 2 in the standard pattern library); Period feature: PF(10 ms)=[0.68, 2, 0.73]. Calculation of the energy flow feature: Circulating current energy: Ecirc(10 ms)=∫ 10ms 10·9ms 1.95 mH·Icirc_comp 2 (τ)dτ / 2 = 3.25 μJ; Power flow: Pcirc(10 ms) = dEcirc(10 ms) / dt=-0.28 W; Energy flow direction index (γ = 0.75): Denergy(10 ms)=sign(-0.28)·|-0.28| 0.75= -0.3853; Energy flow characteristic: EF(10ms)=[3.25μJ, -0.28W, -0.3853]. Parameter sensitivity characteristic extraction: Resistance difference sensitivity: SR(10ms) = 2.35A / Ω; Inductance difference sensitivity: SL(10ms)=0.87A / mH; Capacitance difference sensitivity: SC(10ms) =0.0023A / μF; Parameter sensitivity characteristic: SF(10ms)=[2.35A / Ω, 0.87A / mH, 0.0023A / μF]. Feature weight determination and formation of circulating current phase characteristics: Trajectory geometric feature weight: 0.25; Stability index weight: 0.20; Dynamic evolution feature weight: 0.20; Energy flow feature weight: 0.15; Parameter sensitivity feature weight: 0.20; The final circulating current phase characteristic CFPS(10ms) integrates the above features and their weights.

[0068] 2.6. Calculation of the coupling strength matrix between modules: C(10ms)=[[1.00, 0.65, -0.83, -0.75], [0.65, 1.00, -0.58, -0.92], [-0.83, -0.58, 1.00, 0.71], [-0.75, -0.92, 0.71, 1.00]]; Dynamic trend prediction (taking Module 1 as an example): Construction of a local linear prediction model: Select 10 historical points similar to the current state to construct the mapping Flocal; Matrix coefficients: A = [0.92, 0.05, -0.03, 0.02, -0.01, 0.00], b = -0.0125. Generation of a multi-step prediction sequence: Prediction at t+1ms: Icirc_pred,1(11ms) = -0.5832A; Prediction at t+2ms: Icirc_pred,1(12ms) = -0.5901A; Prediction at t+5ms: Icirc_pred,1(15ms) = -0.6023A. Extraction of circulating current prediction data: Predicted peak value: Ipeak,1 = -0.6053A (t=17ms); Predicted root mean square value: Irms,1 = 0.5897A; Prediction credibility: Conf,1 = 0.92. Final dynamic trend prediction data: Prediction sequence: DTP = [Icirc_pred(11ms), Icirc_pred(12ms),..., Icirc_pred(20ms)]; Trend index: max(dI / dt) = -15.3A / s, ΔE = -0.12μJ, Stability = 0.85.

[0069] Step 3. Calculation of virtual impedance parameters for working condition adaptation.

[0070] 3.1. Calculation of system operating condition index vector: Load change rate: |dIload / dt| = 128 A / s; Grid fluctuation: |dVgrid / dt| = 85 V / s; Active power: Pactive = 7.2 MW; Reactive power: Qreactive = 5.8 MVAr; Harmonic distortion: THD = 0.027; Module coupling degree: C_norm = 0.78; Operating condition index vector: W(10ms) = [128 A / s, 85 V / s, 7.2 MW, 5.8 MVAr, 0.027, 0.78]. Operating condition conversion detection: Operating condition index at the previous moment: W(9ms) = [125 A / s, 83 V / s, 7.15 MW, 5.75 MVAr, 0.026, 0.77]; Change rate: ||W(10ms) - W(9ms)|| / ||W(9ms)|| = 0.015; Threshold θW = 0.05; 0.015 < 0.05, no operating condition conversion.

[0071] 3.2. Calculation of energy flow characteristic data: Circulating current energy: Ecirc(10ms) = 3.25 μJ; Power flow: Pcirc(10ms) = -0.28 W; Energy flow direction index: Denergy(10ms) = -0.3853; Energy flow characteristic data: EF(10ms) = [3.25 μJ, -0.28 W, -0.3853]. Parameter sensitivity data: Resistance difference sensitivity: SR(10ms) = 2.35 A / Ω; Inductance difference sensitivity: SL(10ms) = 0.87 A / mH; Capacitance difference sensitivity: SC(10ms) = 0.0023 A / μF; Parameter sensitivity data: S(10ms) = [2.35 A / Ω, 0.87 A / mH, 0.0023 A / μF]. Complexity characteristics of circulating current signal: Sample entropy calculation (m = 2, r = 0.15 × standard deviation): Proportion of similar subsequences: A = 0.172, B = 0.258; Sample entropy value: SampEn(10ms) = -ln(0.172 / 0.258) = -ln(0.667) = 0.405. Multiscale entropy analysis (τ = 1, 2, 3): MSE(τ = 1) = 0.405, MSE(τ = 2) = 0.387, MSE(τ = 3) = 0.352; Multiscale entropy characteristics: MSE = [0.405, 0.387, 0.352]. Permutation entropy calculation (pattern length m = 3): Probability distribution of different permutation patterns: p(π) = [0.21, 0.18, 0.15, 0.12, 0.11, 0.10, 0.08, 0.05]; Permutation entropy: PermEn(10ms) = -Σ p(π)·log(p(π)) = 1.982. Fuzzy entropy calculation: Average value of fuzzy similarity matrix: A = 0.185, B = 0.264; Fuzzy entropy value: FuzzyEn(10ms) = -ln(0.185 / 0.264) = 0.356. Entropy change rate calculation: Sample entropy at the previous moment: SampEn(9ms) = 0.398; Sample entropy change rate: dSampEn(10ms) = (0.405 - 0.398) / 1ms = 0.007 / ms. Calculation of comprehensive complexity index: Entropy index fusion weight: W_entropy = [0.30, 0.25, 0.20, 0.15, 0.10]; Signal complexity characteristics: CI(10ms) = 0.30×0.405 + 0.25×0.387 + 0.20×1.982 + 0.15×0.356 + 0.10×0.007 = 0.6688.

[0072] 3.3. Determination of Key Control Moments: Predicted peak moment \(t_{peak}=17\mathrm{ms}\); Predicted zero crossing \(t_{zero}=25\mathrm{ms}\) (next cycle); Moment of maximum change rate \(t_{max\_rate}=13\mathrm{ms}\); Moment of maximum energy \(t_{max\_energy}=16\mathrm{ms}\); Weight of key moments \(W_{time}=[0.35, 0.20, 0.25, 0.20]\); Key moment data \(KT = \{t_{peak}:17\mathrm{ms},t_{zero}:25\mathrm{ms},t_{max\_rate}:13\mathrm{ms},t_{max\_energy}:16\mathrm{ms}\}\). Predicted energy characteristics: Future energy flow direction \(Denergy\_pred(15\mathrm{ms})=-0.412\); Predicted energy peak \(Epeak\_pred = 3.42\mu J(t = 16\mathrm{ms})\); Energy change sensitivity \(E\_sensitivity = 2.35\); Predicted energy characteristics \(PE = [Denergy\_pred,Epeak\_pred,E\_sensitivity]\). Predicted parameter sensitivity: Predicted resistance sensitivity \(SR\_pred(15\mathrm{ms}) = 2.42\mathrm{A} / \Omega\); Predicted inductance sensitivity \(SL\_pred(15\mathrm{ms}) = 0.91\mathrm{A} / \mathrm{mH}\); Predicted capacitance sensitivity \(SC\_pred(15\mathrm{ms}) = 0.0025\mathrm{A} / \mu F\); Predicted sensitivity values \(PS=[SR\_pred,SL\_pred,SC\_pred]\). Predicted complexity characteristics: Predicted entropy value change \(CI\_pred(15\mathrm{ms}) = 0.683\); Predicted state transition point \(T\_change\_pred = 22\mathrm{ms}\); Predicted complexity characteristics \(PC=[CI\_pred,T\_change\_pred]\). Impedance - circulating current mapping model: Establish the mapping relationship \(Z\_I\_model\) based on historical data; For state 1, the relationship between the circulating current peak and the virtual impedance is: \(Ipeak = 0.98\times Ipeak\_0\times\exp(-0.27\times Z)\); For module 1, the prediction relationship at \(t = 15\mathrm{ms}\): \(Icirc\_pred,1(15\mathrm{ms})=-0.6023\mathrm{A}\times\exp(-0.27\times Z)\). Design of multi - objective optimization function: Circulating current suppression objective: \(J_{circ}=\sum|Icirc\_pred,i|\) 2 = 0.6023 2 + 0.1323 2 + 0.1853 2 + 0.6353 2 =0.7932; Dynamic response objective: \(J_{resp}=\sum|dIcirc\_pred,i / dt|\) 2 = 15.3 2 + 6.2 2 + 7.5 2+ 16.1 2 =537.59; Working condition adaptive weight: α=0.65 (dynamically calculated according to working condition index W); Comprehensive objective function: J = 0.65×0.7932 + (1-0.65)×537.59 / 1000= 0.7038. Optimal virtual impedance solution (taking t=15ms as an example): Objective function: J(Z)=0.65×(Σ|Icirc_pred,i(15ms)| 2 )+0.35×(Σ|dIcirc_pred,i(15ms) / dt| 2 ) / 1000; Constraints: 0.05Ω≤Z≤1.0Ω, |dZ / dt|≤5Ω / s; Optimization result: Z_optimal(15ms) =0.52Ω. Feedforward virtual impedance parameter sequence (future 10ms window): Zv_ff(11ms) = 0.46Ω; Zv_ff(12ms) = 0.48Ω; Zv_ff(13ms) = 0.50Ω; Zv_ff(14ms) = 0.51Ω; Zv_ff(15ms) = 0.52Ω; Zv_ff(16ms) = 0.54Ω; Zv_ff(17ms) = 0.56Ω; Zv_ff(18ms) = 0.55Ω; Zv_ff(19ms) = 0.53Ω; Zv_ff(20ms) = 0.50Ω.

[0073] 3.4. Key features extraction of the circulating current phase characteristic CFPS: Trajectory geometric feature: GF(10ms) = [25.7, 158.3]; Stability index: SI(10ms) = 0.87; Dynamic evolution feature: DF(10ms) = [294.8, 4327.6, -1215.4]. Adaptive PID controller design: Adjust KP based on the trajectory curvature: KP = 0.5×(1 + 0.2×tanh(25.7 - 20)) = 0.6; Adjust KI based on the stability index: KI = 0.1×(1 + 0.5×(0.87 - 0.5)) = 0.1185; Adjust KD based on the velocity vector: KD = 0.05×(1 + 0.3×tanh(294.8 - 200)) = 0.065; Adaptive PID parameters: KPID = [0.6, 0.1185, 0.065]. Fine-tuning of parameter sensitivity features and periodic features: Inductance sensitivity threshold: SL_threshold = 0.8A / mH; SL(10ms) = 0.87A / mH > 0.8A / mH, adjustment is required; KP_adj = KP×(1 + 0.15×(SL - SL_threshold) / SL_threshold) = 0.6×(1 + 0.15×(0.87 - 0.8) / 0.8) = 0.6079; Periodic intensity adjustment: KI_adj = KI×(1 - 0.2×(PI - 0.5)) = 0.1185×(1 - 0.2×(0.68 - 0.5)) = 0.1142; Preset adjustment based on the circulating pattern CP = 2: KD_adj = KD×1.05 = 0.065×1.05 = 0.0683; Feature matching parameters: FMP = [1.0131, 0.964, 1.05]. Calculation of the base impedance value: Circulating current of module 1 at t = 10ms: Icirc_comp,1(10ms) = -0.5756A; Integral term (assuming a 10ms window in the past): ∫Icirc_comp,1(τ)dτ = -5.2mA·s; Differential term: dIcirc_comp,1(10ms) / dt = -14.8A / s; Basic PID control: Zv_base,1 = 0.6079×(-0.5756) + 0.1142×(-5.2×10 -3 ) + 0.0683×(-14.8) = -1.3613Ω; Since the virtual impedance is positive, take the absolute value: Zv_base,1 = 1.3613Ω. State feedback enhancer design: Extract the subset of state variables: Xsub = [Icirc_comp, dIcirc_comp / dt, d 2 Icirc_comp / dt2 ; State feedback gain matrix: K_state = [0.25, 0.02, 0.001]; State feedback component: Zv_state,1 = 0.25×(-0.5756) + 0.02×(-14.8) + 0.001×(-252) = -0.6919 Ω; Take the absolute value: Zv_state,1 = 0.6919 Ω. Resonant suppression compensator design: Identify the resonant frequency: f_res = 270 Hz; Resonant suppressor parameters: A_res = 0.18 Ω, φ_res = 180°; Suppression component at t = 10 ms: Zv_res,1 = 0.18×sin(2π×270×0.01 + π) = 0.1712 Ω. Stability enhancer design: Stability threshold: λL_threshold = 0.05; λL(10 ms) = -0.032 < λL_threshold, no enhancement required; Stability enhancement component: Zv_stab,1 = 0 Ω. Feedback virtual impedance parameter integration: Component weights (dynamically adjusted based on CFPS): W_fb = [0.45, 0.30, 0.25, 0]; Feedback virtual impedance calculation: Zv_fb,1 = 0.45×1.3613 + 0.30×0.6919 + 0.25×0.1712 + 0×0 = 0.863 Ω.

[0074] 3.5 Physical Limitations of Power Electronic Devices: Maximum allowable impedance value: Zmax(10ms) = 2.0Ω; DC bus voltage: Vdc = 3000V; Rated current: Irated = 100A; Switching period: Tsw = 1ms; Dynamic safety factor: η(10ms) = 0.5; Maximum rate of change limitation: |dZv / dt|_max = 0.5×3000 / (100×0.001) = 15Ω / s. The feedforward virtual impedance parameter after constraint: |Zv_ff(15ms)| = 0.52Ω < Zmax(10ms) = 2.0Ω, satisfying the constraint; |dZv_ff / dt|_max = |(0.52 - 0.51) / 0.001| = 10Ω / s < 15Ω / s, satisfying the constraint; Constrained virtual impedance parameter: Zv_ff_constrained = Zv_ff = 0.52Ω (t = 15ms). The feedback virtual impedance parameter after constraint: |Zv_fb,1| = 0.863Ω < Zmax(10ms) = 2.0Ω, satisfying the constraint; |dZv_fb / dt|_max = |(0.863 - 0.845) / 0.001| = 18Ω / s > 15Ω / s, exceeding the constraint; Constraint handling: Zv_fb_constrained,1 = 0.845 + 15×0.001 = 0.86Ω.

[0075] Step 4: Hybrid Control and Multi-module Coordination for Stability Assurance.

[0076] 4.1 System State Vector Definition: Circulating current vector: Icirc_comp(10ms) = [-0.5756A, 0.1256A, -0.1750A, 0.6274A] T ; Initial hybrid weight vector: ω(10ms) = [0.7, 0.7, 0.7, 0.7] T ; Ideal weight vector: ω* = [0.65, 0.65, 0.65, 0.65] T . Lyapunov function parameter settings: Equivalent inductance matrix: L = diag(1.95mH, 1.95mH, 1.95mH, 1.95mH); Weight adjustment coefficient matrix: kω = diag(0.5, 0.5, 0.5, 0.5); Positive definite weight matrix: Q = [[1, 0.2, 0.2, 0.2], [0.2, 1, 0.2, 0.2], [0.2, 0.2, 1, 0.2], [0.2, 0.2, 0.2, 1]]. Lyapunov function calculation: Circulating current energy term: (1 / 2)L·Icirc_compT ·Icirc_comp = (1 / 2)×1.95×10 -3 ×((-0.5756) 2 + 0.1256 2 +(-0.1750) 2 + 0.6274 2 )=7.385×10 -4 J; Weight bias term: (ω - ω*) T = [0.05, 0.05, 0.05, 0.05] T ; (ω - ω*) T ·Q·(ω - ω*) = [0.05, 0.05, 0.05, 0.05]·[[1, 0.2, 0.2, 0.2], [0.2, 1, 0.2, 0.2], [0.2, 0.2, 1, 0.2], [0.2, 0.2, 0.2, 1]]·[0.05, 0.05, 0.05, 0.05] T =0.016; Weight energy term: (1 / 2)kω·(ω - ω*) T ·Q·(ω - ω*) = (1 / 2)×0.5×0.016 = 0.004; Lyapunov function value: V(10ms) = 7.385×10 -4 + 0.004 = 0.00474 J. Lyapunov function time derivative: Estimated circulating current change rate: dIcirc_comp(10ms) / dt = [-14.8, 3.2, -4.5, 15.2] T A / s; Circulating current energy change rate: Icirc_comp T ·L·(dIcirc_comp / dt) = [-0.5756, 0.1256, -0.1750, 0.6274]·1.95×10 -3 ×[-14.8, 3.2, -4.5, 15.2] T =0.03753 W; Weight change rate (tentative): dω / dt = [-0.01, -0.01, -0.01, -0.01] T s -1 ; Weight energy change rate: kω·(ω - ω*) T·Q·(dω / dt) = 0.5×[0.05, 0.05, 0.05, 0.05]·[[1, 0.2, 0.2, 0.2],[0.2, 1, 0.2, 0.2], [0.2, 0.2, 1, 0.2], [0.2, 0.2, 0.2, 1]]·[-0.01, -0.01, -0.01, -0.01] T = -0.0016 W; The time derivative of the Lyapunov function: dV / dt = 0.03753 + (-0.0016) = 0.03593 W. Stability criterion: Stability index: SI = -dV / dt = -0.03593 < 0, the system is unstable; Stability evaluation data: SE = {V(10ms)=0.00474J, dV / dt=0.03593W, SI=-0.03593}.

[0077] 4.2. Read the stability evaluation data SE: Lyapunov function value: V(10ms) = 0.00474J; Time derivative: dV / dt = 0.03593W; Stability index: SI = -0.03593. Weight adjustment guided by Lyapunov stability: Module 1 parameters (at 10ms): Icirc_comp,1=-0.5756A; Zv_fb,1 = 0.86Ω; Zv_ff,1 = 0.46Ω (current moment); |Zv_fb,1-Zv_ff,1|=|0.86-0.46|=0.4Ω; d(Zv_fb,1-Zv_ff,1) / dt=d(0.4) / dt=0.05Ω / ms; sign(Icirc_comp,1·d(Zv_fb, 1-Zv_ff,1) / dt) = sign((-0.5756)×0.05)=-1; Learning rate parameter γ = 0.2; Exponential parameter α = 0.5; Damping coefficient β = 0.05. Weight adjustment calculation: dω1 / dt = -0.2×(-1)×|-0.5756|×|0.4| 0.5 - 0.05×(0.7-0.65); dω1 / dt = 0.2×0.5756×0.6325 -0.05×0.05; dω1 / dt = 0.07276 - 0.0025 = 0.07026 s -1 . Similarly calculate dω / dt for other modules, weight adjustment vector: dω / dt = [0.07026, -0.00582,0.01621,-0.07965] T s -1Output base mixing weights: Updated weights (Δt = 1ms): ω1(11ms) = ω1(10ms) + dω1 / dt × Δt = 0.70070; ω2(11ms) = 0.7 - 0.00582 × 0.001 = 0.69994; ω3(11ms) = 0.7 + 0.01621 × 0.001 = 0.70016; ω4(11ms) = 0.7 - 0.07965 × 0.001 = 0.69920; Base mixing weights: ω(11ms) = [0.70070, 0.69994, 0.70016, 0.69920] T 。

[0078] 4.3. Circulating current frequency division decomposition: High-frequency component (>500Hz): Icirc_hf,1(10ms) = -0.0836A; Medium-frequency component (50 - 500Hz): Icirc_mf,1(10ms) = -0.4632A; Low-frequency component (<50Hz): Icirc_lf,1(10ms) = -0.0288A. Hierarchical Lyapunov function construction: High-frequency layer: V_hf(10ms) = 6.812×10 -6 J; Medium-frequency layer: V_mf(10ms) = 4.421×10 -4 J; Low-frequency layer: V_lf(10ms) = 1.611×10 -6 J; Frequency division layer weight calculation: High-frequency learning rate: γ_hf = 0.5; Medium-frequency learning rate: γ_mf = 0.2; Low-frequency learning rate: γ_lf = 0.05; High-frequency weight change: dω_hf,1 / dt = 0.1536 s -1 ; Medium-frequency weight change: dω_mf,1 / dt = 0.0582 s -1 ; Low-frequency weight change: dω_lf,1 / dt = 0.0103 s -1; High-frequency optimal weight: ω_hf,1(11ms) = 0.7 + 0.1536×0.001 = 0.7154; Medium-frequency optimal weight: ω_mf,1(11ms) = 0.7 + 0.0582×0.001 = 0.7058; Low-frequency optimal weight: ω_lf,1(11ms) = 0.7 + 0.0103×0.001 = 0.7010. Dynamic frequency band weight calculation: Calculated based on the frequency band energy distribution: High-frequency energy ratio: E_hf / E_total = 0.152; Medium-frequency energy ratio: E_mf / E_total = 0.831; Low-frequency energy ratio: E_lf / E_total = 0.017. Frequency band weight: λ_hf(10ms) = 0.152; λ_mf(10ms) = 0.831; λ_lf(10ms) = 0.017. Multi-scale hybrid weight calculation: Multi-scale weight of module 1: ω_ms,1(11ms) = 0.152×0.7154 + 0.831×0.7058 + 0.017×0.7010 = 0.7071; Similarly calculate the multi-scale hybrid weights of other modules; Final multi-scale hybrid weight vector: ω_ms(11ms)=[0.7071, 0.6998, 0.7003, 0.6991] T . Global stability verification: Synthetic Lyapunov function: V_total = V_hf + V_mf + V_lf = 6.812×10 -6 +4.421×10 -4 +1.611×10 -6 =4.505×10 -4 J; Time derivative calculation: dV_total / dt = dV_hf / dt + dV_mf / dt + dV_lf / dt = -1.253×10 -4 -9.872×10 -4 -8.35×10 -6 =-1.1209×10 -3 W; Result: dV_total / dt < 0, satisfying the stability condition.

[0079] 4.4. Read feedforward and feedback virtual impedance parameters: The feedback virtual impedance of Module 1: Zv_fb,1(11ms) = 0.865Ω; the feedforward virtual impedance of Module 1: Zv_ff,1(11ms) = 0.46Ω; the multi-scale mixing weight of Module 1: ω_ms,1(11ms) = 0.7071. Calculate the comprehensive virtual impedance value of Module 1: Zv,1(11ms)=0.7071×0.865+(1 - 0.7071)×0.46; Zv,1(11ms)=0.6116 + 0.1347 = 0.7463Ω. Similarly calculate the comprehensive virtual impedance of other modules: Zv,2(11ms)=0.6998×0.73+(1 - 0.6998)×0.42 = 0.637Ω; Zv,3(11ms)=0.7003×0.78+(1 - 0.7003)×0.44 = 0.6781Ω; Zv,4(11ms) = 0.6991×0.83+(1 - 0.6991)×0.45 = 0.7157Ω. Verify the stability of the fused impedance: Calculate the change rate of the Lyapunov function: dV / dt|Zv = -1.237×10 -3 W < 0, satisfying the stability condition. The final comprehensive virtual impedance value: Zv(11ms) = [0.7463, 0.637, 0.6781, 0.7157] T Ω.

[0080] 4.5. Calculation of module voltage correction value: Module 1 voltage correction value: ΔV1(11ms) = -Zv,1(11ms)·Icirc_comp,1(11ms)=-0.7463×(-0.5803) = 0.4331V; Module 2 voltage correction value: ΔV2(11ms) =-0.637×0.1267 = -0.0807V; Module 3 voltage correction value: ΔV3(11ms) = -0.6781×(-0.1762) =0.1195V; Module 4 voltage correction value: ΔV4(11ms) = -0.7157×0.6298 = -0.4507V. Power balance coordination optimization: Module 1 power: P1=23.5A×750.2V = 17.63kW; Module 2 power: P2 = 24.2A×742.8V=17.98kW; Module 3 power: P3=23.9A×748.1V = 17.88kW; Module 4 power: P4= 24.7A×745.3V =18.41kW; Average power: P*=17.975kW; Module 1 balance compensation: ΔV1 = 0.05×(P*-P1) / P* = 0.05×(17.975-17.63) / 17.975= 0.00096V; Similar calculations are made for the balance compensation of other modules; Coordinated voltage correction value: ΔV1_coord(11ms)=0.4331+ 0.00096 = 0.4341V; ΔV2_coord(11ms) = -0.0807 -0.00014=-0.0808V; ΔV3_coord(11ms)= 0.1195 + 0.00026 = 0.1198V; ΔV4_coord(11ms) = -0.4507 - 0.00121 = -0.4519V. Final control voltage command: Reference voltage (t=11ms): Vref = [751.3V,751.3V, 751.3V, 751.3V] T; Module control voltage: V1_cmd(11ms) = 751.3 + 0.4341 = 751.7341V; V2_cmd(11ms) = 751.3 - 0.0808 = 751.2192V; V3_cmd(11ms) = 751.3 + 0.1198 = 751.4198V; V4_cmd(11ms) = 751.3 - 0.4519 = 750.8481V. PWM control signal generation: PWM modulation strategy: Phase-shifted carrier SPWM; Carrier frequency: 1000Hz; DC bus voltage: Vdc = [3050V, 3030V, 3070V, 3020V]; Modulation ratio of module 1: m1 = V1_cmd / Vdc,1 = 751.7341 / 3050 = 0.2465; Calculate the duty cycle and phase of the PWM wave of each module to generate the final PWM control signal.

[0081] To verify the effectiveness of this embodiment, a comparative test is carried out with the traditional control method. The traditional methods include: Fixed virtual impedance method (FVI): Use a fixed virtual impedance value for control; Single PID feedback method (SPID): Only use PID-based feedback control; Conventional feedforward-feedback hybrid method (CFFB): Use fixed weights to fuse feedforward and feedback control. The test is carried out under the same hardware platform and working conditions, and the following indicators are analyzed by comparison; The peak value of the circulating current in this embodiment is reduced to 12.5% of the original, and the root mean square value is reduced by 87.3%; The peak value of the circulating current of the FVI method is reduced to 43.2% of the original, and the root mean square value is reduced by 59.7%; The peak value of the circulating current of the SPID method is reduced to 35.8% of the original, and the root mean square value is reduced by 64.2%; The peak value of the circulating current of the CFFB method is reduced to 21.3% of the original, and the root mean square value is reduced by 79.5%. The dynamic response performance (stable time after step load change) of this embodiment is 4.2ms, that of the FVI method is 15.7ms, that of the SPID method is 9.3ms, and that of the CFFB method is 7.8ms. The adaptability of this embodiment to grid disturbances (circulating current fluctuation during 15% voltage sag) is 7.3%, that of the FVI method is 28.5%, that of the SPID method is 19.2%, and that of the CFFB method is 14.7%. The parameter robustness of this embodiment (system stability under the condition of ±15% deviation of module parameters) is all stable, and the circulating current increases by 10.2%; The FVI method is unstable when the parameter deviation is more than ±12%; When the parameter deviation of the SPID method is ±15%, the circulating current increases by 42.5%; When the parameter deviation of the CFFB method is ±15%, the circulating current increases by 27.8%.

[0082] The present invention realizes the precise matching between the phase space reconstruction parameters and the real-time operating state of the power electronic system, enabling the phase space to capture the dynamic characteristics of the circulating current more accurately, especially in scenarios where the switching state changes rapidly. For high-voltage SVG systems with high-frequency switching, it can adaptively cope with the system dynamics changes under different switching combinations, fundamentally improving the effectiveness of subsequent control strategies. Especially in high-power modular systems, the circulating current suppression efficiency is increased by more than 3 times. By constructing a multi-physical enhanced state space that integrates electrical-thermal-switching characteristics, the limitations of the traditional single-current state space are fundamentally overcome, comprehensively capturing the coupling relationships between different physical quantities in the system. During the actual operation of high-voltage SVG, it can simultaneously reflect the influence of the junction temperature difference of power devices on the circulating current and the system asymmetry caused by voltage imbalance, enabling the control system to anticipate and respond to the changes in the circulating current under the combined action of multiple physical factors. Especially in high-voltage high-power scenarios, the system stability margin is improved. By calculating the geometric characteristics of the phase trajectory, a complete set of circulating current phase characteristics CFPS is established, directly reflecting the inherent characteristics of the circulating current dynamics system. Different from traditional methods that only focus on the surface characteristics of amplitude and spectrum, the present invention deeply analyzes the geometric shape and stability of the phase space trajectory, capable of identifying the critically stable states and nonlinear dynamic behaviors in the circulating current. In practical applications, it enables the system to identify precursor signals before the circulating current becomes unstable. Especially in scenarios of sudden load changes, compared with traditional methods, it can initiate preventive control earlier, suppressing the circulating current peak below the safety threshold, improving the reliability and robustness of the system. By simultaneously calculating multi-dimensional entropy indicators, a comprehensive signal complexity analysis framework is established, achieving the precise quantification of the non-stationary dynamic characteristics of the circulating current. Beyond traditional statistical feature analysis, it can deeply reveal the inherent complexity and predictability of the circulating current from the perspective of information theory. Especially under conditions of grid harmonics and rapid load fluctuations, the complexity characteristics based on entropy analysis can accurately identify the system behavior patterns, predict the drastic changes in the circulating current in advance, providing sufficient response time for the control system; directly improving the pertinence of virtual impedance parameter optimization and enhancing the safety margin of system operation. By constructing a stability-oriented gain modulation law, stable control with strict mathematical proof is achieved, theoretically ensuring the stable convergence of the system and solving the fundamental problem that traditional empirical tuning methods lack stability guarantee. During the actual operation of high-voltage SVG systems, it ensures that the derivative of the Lyapunov function of the system is always negative under various operating conditions, thus theoretically guaranteeing the stability of the system. It is especially suitable for grid fluctuation scenarios in high-voltage transmission systems, improving the reliability and robustness of SVG systems. By decomposing the circulating current into high-frequency, medium-frequency, and low-frequency components, designing independent control layers for different frequency bands, and dynamically adjusting the control weights based on the frequency band energy distribution, precise control for disturbances at different time scales is achieved.It overcomes the limitation that it is difficult for a traditional single controller to handle fast switching transients and slow parameter drifts simultaneously, enabling the control system to specifically process circulating current components of different frequencies. In practical applications, the high-frequency control layer can quickly suppress the circulating current spikes caused by switching transients, the medium-frequency control layer processes the fluctuations caused by load changes, and the low-frequency control layer compensates for long-term drifts such as device aging. The test results show that, compared with single-band control, the present invention improves the efficiency while maintaining effective compensation for low-frequency drifts; it can maintain the best performance under various disturbances on different time scales, comprehensively improving the dynamic response ability and steady-state accuracy of the SVG system. The present invention effectively solves the problem of circulating current caused by parameter differences between modules, improving the reliability, stability, and power quality of system operation.

[0083] 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. Modular high-voltage SVG circulating current suppression method, characterized in that, Including: Obtain module data, perform preprocessing to obtain preprocessed data, including circulating current data, voltage data, and status data; Construct an enhanced phase space based on the preprocessed data to generate feature data, including circulating current phase features and dynamic trend prediction data; Based on the feature data, combined with system physical constraints, calculate virtual impedance parameters, including feedforward and feedback virtual impedance parameters; Use the virtual impedance parameters for gain modulation to generate a control signal and output control data, including the module control voltage correction value and the PWM control signal.

2. The method according to claim 1, characterized in that, Generate feature data, including: Analyze the status data, split the circulating current data into a predetermined number of subsequences according to the switch status to form segmented circulating current data; Read the system topology structure and switch status, dynamically calculate the time delay and embedding dimension for the subsequences to generate reconstruction parameters; Extract key feature points from the circulating current data, and combined with the reconstruction parameters, construct an enhanced state space, fuse the time-delay sequence in the circulating current data with physical quantities to form an enhanced state vector, and calculate the geometric characteristics of the phase trajectory based on this to generate circulating current phase features; Evaluate the coupling relationship of the circulating current between modules to form a module coupling matrix, and combined with the circulating current phase features, generate dynamic trend prediction data.

3. The method according to claim 2, characterized in that, Generate reconstruction parameters, including: Extract the topological connection matrix at the current switching moment from the status data, and calculate the equivalent circuit parameters under the current topology based on this; Analyze the frequency characteristics of the segmented circulating current data, determine the main frequency component and bandwidth, and calculate the optimal time delay parameter combined with the equivalent circuit parameters; Evaluate the complexity of the segmented circulating current data, combined with the total number of system modules and the number of currently active modules, determine the optimal embedding dimension; Associate the optimal time delay parameter and embedding dimension with the topological connection matrix to form reconstruction parameters.

4. The method according to claim 2, wherein Form an enhanced state vector, including: Construct a basic time-delay embedding vector according to the circulating current data and reconstruction parameters; Extract key voltage features from the voltage data; Calculate the dynamic characteristics of the circulating current data; Based on the pre-stored system temperature data, calculate the temperature characteristics of key power devices, including the junction temperature difference and the temperature change rate; Weightedly fuse the basic time-delay embedding vector with the key voltage features, dynamic characteristics, and temperature characteristics through a multi-physical quantity weight fusion matrix to form an enhanced state vector.

5. The method according to claim 2, wherein Generate circulating current phase features, including: Read and based on the enhanced state vector, in sequence: calculate the trajectory curvature and twist to form trajectory geometric features; evaluate the local stability of the phase trajectory to form a stability index; calculate the phase space velocity vector and acceleration vector to form dynamic evolution features; Based on the circulating current data and status data, form an energy flow feature; Extract the sensitivity of the circulating current data to module parameter differences to form a parameter sensitivity feature; Determine the weights of the trajectory geometric features, stability index, dynamic evolution features, energy flow features, and parameter sensitivity features to form circulating current phase features.

6. The method according to claim 2, wherein Generate dynamic trend prediction data, including: Based on the enhanced state vector and circulating current phase features, select a set of points similar to the current state from historical data to construct a local linear prediction model; Based on the local linear prediction model, the current state is used to predict the state at the next moment, and then the prediction result is used as the new initial condition to predict the more distant future, generating a multi-step prediction sequence; Extract the predicted circulating current data from the multi-step prediction sequence, calculate the key dynamic indicators of the predicted circulating current, and integrate them into the dynamic trend prediction data.

7. The method according to claim 1, wherein Calculate the virtual impedance parameters, including: Analyze the energy distribution and flow characteristics of the circulating current data within the switching period to generate energy flow characteristic data; Evaluate the sensitivity of the circulating current data to the differences in module parameters to form parameter sensitivity data; Calculate the sample entropy and entropy change rate of the enhanced phase space, and analyze the complexity characteristics of the circulating current; combine it with the dynamic trend prediction data, energy flow characteristic data, and parameter sensitivity data to calculate the optimal virtual impedance value at the future moment and form the feedforward virtual impedance parameters; Based on the phase characteristics of the circulating current, construct an adaptive feedback control strategy, calculate the virtual impedance value of the real-time response, and generate the feedback virtual impedance parameters.

8. The method according to claim 7, characterized in that, Analyze the complexity characteristics of the circulating current, including: Read the circulating current data and perform the following operations: Segment it into subsequences of length m and m + 1, calculate the proportion of similar subsequences, determine the conditional probability, and obtain the sample entropy of the original signal; where m is a natural number greater than 0; Perform coarse-graining processing on it, calculate the sample entropy for each coarsened sequence, form an entropy-scale curve, and generate multi-scale entropy characteristics; Analyze the sequence pattern and fuzzy boundary, and calculate the permutation entropy and fuzzy entropy; Calculate the change rates of the sample entropy of the original signal, multi-scale entropy characteristics, permutation entropy, and fuzzy entropy, and perform weighted fusion based on this to generate the complexity characteristics of the circulating current.

9. The method according to claim 7, wherein Form the feedforward virtual impedance parameters, including: Based on the dynamic trend prediction data, analyze the dynamic characteristics of the predicted circulating current to determine the key control moments; Based on the energy flow characteristic data, parameter sensitivity data, and complexity characteristics of the circulating current, predict the changes in the energy flow characteristics, parameter sensitivity, and complexity of the circulating current at the future moment to form the predicted characteristic data; For the determined key control moments, apply the predicted characteristic data and the pre-constructed impedance optimization model to solve the optimal virtual impedance value within the prediction time window, generate a virtual impedance time series, and form the feedforward virtual impedance parameters.

10. The method according to claim 7, wherein Generate the feedback virtual impedance parameters, including: Read the phase characteristics of the circulating current, and: Extract the trajectory geometric characteristics, stability indicators, and dynamic evolution characteristics from it to construct the PID control parameters; Extract the parameter sensitivity characteristics and period characteristics from it, adjust the PID control parameters and evaluate the stability indicators to form the characteristic matching parameters, resonance suppression components, and stability enhancement components; generate the basic impedance value through the PID control parameters and the characteristic matching parameters; Dynamically adjust the pre-configured state feedback gain matrix according to it to form the state feedback component; Perform weighted fusion on the basic impedance value, state feedback component, resonance suppression component, and stability enhancement component to form the feedback virtual impedance parameters.

Images