Modular high-voltage SVG circulating current suppression method

Abstract

This invention discloses a modular high-voltage SVG circulating current suppression method, comprising: acquiring and preprocessing module circulating current, voltage, and state data; constructing an enhanced phase space that integrates multiple physical quantities, extracting circulating current phase characteristics, and generating dynamic trend prediction data; calculating optimal feedforward and feedback virtual impedance parameters based on the phase characteristics and trend predictions; and applying a multi-timescale hybrid control strategy to generate module control voltage correction values ​​and PWM control signals. Through topology-aware phase space reconstruction, multidimensional entropy analysis, and stability-guaranteed hybrid control, this method achieves efficient circulating current suppression and improves the stability and robustness of the SVG system.

Description

Technical Field

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

[0002] With the increasing integration of intelligent power grids and renewable energy, high-voltage static VAR generators (SVGs) have become a core device for reactive power compensation in modern power systems, thanks to their fast response, dynamic compensation, and harmonic suppression. In large-capacity applications, SVGs typically employ a modular multi-level structure, expanding capacity and increasing voltage levels by connecting multiple power units in series and parallel. However, circulating currents, caused by inconsistent parameters and varying switching characteristics between modules, have become a key bottleneck restricting the performance and reliability of high-voltage SVGs, severely impacting system safety and power quality. Therefore, developing efficient circulating current suppression technologies is of great engineering value and research significance.

[0003] Current high-voltage SVG circulating current suppression technologies mainly include three categories: the first is a method based on hardware improvement, such as adding balancing inductors and physical damping networks, which limit the circulating current path through additional hardware circuits, but increase system costs and losses; the second is a method based on improved PWM modulation, such as carrier phase shift technology and multi-carrier cooperative sorting strategy, which reduces the generation of circulating current by optimizing switching timing, but has a limited scope of application; the third is a method based on virtual impedance, which simulates impedance characteristics through software algorithms and can suppress circulating current without adding physical devices. Among them, the fixed virtual impedance method and simple feedback regulation method are more widely used, with the advantages of simple implementation and no increase in hardware costs.

[0004] However, existing virtual impedance methods still have significant shortcomings. This is particularly true for high-voltage, large-capacity SVG systems, which have numerous modules, complex topologies, and strong physical coupling. Traditional virtual impedance methods struggle to achieve precise control. First, traditional virtual impedance methods are often based on simplified linear model designs and lack in-depth analysis of the nonlinear dynamic characteristics of SVG, resulting in limited suppression effectiveness under complex operating conditions. Second, existing methods often employ fixed parameters or simple adaptive adjustments, making them ineffective in addressing the rapid changes in circulating current caused by grid disturbances and sudden load changes. Summary of the Invention

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

[0006] The technical solution is a modular high-voltage SVG circulating current suppression method, including:

[0007] Acquire module data, pre-process it, and obtain pre-processed data, including circulating current data, voltage data, and status data;

[0008] Constructing an enhanced phase space based on preprocessed data to generate characteristic data, including circulating current phase characteristics and dynamic trend prediction data;

[0009] Based on the characteristic data and combined with the physical constraints of the system, the virtual impedance parameters are calculated, including feedforward and feedback virtual impedance parameters;

[0010] The virtual impedance parameters are used for gain modulation to generate control signals and output control data, including module control voltage correction values ​​and PWM control signals.

[0011] Beneficial effects: The present invention can more accurately capture the dynamic characteristics of the circulating current and improve the efficiency of circulating current suppression; it can identify precursor signals before the circulating current becomes unstable, start preventive control earlier, and suppress the peak value of the circulating current below the safety threshold, thereby improving the reliability and robustness of the system. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] Figure 1 A flowchart of the steps of a modular high-voltage SVG circulating current suppression method provided in an embodiment of the present application.

[0013] Figure 2 A flowchart of the steps for generating feature data provided in an embodiment of the present application.

[0014] Figure 3 A flowchart of the steps for generating reconstruction parameters provided in an embodiment of the present application.

[0015] Figure 4 A flowchart of the steps for forming an enhanced state vector provided in an embodiment of the present application.

[0016] Figure 5 A flowchart of the steps for generating circulating current phase characteristics provided in an embodiment of the present application. DETAILED DESCRIPTION

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

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

[0019] The study found that the existing control strategy is difficult to take into account both the system dynamic response speed and steady-state suppression accuracy at the same time, and it is easy to cause control system instability when the module parameters differ greatly.

[0020] like Figure 1 As shown, a modular high-voltage SVG circulating current suppression method is proposed, including the following steps:

[0021] Acquire module data, pre-process it, and obtain pre-processed data, including circulating current data, voltage data, and status data;

[0022] Constructing an enhanced phase space based on preprocessed data to generate characteristic data, including circulating current phase characteristics and dynamic trend prediction data;

[0023] Based on the characteristic data and combined with the physical constraints of the system, the virtual impedance parameters are calculated, including feedforward and feedback virtual impedance parameters;

[0024] The virtual impedance parameters are used for gain modulation to generate control signals and output control data, including module control voltage correction values ​​and PWM control signals.

[0025] Specifically, the module data includes module current data, module voltage data, switch state data, DC bus voltage data, and system temperature data. The module data is subjected to noise filtering and data standardization, and preprocessed circulating current data, preprocessed voltage data, and preprocessed state data are output. Using the preprocessed circulating current data, preprocessed voltage data, and preprocessed state data, an enhanced phase space based on the characteristics of the power electronics system is constructed, outputting circulating current phase characteristics and dynamic trend prediction data. Based on the circulating current phase characteristics and dynamic trend prediction data, combined with the system's physical constraints, the optimal virtual impedance parameters are calculated, and 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 implemented, generating coordinated control signals for each module, and outputting module control voltage correction values ​​and PWM control signals.

[0026] According to one aspect of the present application, obtaining preprocessed data includes:

[0027] The raw current data, raw voltage data and switch status data are collected from each module at a set sampling frequency (above 10kHz), and the DC bus voltage data and temperature data are also obtained.

[0028] Switch-state-aware filtering technology is used to process raw current data, identifying switching transients and eliminating the corresponding noise to generate filtered current data. This technology adaptively removes switching noise based on the switching moment and state change characteristics, rather than using traditional fixed-parameter filtering.

[0029] Analyze the current distribution between modules, calculate the deviation of each module current from the average current, and extract the circulating current component;

[0030] The circulating current component is compensated and corrected based on module characteristic differences (including gate voltage, junction temperature, and DC voltage) to obtain a more accurate compensated circulating current. The impact of actual hardware parameter differences on the circulating current is also taken into account.

[0031] Apply a filtering algorithm adapted to the characteristics of the power electronic system to the raw voltage data, perform denoising and standardization, and output preprocessed voltage data;

[0032] The switch status data is sorted and formatted, and a standard data structure containing the relationship between switch status and time is constructed to form pre-processed status data for subsequent analysis and use.

[0033] Specifically, raw current data Iraw, raw voltage data Vraw, and switch state data Sraw are collected from each module at a set sampling frequency fs (above 10kHz). DC bus voltage data Vdc and temperature data Tj are also obtained. Switch-state-aware filtering is applied to the raw current data Iraw to eliminate switching transient noise, resulting in filtered current data If: If(t) = Iraw(t) - Σ α_k·fsw(t-tk, Sraw(tk), Sraw(tk+1)); where tk is the switching instant, fsw is the switching transient model function, α_k is the adaptive weight, and t is the current time. The current difference between modules is calculated to extract the circulating current component: Icirc,i(t) = Ii(t) - (1 / N)·Σ Ij(t); where N is the total number of modules in parallel and Ii(t) is the current of the i-th module. The circulating current component Icirc is compensated based on module parameter differences, resulting in the compensated circulating current: Icirc_comp,i(t)=Icirc,i(t)·(1+ηi(Vgs,i,Tj,i,Vdc)). ηi is the compensation coefficient that takes into account the gate voltage Vgs,i, the junction temperature Tj,i, and the DC voltage Vdc. The raw voltage data Vraw is filtered and normalized to obtain preprocessed voltage data Vp. The switch state data Sraw is converted to a standard format to form preprocessed state data Sp, which contains the specific switching state and timing information for each module.

[0034] like Figure 2 As shown, according to one aspect of the present application, generating feature data includes:

[0035] Analyze the state data and divide the circulating current data into predetermined subsequences according to the switch state to form segmented circulating current data;

[0036] Read the system topology and switch states, dynamically calculate the time delay and embedding dimension for the subsequence, and generate reconstruction parameters;

[0037] Key feature points are extracted from the circulating current data and combined with reconstruction parameters to construct an enhanced state space. The time delay sequence in the circulating current data is fused with physical quantities to form an enhanced state vector. Based on the enhanced state vector, the geometric characteristics of the phase trajectory are calculated to generate the circulating current phase characteristics.

[0038] The coupling relationship of the circulating currents between modules is evaluated to form a module coupling matrix, and the dynamic trend prediction data is generated by combining the phase characteristics of the circulating currents.

[0039] Specifically, the preprocessed state data is analyzed to identify different switching combinations in the system. The compensated circulating current is then segmented into multiple subsequences based on the switching states, generating segmented circulating current data. This segmentation method differs from traditional fixed-window segmentation in that it adaptively segments the system based on the actual operating state. The optimal time delay parameter and embedding dimension are dynamically calculated for each circulating current subsequence based on the current SVG system topology and switching states. The calculation process considers factors such as equivalent inductance, equivalent resistance, number of modules, and number of active modules to generate a reconstruction parameter set. This dynamic parameter selection method overcomes the limitations of traditional fixed-parameter reconstruction. A variable-length adaptive sliding window technique is implemented to determine the minimum necessary analysis time window based on the real-time characteristics of the circulating current and generate adaptive window data. This method dynamically adjusts the analysis window size based on the periodicity and rate of change of the circulating current, improving computational efficiency. Key turning points, extreme points, and rapid change points are intelligently extracted from the compensated circulating current, rather than using all sampled points, to form a set of key feature points. This selective sampling method reduces computational complexity while preserving the key dynamic characteristics of the circulating current. An enhanced state space is constructed, combining the time-delayed sequence of compensated circulating currents with key power electronics quantities (such as voltage and current rate of change, and temperature differences) to form an enhanced state vector. This multi-quantity fusion state space surpasses traditional single-variable reconstruction methods and more comprehensively captures system characteristics. Local geometric properties of the phase trajectory, including trajectory curvature, local stability index, and phase space velocity variation, are calculated. Combined with energy flow characteristics and parameter sensitivity, circulating current phase signatures are generated. These signatures directly reflect the dynamic nature of circulating currents, rather than focusing solely on time- or frequency-domain characteristics. The interrelationships between circulating currents in different modules are evaluated, and the coupling strength and phase differences between modules are calculated to form a module coupling matrix, characterizing the dynamic interactions between modules. Based on the circulating current phase signatures and the module coupling matrix, a recursive prediction technique is used to predict circulating current trends and variations at multiple future time points, generating dynamic trend forecast data. This prediction technique not only considers historical data but also incorporates system physical characteristics, improving prediction accuracy.

[0040] In one embodiment of the present application, the compensation circulating current Icirc_comp is divided into a set of subsequences Icirc_seg related to the switching state according to the preprocessed state data Sp: for the time t, the switching state matrix Sp(t) is defined; the time points with the same or similar switching states are classified: Tk={t|Sp(t)≈Sk}; and a subsequence is generated for each group of time points: Icirc_seg,k={Icirc_comp(t)|t∈Tk}. Based on the SVG topology and the current switching state, the optimal time delay τk and embedding dimension mk for each subsequence are calculated: τk = fτ(Leq(Sk),Req(Sk),freq_char(Icirc_seg,k)); mk = fm(N_module,N_active(Sk), complexity(Icirc_seg,k)). 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 active modules, complexity is the signal complexity assessment function, fτ and fm are the functions used to calculate the time delay and embedding dimension, respectively, Sk is the switching state, and k is the state index. Using a variable-length sequence adaptive sliding window technique, the minimum necessary time window Lwindow is determined: Lwindow(t) = max{Lmin,β·fperiod(Icirc_comp(t),Sp(t))}. Lmin is the minimum window length, fperiod estimates the main period of the current circulating current, and β is the safety factor. A set of key feature points Kpoints is extracted from Icirc_comp instead of using all sampling points: Kpoints={t|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.

[0041] Construct an enhanced state space vector X(t) by incorporating key power electronics 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 preprocessing 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 time, and k is the state index. Calculate the local geometric characteristics of the circulating current phase trajectory to generate the circulating current phase signature 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 characteristic CFPS = [κ(t),λL(t), v(t), energy_flow(t), param_sensitivity(t)]; where energy_flow(t) is the energy flow characteristic, 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(τ)dτ]; C(t) = [C_{i,j}(t)]_{N×N}. The circulating current phase feature CFPS and the coupling strength matrix C(t) are combined to apply the recursive prediction technique to predict the future circulating current trend and generate dynamic trend prediction data DTP: predict the next moment state: X(t+Δt) = F(X(t), CFPS, C(t)); predict the multi-step state: X(t+kΔt) = G(X(t+(k-1)Δt), X(t+(k-2)Δt), ..., X(t)); extract the predicted circulating current: Icirc_pred(t+kΔt) = extract_I(X(t+kΔt)); where F() represents the state function for predicting the next moment, G() represents the multi-step state function, and extract_I() represents the circulating current extraction function. Dynamic trend prediction data DTP = [Icirc_pred(t+Δt), Icirc_pred(t+2Δt), ..., Icirc_pred(t+nΔt), trend_indicators(t)]; where trend_indicators(t) contains key indicators of the circulating current change trend.

[0042] like Figure 3 As shown, according to one aspect of the present application, generating reconstruction parameters includes:

[0043] Extract the topological connection matrix at the current switching moment from the state data, and calculate the equivalent circuit parameters under the current topology based on the topological connection matrix;

[0044] Analyze the frequency characteristics of segmented circulating current data, determine the main frequency component and bandwidth, and calculate the optimal time delay parameters based on equivalent circuit parameters;

[0045] Evaluate the complexity of the segmented circulating current data and determine the optimal embedding dimension based on the total number of system modules and the number of currently active modules;

[0046] The optimal time delay parameter and embedding dimension are associated with the topological connectivity matrix to form the reconstruction parameters.

[0047] Specifically, the pre-processed state data Sp is read, the topological connection matrix Tconn at the current switching moment is extracted, and the connection relationship and switching state of each module with the system are recorded. This matrix encodes the four switching states of each H-bridge module to form a topological description. Based on Tconn, the equivalent circuit parameters under the current topology are calculated, including the module equivalent inductance Leq, equivalent resistance Req, and equivalent capacitance Ceq. The calculation process takes into account the series and parallel relationship between modules. Unlike the traditional fixed parameter method, these parameters change in real time with the switching state. The frequency characteristics of the segmented circulating current data Icirc_seg are analyzed, and the main frequency component f_main, bandwidth BW, and harmonic distribution HarmDist are calculated. Different from the traditional FFT analysis, short-term adaptive spectrum estimation is used here to adapt to the non-stationary characteristics of the circulating current. According to the principles of information theory, the optimal time delay parameter τk is calculated based on the main frequency component and bandwidth to ensure maximum independence of data points in phase space: for high-frequency-dominated circulating currents (f_main>f_threshold), a shorter delay is set to capture rapid changes; for low-frequency-dominated circulating currents (f_main ≤ f_threshold), a longer delay is set to capture the complete cycle characteristics; the delay parameter is dynamically adjusted: τk = α·(1 / f_main)·(Leq / Leq_ref) β , where α and β are parameters preset according to the system response characteristics, and f_threshold is the threshold of the main frequency component.

[0048] The complexity of the segmented circulating current data Icirc_seg is evaluated, and the complexity index CI is calculated, including waveform complexity, irregularity, and entropy. Unlike traditional signal complexity assessments, this method specifically focuses on the impact of switching transitions and state transitions in power electronics systems on complexity. The optimal embedding dimension mk is determined based on the total number of system modules N_module, the number of currently active modules N_active, and the complexity index CI. The base dimension is determined as: m_base = N_active + 1, accounting for inter-module degrees of freedom. Complexity compensation is performed as: m_comp = ┌CI·γ┐, where γ is the compensation coefficient and ┌ ┐ indicates rounding up. The final embedding dimension is: mk = min(m_base + m_comp, m_max), ensuring that the maximum allowed dimension m_max is not exceeded. The calculated τk and mk are associated with the current topological state Tconn and stored as a reconstruction parameter set RP, which contains the time delay, embedding dimension, and the corresponding topological state identifier. This parameter set will be used in subsequent phase space construction. Dynamically evaluate the performance of the reconstruction parameter set RP, calculate the parameter sensitivity index PS, and make fine-tuning when 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 changing trend of the parameter sensitivity index PS. If PS continues to increase, it indicates that the current reconstruction parameters are gradually becoming incompatible with system changes, triggering parameter recalculation: If PS>α·PS_previous and persists for N cycles, then: re-execute the above steps to update the reconstruction parameter set RP; use the weighted average of the new and old parameters for a smooth transition; output the optimized reconstruction parameter set RP_opt for subsequent phase space construction. PS_previous is the parameter sensitivity index of the previous cycle.

[0049] like Figure 4 As shown, according to one aspect of the present application, forming an enhanced state vector includes:

[0050] Construct a basic time delay embedding vector based on circulating current data and reconstruction parameters; extract key voltage features from voltage data; calculate the dynamic characteristics of circulating current data; and calculate the temperature characteristics of key power devices, including junction temperature difference and temperature change rate, based on pre-stored system temperature data.

[0051] The basic time delay embedding vector is weightedly fused with key voltage characteristics, dynamic characteristics and temperature characteristics through a multi-physical quantity weight fusion matrix to form an enhanced state vector.

[0052] Specifically, read the pre-processed circulating 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 a basic time delay embedding vector Xbase, which contains the current and historical circulating current data: extract the values ​​of the current time t and the past (mk-1) delay times from Icirc_comp to form a 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 circulating current. Read the pre-processed voltage data Vp and extract key voltage features, including the module output voltage Vout, the voltage change rate dV / dt and the voltage imbalance Vunbal. Different from traditional methods, here we pay special attention to the correlation between voltage characteristics and circulating current. Calculate the dynamic characteristic indicators of the pre-processed circulating current data Icirc_comp, including the change rate dI / dt, acceleration d 2 I / dt 2 and its changing trend Itrend. These derivative characteristics reflect the dynamic behavior of the circulating current, a feature often overlooked in traditional state space. System temperature data is read to calculate the junction temperature difference ΔTj and temperature change rate dTj / dt between key power devices. Electrothermal coupling effects are incorporated into the state space to extract the dynamic characteristic parameters of the switching devices, including switching delay Tdelay, rise / fall times Tr / Tf, and on-state voltage Von. These reflect the actual switching behavior of the power devices and influence the circulating current distribution.

[0053] Construct a multi-physical quantity weight fusion matrix W_fusion, and adaptively adjust the weights of each physical quantity according to the current working conditions: Under steady-state conditions: increase the weights of temperature and voltage imbalance; Under dynamic conditions: increase the weights of current change rate and acceleration; Under disturbance conditions: increase the weights of switch characteristic parameters. This adaptive weight distribution mechanism is not available in traditional state space methods. The basic delay embedding vector Xbase is weightedly fused with voltage characteristics, current dynamic characteristics, temperature characteristics and switch characteristics to construct an 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; W_fusion is applied to perform weighted fusion of multiple physical quantities, outputting a complete enhanced state vector X(t) as the basis for subsequent feature extraction and prediction. The enhanced state vector X(t) is normalized to ensure dimensional consistency among different physical quantities, and a standardized enhanced state vector Xnorm(t) is output. This normalization takes into account the typical ranges and dynamic characteristics of different physical quantities.

[0054] like Figure 5 As shown, according to one aspect of the present application, generating a circulating current phase feature includes:

[0055] Read and, based on the enhanced state vector, sequentially: calculate the trajectory curvature and torsion to form the trajectory geometric characteristics; evaluate the local stability of the phase trajectory to form a stability index; calculate the phase space velocity vector and acceleration vector to form the dynamic evolution characteristics;

[0056] Based on the circulating current data and state data, energy flow characteristics are formed;

[0057] Extract the sensitivity of circulating current data to module parameter differences to form parameter sensitivity characteristics;

[0058] The weights of trajectory geometric characteristics, stability index, dynamic evolution characteristics, energy flow characteristics and parameter sensitivity characteristics are determined to form the circulating current phase characteristics.

[0059] Specifically, a series of normalized augmented state vectors Xnorm(t) are read to construct a phase trajectory Xtraj, which contains state evolution information at consecutive time points. This trajectory captures the dynamic behavior of the system in the augmented state space. Local geometric features of the phase trajectory are calculated to analyze its shape and curvature: A continuous sequence of points in the phase trajectory is read to calculate the trajectory curvature κ(t), which characterizes the curvature and reflects the nonlinearity of the system's state changes; the trajectory distortion τ(t) is calculated to characterize the degree of distortion in high-dimensional space; and the trajectory geometric features GF are output. These features surpass traditional linear analysis methods and can capture the nonlinear dynamics of the system. The local stability of the phase trajectory is evaluated to quantify the system's sensitivity to perturbations: a reference point and its neighboring points are selected in the phase trajectory, the evolution trajectories of these points are tracked, and the trajectory separation index is calculated. The local Lyapunov exponent λL(t) is extracted to characterize the local stability of the system; and the stability index SI is output. These indices can predict the system's behavior under perturbations and provide forward-looking information for control.

[0060] Calculate the time evolution characteristics of the phase trajectory to capture the rate of change of the system state: read the state of continuous time points in the phase trajectory, calculate the phase space velocity vector v(t)=dX(t) / dt, and characterize the rate of change of the state; calculate the phase space acceleration vector a(t)=d 2 X(t) / dt 2, characterizing the acceleration characteristics of state changes; outputting dynamic evolution features DF, which reflect the time-varying pattern of the system state. Analyzing the periodicity and reproducibility of phase trajectories to identify regular patterns in the system: constructing a delay coordinate diagram to analyze the self-similarity of phase trajectories; calculating recurrence graph indicators to quantify the reproducibility of trajectories; extracting the periodic intensity PI and the cyclic pattern CP; outputting periodic features PF, which help identify regular components in the circulating current. Calculating the energy flow characteristics of the circulating current within the switching cycle: reading preprocessed circulating current data Icirc_comp and preprocessed state data Sp, dividing the energy analysis window according to the switching cycle, calculating the circulating current energy Ecirc(t) and power flow Pcirc(t); determining the energy flow direction indicator; and outputting the energy flow features EF, which establish a connection between phase dynamics and physical energy. Extract the sensitivity of circulating current to module parameter differences and analyze the impact of parameter changes on the system: Build a parameter change response model based on the enhanced state vector X(t) and calculate the sensitivity to resistance differences SR(t), inductance differences SL(t), and capacitance differences SC(t). Comprehensively form the parameter sensitivity vector S(t) and output the parameter sensitivity features SF, which reflect the response characteristics of the circulating current to system parameter changes.

[0061] All of the above features are integrated to construct a complete circulating current phase signature CFPS: the trajectory geometry feature GF, stability index SI, dynamic evolution feature DF, period feature PF, energy flow feature EF and parameter sensitivity feature SF are read; the weight of each feature is determined through feature importance assessment; the final circulating current phase signature CFPS is synthesized, which contains comprehensive information on the circulating current dynamics; this feature set will be used for subsequent trend prediction and control parameter calculation.

[0062] According to one aspect of the present application, generating dynamic trend prediction data includes:

[0063] Based on the enhanced state vector and circulating current phase characteristics, a point set similar to the current state is selected from the historical data to build a local linear prediction model;

[0064] Based on the local linear prediction model, the current state is used to predict the state at the next moment, and the prediction results are used as new initial conditions to predict the further future, generating a multi-step prediction sequence;

[0065] Circulating current prediction data are extracted from the multi-step prediction sequence, and the key dynamic indicators of the predicted circulating current are calculated and integrated into dynamic trend prediction data.

[0066] Specifically, the enhanced state vector X(t) and circulating current phase characteristics CFPS of the historical time series are read to construct the initial conditions for the prediction. This ensures that the prediction is based on the complete state information of the system. A local linear prediction model in phase space is constructed to build a linear approximation near the current state: a set of points Xsim similar to the current state is selected from the historical data; based on the evolution trajectory of these points, a local linear map Flocal is constructed; and a local prediction model LPM is output, which adapts to the local linearization of the nonlinear characteristics of the circulating current. A recursive rolling prediction strategy is implemented to gradually extend the prediction range: the current state X(t) and the local prediction model LPM are read; the next state X(t+Δt) = Flocal(X(t)) is predicted; the prediction results are used as new initial conditions to predict the further future X(t+2Δt); this process is repeated to generate a multi-step prediction sequence Xpred. This recursive prediction method outperforms traditional single-step prediction and can capture long-term trends. An error correction mechanism based on phase characteristics is introduced to improve the accuracy of recursive predictions: the systematic error patterns in historical predictions are analyzed; the prediction correction factor CF is adjusted based on the circulating current phase characteristics CFPS; the correction factor is applied to optimize the prediction results: Xpred_corr = Xpred + CF; and the corrected prediction sequence Xpred_corr is output, which reduces the error accumulation in the recursive prediction.

[0067] A multi-model prediction fusion approach combines the strengths of different prediction strategies: multiple prediction models, including a local linear model, a nonlinear dynamic model, and a pattern matching model, are run in parallel. The weights Wmodel of each model are dynamically adjusted based on historical accuracy and the current circulating current phase characteristics (CFPS). The prediction results of each model are weighted and fused to form an integrated prediction (Xpred_ensemble). This multi-model fusion approach avoids the limitations of a single model and improves prediction robustness. Circulating current trend information is extracted from the predicted state sequence: the circulating current prediction component Icirc_pred is extracted from Xpred_ensemble. The predicted circulating current trend is analyzed, and the peak value Ipeak, root mean square value Irms, and frequency characteristic Ifreq are calculated. The reliability of the prediction is evaluated, generating a confidence index (Conf). The circulating current prediction data (IP) is output, containing the predicted amplitude and characteristics of the future circulating current. Key dynamic indicators of the predicted circulating current, including the maximum rate of change (max(dI / dt), energy change trend (ΔE), and stability expectation (Stability), are calculated to form trend indicator data (TI). The circulating current prediction data (IP) and trend indicator data (TI) are integrated to form a complete dynamic trend prediction data (DTP), providing forward-looking information for subsequent virtual impedance parameter calculation.

[0068] According to one aspect of the present application, calculating virtual impedance parameters includes:

[0069] Analyze the energy distribution and flow characteristics of circulating current data within the switching cycle to generate energy flow characteristic data;

[0070] Evaluate the sensitivity of circulating current data to module parameter differences to form parameter sensitivity data;

[0071] Calculate the sample entropy and entropy change rate of the enhanced phase space and analyze the complexity characteristics of the circulating current;

[0072] The circulating current complexity characteristics are combined with 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;

[0073] An adaptive feedback control strategy is constructed based on the phase characteristics of the circulating current, the real-time response virtual impedance value is calculated, and the feedback virtual impedance parameters are generated.

[0074] Specifically, the system comprehensively analyzes preprocessed circulating current data, preprocessed voltage data, and the module coupling matrix to calculate a vector of operating condition indicators representing the current system operating state. These indicators include load change rate, grid fluctuation, active power, reactive power, harmonic distortion, and module coupling, generating operating condition indicator data. The system monitors the changing trends of these indicators. When significant changes are detected (such as sudden load changes or grid disturbances), these indicators are marked as operating condition transition points, triggering corresponding adjustments to the control strategy and outputting an operating condition transition flag. This proactive operating condition identification method enables proactive response to system state changes. The system analyzes the energy distribution and flow characteristics of the circulating current within the switching cycle, calculates the circulating current energy and its rate of change, determines the direction and intensity of energy flow, and generates energy flow characteristic data. This data reflects the dynamic energy characteristics of the circulating current and serves as an important basis for optimizing virtual impedance. The system also assesses the sensitivity of the circulating current to inter-module parameter differences, calculating the sensitivity to resistance, inductance, and capacitance differences, generating parameter sensitivity data. This sensitivity analysis helps identify the primary factors influencing the circulating current and optimize control parameters accordingly. Complexity metrics of the circulating current phase space, including sample entropy and entropy change rate, are calculated to extract signal complexity features. These metrics reflect the nonstationarity and uncertainty of the circulating current and guide the dynamic adjustment of the control strategy. Based on dynamic trend prediction data, energy flow characteristics, parameter sensitivity data, and signal complexity characteristics, the optimal virtual impedance value for the future is comprehensively calculated to form a feedforward virtual impedance parameter. This parameter is optimized for the predicted circulating current trend and provides preventive control. An adaptive feedback control strategy is designed based on the current circulating current phase characteristics to calculate the real-time virtual impedance value and generate the feedback virtual impedance parameter. This parameter is responsible for real-time adjustment of the current circulating current state to ensure control accuracy. Taking into account the physical limitations of power electronic devices, safety constraints are imposed on the virtual impedance parameter to ensure that the parameter value and rate of change are within the system's tolerance range. The constrained virtual impedance parameter is then output. This ensures the feasibility of the control parameters and the safe operation of the system.

[0075] In one embodiment of the present application, based on the pre-processed circulating current data, pre-processed 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)]; wherein C_norm(t) is the norm of the coupling matrix, which characterizes the degree of coupling between modules. Detect significant changes in operating conditions and trigger control strategy adjustments when necessary: ​​if||W(t) - W(t-Δt)||>θW·||W(t-Δt)||: mark as the operating condition transition point, and adjust the control parameter calculation strategy; wherein θW is the operating condition change threshold. Calculate the circulating current energy characteristic Ecirc(t) and power flow characteristic Pcirc(t) within the switching cycle: Ecirc(t) = ∫ t t+Tsw L(τ)·Icirc_comp 2 (τ)dτ / 2; Pcirc(t)=dEcirc(t) / dt; Energy flow direction indicator: Denergy(t)=sign(Pcirc(t))·|Pcirc(t)| γWhere Tsw is the switching period, L(τ) is the time-varying equivalent inductance, and γ is the energy exponent parameter. Extract the sensitivity vector S(t) of the circulating current to 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 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 points with a distance less than r in the m+1 and m-dimensional phase spaces, respectively. Based on the dynamic trend prediction data DTP, energy characteristics Ecirc(t), sensitivity vector S(t) and entropy change rate dSampEn(t), the feedforward virtual impedance parameter Zv_ff is calculated: for 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 time; k is the prediction step; Δt is the time interval; h is the mapping function, which maps the predicted circulating current characteristics to the optimal virtual impedance value. Based on the current circulating current phase characteristics (CFPS) and system state, the feedback virtual impedance parameter Zv_fb is calculated as follows: Zv_fb(t)=KP·Icirc_comp(t)+KI·∫Icirc_comp(τ)dτ+KD·dIcirc_comp(t) / dt; where KP is the proportional coefficient, KI is the integral coefficient, and KD is the differential coefficient. It is dynamically adjusted based on the CFPS and the operating condition indicator W(t); Icirc_comp is the compensation circulating current; t is time; and τ is the integral variable. The physical limitations of power electronics devices are applied to ensure that the virtual impedance parameter meets the realizability constraints: |Zv_ff(t)|≤Zmax(t); |dZv_ff(t) / dt|≤η(t)·Vdc / (Irated·Tsw); where Zmax(t) is the maximum allowable impedance, η(t) is the dynamic safety factor, Vdc is the DC bus voltage, and Irated is the rated current.

[0076] According to one aspect of the present application, the process of generating operating condition indicator data specifically includes: reading preprocessed circulating current data, calculating the circulating current amplitude, phase, frequency, and waveform factor, and forming basic circulating current characteristics. Analyzing preprocessed voltage data, extracting grid voltage variation characteristics, including the voltage change rate dVgrid / dt, fluctuation index, and imbalance, and outputting grid status characteristics. Reading preprocessed status data and the module coupling matrix, evaluating the coupling state between modules: calculating the norm of the coupling matrix to characterize the overall coupling strength; extracting the maximum coupling value and the corresponding module pair; analyzing the uniformity of the coupling distribution and calculating the coupling imbalance index; and outputting module coupling characteristics, which reflect the interaction between modules within the system. Calculating system power characteristics, including active power Pactive, reactive power Qreactive, and fluctuating power Sosc: reading preprocessed voltage data Vp and preprocessed circulating current data Icirc_comp; calculating instantaneous power S(t) = Vp(t)·Icirc_comp(t); decomposing the active, reactive, and fluctuating components; and outputting power characteristic data PD, which characterizes the energy transfer state of the system. Analyze the system's harmonic state and calculate total harmonic distortion (THD) and major harmonic components (Harm_n): Harmonic analysis is performed on preprocessed circulating current data; the ratio of each harmonic amplitude to the fundamental is calculated; the dominant harmonics and their phase relationships are determined; and harmonic characteristic data (HD) is output to characterize the system's nonlinear distortion state. Evaluate system load characteristics, including load change rate (dIload / dt), load type indicators, and load stability: Analyze load-related patterns in preprocessed circulating current data; determine load type based on power factor and harmonic characteristics; assess the regularity and predictability of load changes; and output load characteristic data (LD) to characterize the system's external load conditions. Construct a working condition evaluation weight matrix (W_eval) to dynamically adjust the weights of various features based on system operational priorities: In high-precision circulating current control mode, the weights of circulating current characteristics and module coupling characteristics are increased; in fast dynamic response mode, the weights of load characteristics and grid status characteristics are increased; and in harmonic suppression mode, the weights of harmonic characteristic data are increased. This adaptive weight adjustment mechanism is not available in traditional working condition evaluation methods. Integrate all features to construct the 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; apply W_eval for weighted fusion; output the final operating condition index data W(t), which comprehensively represents the current operating condition of the system; this index will be used for subsequent control parameter calculation and strategy adjustment.

[0077] According to one aspect of the present application, analyzing the complexity characteristics of circulating current includes:

[0078] To read the circulating current data, perform the following operations:

[0079] The circulating current data is segmented into subsequences of length m and m+1, the proportion of similar subsequences is calculated, the conditional probability is determined, and the entropy of the original signal sample is obtained; where m is a natural number greater than 0;

[0080] The circulating current data is coarse-grained, and the sample entropy of each coarse-grained sequence is calculated to form an entropy-scale curve and generate a multi-scale entropy feature;

[0081] Based on circulating current data, the sequence patterns and fuzzy boundaries are analyzed, and the permutation entropy and fuzzy entropy are calculated;

[0082] The change rates of the original signal sample entropy, multi-scale entropy features, permutation entropy and fuzzy entropy are calculated, and the entropies are weighted and fused based on the change rates to generate the circulating current complexity features.

[0083] Specifically, preprocessed circulating current data Icirc_comp is read within a time window to ensure that it contains sufficient dynamic information. The time window length is adaptively adjusted based on system characteristics, rather than the traditional fixed window length. Sample entropy is calculated to quantify the complexity and irregularity of the time series: Icirc_comp is segmented into subsequences of length m and m+1; the proportion of similar subsequences is calculated to determine the conditional probability; the sample entropy value SampEn = -ln(A / B), where A is the m+1-dimensional similarity ratio and B is the m-dimensional similarity ratio, is calculated; and the sample entropy value is output to characterize the randomness and unpredictability of the circulating current. Multiscale entropy analysis is implemented to assess complexity at different time scales: coarse-graining processing with different scaling factors τ is applied to Icirc_comp; sample entropy is calculated for each coarse-grained sequence; a multiscale entropy curve is generated; and the multiscale entropy feature MSE is output to characterize the complex dynamics of the circulating current at different time scales. Calculate the permutation entropy to capture the sequential pattern of the time series: split Icirc_comp into overlapping subsequences; analyze the relative size order of the data points in each subsequence; count the frequency of occurrence of different permutation patterns; calculate the permutation entropy PermEn=-Σ p(π)·log(p(π)), where p(π) is the probability of permutation pattern π; output the permutation entropy value PermEn, which represents the sequence complexity of the circulating current.

[0084] Fuzzy entropy analysis is performed to address noise and fuzzy boundaries in the circulating current: A fuzzy membership function is applied to Icirc_comp to calculate sequence similarity; a fuzzy entropy value is calculated based on the fuzzy similarity; and the fuzzy entropy value, FuzzyEn, is output, providing a more noise-robust complexity assessment. The entropy change rate dEntropy is calculated to characterize the dynamic characteristics of system state transitions: The time variations of SampEn, MSE, PermEn, and FuzzyEn are tracked; the change rates of each entropy value (dSampEn, dMSE, dPermEn, and dFuzzyEn) are calculated; and the entropy change rate vector dEntropy is synthesized to represent the dynamic changes in system complexity. Entropy gradient analysis is performed to identify system state transition points: The spatial and temporal gradients of dEntropy are calculated; points with significant gradient changes are identified and marked as transition points; the intervals near the transition points are analyzed in detail; and state transition data ST is output, containing the location and characteristics of the transition points. Integrate various entropy indicators to construct a comprehensive complexity evaluation model: design entropy feature fusion weights and adjust the weights of various entropy indicators according to the current circulating current characteristics; fuse SampEn, MSE, PermEn, FuzzyEn and dEntropy; generate a comprehensive complexity index CI; output signal complexity feature CI to comprehensively characterize the complex dynamic characteristics of the circulating current; and use it for subsequent control parameter optimization.

[0085] According to one aspect of the present application, forming a feedforward virtual impedance parameter includes:

[0086] Predict the dynamic characteristics of circulating current based on dynamic trend prediction data analysis and determine the critical control moments;

[0087] Based on the energy flow characteristic data, parameter sensitivity data and circulating current complexity characteristics, the change amount of energy flow characteristics, parameter sensitivity and circulating current complexity at a future moment is predicted to form prediction characteristic data;

[0088] For the determined critical control moments, the predicted characteristic data and the pre-built impedance optimization model are applied to solve the optimal virtual impedance value within the predicted time window, generate the virtual impedance time series, and form the feedforward virtual impedance parameters.

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

[0090] The signal complexity characteristic CI is read to analyze the nonstationary characteristics of the circulating current: the entropy change CI_pred at future moments is predicted; possible state transition points are identified; control strategies are designed for different complexity intervals; and the predicted complexity characteristic PC is output to provide a basis for managing nonstationary circulating currents. A mapping model is constructed to link the virtual impedance parameters to the circulating current: Based on historical data, a relationship model Z_I_model is established between the virtual impedance and the circulating current response. This model considers the nonlinear mapping relationship under different operating conditions, surpassing the traditional linear impedance model. The output impedance mapping model ZIM provides a basis for solving the optimal parameters. A multi-objective optimization function is constructed to balance circulating current suppression and dynamic response: the circulating current suppression objective function J_circ = f(Icirc_pred); the dynamic response objective function J_resp = g(dIcirc_pred / dt); and a comprehensive objective function J = α·J_circ + (1-α)·J_resp is constructed, where α is a weight dynamically adjusted according to the operating conditions. The output optimization objective function OF guides the solution of the optimal parameters. Apply a predictive optimization algorithm to solve the optimal virtual impedance parameters for future moments: read the impedance mapping model ZIM and the optimization objective function OF; solve the virtual impedance value corresponding to the minimum value of the objective function for each key moment in the prediction time window; generate a virtual impedance time series Z_seq; apply smoothing constraints to ensure the continuity of impedance changes; and output the feedforward virtual impedance parameter Zv_ff, which contains the optimal impedance values ​​for multiple future moments.

[0091] According to one aspect of the present application, generating feedback virtual impedance parameters includes:

[0092] Read the circulating current phase characteristics and:

[0093] Extract trajectory geometry, stability index and dynamic evolution characteristics from circulating current phase characteristics to construct PID control parameters;

[0094] Extract parameter sensitivity characteristics and period characteristics from the circulating current phase characteristics, adjust PID control parameters and evaluate stability indicators to form characteristic matching parameters, resonance suppression components and stability enhancement components; generate basic impedance values ​​through PID control parameters and characteristic matching parameters;

[0095] Dynamically adjust a pre-configured state feedback gain matrix according to the phase characteristics of the circulating current to form a state feedback component;

[0096] The basic impedance value, state feedback component, resonance suppression component and stability enhancement component are weighted and fused to form the feedback virtual impedance parameter.

[0097] Specifically, the current circulating current phase signature (CFPS) is read to extract key features required for real-time control, including the trajectory geometry (GF), stability index (SI), and dynamic evolution (DF). Based on these phase signatures, an adaptive PID controller is designed to dynamically adjust control parameters. The proportional coefficient (KP) is adjusted based on the trajectory curvature (κ(t)). When the curvature is large, KP is increased to improve response speed, while when the curvature is small, KP is decreased to avoid overshoot. The integral coefficient (KI) is adjusted based on the local Lyapunov exponent (λL(t)). When stability is high, KI is increased to eliminate steady-state errors, while when stability is low, KI is decreased to prevent integral saturation. The differential coefficient (KD) is adjusted based on the phase space velocity (v(t). When the velocity is high, KD is increased to provide damping, while when the velocity is small, KD is decreased to avoid noise amplification. Finally, the adaptive PID parameter (KPID) is output. This phase-characteristics-based parameter adjustment mechanism differs from traditional error-driven parameter tuning. The parameter sensitivity characteristics SF and period characteristics PF in the circulating current phase characteristics CFPS are read, and the control parameters are fine-tuned according to different circulating current characteristics: the high-sensitivity parameters are optimized: when the inductance sensitivity SL(t) is greater than the threshold, the ratio of KP and KD is adjusted to improve the inductance control accuracy; the control frequency is adjusted according to the periodicity intensity PI: when the periodicity is strong, the control update frequency and the circulating current period are synchronized; the control parameters are preset based on the circulation pattern CP: when the typical pattern is identified, the pre-optimized parameter combination is applied; the characteristic matching parameter FMP is output to achieve more precise parameter adjustment. The feedback virtual impedance base value is constructed to achieve immediate response to the current circulating current: the basic PID control law: Zv_base = KPID(0)·Icirc_comp(t) + KPID(1)·∫Icirc_comp(τ)dτ + KPID(2)·dIcirc_comp(t) / dt; the characteristic matching parameter FMP is applied for fine-tuning: Zv_base = Zv_base·FMP; the base impedance value Zv_base is output as the basic component of feedback control.

[0098] The enhanced state vector X(t) is read and a state feedback enhancer is designed to supplement basic PID control: a subset of state variables directly related to the circulating current, Xsub, is extracted; the state feedback control variable, Zv_state, is calculated: Zv_state = K_state · Xsub, where K_state is the state feedback gain matrix; K_state is dynamically adjusted based on the circulating current phase signature, CFPS, to achieve adaptive state feedback; and the state feedback component, Zv_state, is output to enhance direct response to system conditions. A resonance suppression compensator is introduced to provide additional damping for resonant components in the circulating current: the periodic signature, PF, in the circulating current phase signature, CFPS, is analyzed to identify potential resonant frequencies, f_res; a targeted resonance suppressor, Zv_res, is designed: Zv_res = A_res sin(2πf_res·t + φ_res); the suppressor amplitude, A_res, and phase, φ_res, are dynamically adjusted to create a 180° phase difference with the resonant components of the circulating current; the resonance suppression component, Zv_res, is output to effectively suppress circulating current oscillations of a specific frequency. Based on the stability index SI in the circulating current phase characteristic CFPS, a stability enhancer is designed. When the local Lyapunov exponent λL(t) approaches a critical value, an additional stabilizing component is added. The stabilizing control variable Zv_stab is calculated as: K_stab · (λL_threshold - λL(t)). When λL(t)>λL_threshold, the stability enhancement component Zv_stab is output to prevent the system from becoming unstable under critical conditions. Integrate each control component to generate the final feedback virtual impedance parameters: read the base 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 weight of each component W_fb based on the circulating current phase characteristics 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.

[0099] According to one aspect of the present application, outputting control data includes:

[0100] A Lyapunov energy function based on circulating current and control weights is constructed to evaluate system stability and control effectiveness, generating stability assessment data. This energy-based stability assessment method can theoretically guarantee the stable operation of the system.

[0101] A stability-oriented gain modulation mechanism is constructed to dynamically adjust the control weight according to the circulating current state and virtual impedance difference, ensuring that the system evolves in a stable direction and outputting the basic hybrid weight.

[0102] Implementing a multi-timescale control strategy, designing control layers with different response speeds for high-frequency disturbances, medium-frequency load changes, and low-frequency parameter drifts, calculating the weights of each layer and combining them into a final hybrid weight to form a multi-scale hybrid weight. This hierarchical control structure can simultaneously respond to system changes at different frequencies.

[0103] Based on multi-scale hybrid weights, the feedforward virtual impedance parameters and feedback virtual impedance parameters are weighted and fused to calculate the final integrated virtual impedance value. This feedforward-feedback hybrid control structure combines the advantages of predictive control and real-time feedback.

[0104] Based on the comprehensive virtual impedance value, the voltage compensation amount required for each module is calculated to form a preliminary module voltage correction value;

[0105] Considering the power balance and thermal stress distribution among modules, the module voltage correction value is coordinated and optimized to ensure that each module maintains reasonable power distribution while suppressing the circulating current and outputs a coordinated voltage correction value;

[0106] Combine the coordinated voltage correction value with the reference voltage to generate the final control voltage instruction of each module and form the module control voltage correction value;

[0107] Based on the module control voltage correction value, an appropriate PWM modulation strategy is applied, and factors such as switching frequency and dead time are considered to generate a PWM control signal to drive the power device.

[0108] Specifically, a Lyapunov energy function based on circulating current and control weight is constructed to evaluate system stability and control effect: the system state vector X(t)=[Icirc_comp,ω] is defined, which includes the compensation circulating current and feedforward-feedback hybrid weight; and the Lyapunov candidate function is constructed: 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 between modules; ω is the feedforward-feedback hybrid weight vector, ω* is the ideal weight vector; kω is the weight adjustment coefficient matrix, and Q is the positive definite weight matrix, reflecting the energy cost of weight deviation; T Represents the transpose of a vector or matrix. Calculates 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 that dV / dt<0: Construct a stability indicator SI=-dV / dt. When SI is positive and greater than the threshold, the system evolves in a stable direction; output stability evaluation data SE={V(t),dV / dt, SI} as a direct basis for weight modulation. Construct a gain modulation mechanism based on Lyapunov stability to ensure stable convergence of the system: read the stability evaluation data, and design a dynamic weight adjustment law based on Lyapunov stability theory to ensure that dV / dt<0: dω / dt=-γ·▽ω(dV / dt); where γ is a positive learning rate parameter, dynamically adjusted to accelerate convergence; ▽ω(dV / dt) is the gradient of dV / dt with respect to ω, pointing to the direction in which 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 influence 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| α α is an exponential parameter that is dynamically adjusted based on the circulating current state. A stability constraint is applied. When the stability indicator SI falls below a threshold, a damping term is added: 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; and β is the damping coefficient. The basic hybrid weight ω(t) is output as the initial weight for multi-timescale control.

[0109] Implement a multi-timescale control strategy and perform hierarchical weight optimization based on the stability evaluation data SE: read the basic mixing weights and 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); and construct the corresponding hierarchical Lyapunov function: 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 disturbances, 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, ensuring the fastest high-frequency response speed. The weights of each frequency band are combined to form a multiscale 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, and ω_lf are the optimal weights for each frequency band. The global Lyapunov function stability under the combined weights is verified: dV_total / dt = dV_hf / dt + dV_mf / dt + dV_lf / dt<0. The final multiscale hybrid weight ω_ms(t) is output for feedforward and feedback virtual impedance fusion. Based on multiscale mixing weights, the feedforward and feedback virtual impedances are fused to ensure global stability: the multiscale mixing weight ω_ms(t), the feedforward virtual impedance parameter Zv_ff, and the feedback virtual impedance parameter Zv_fb are read; weighted fusion is applied: Zv(t) = ω_ms(t)·Zv_fb(t)+(1-ω_ms(t))·Zv_ff(t); where Zv is the integrated virtual impedance value and ω_ms is the multiscale mixing weight. The Lyapunov function change rate of the fused impedance is calculated: dV / dt|Zv = Icirc_comp T ·L·(dIcirc_comp / dt)|Zv; Verify the stability condition dV / dt|Zv<0. If not met, 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 integrated virtual impedance value Zv(t) to ensure system stability and dynamic performance.

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

[0111] In another embodiment of the present application, the output control data may also be: constructing a Lyapunov energy function for stability assurance: V(Icirc_comp,ω) = (1 / 2)L·Icirc_comp 2 + (1 / 2)kω·(ω-ω*) 2where ω is the feedforward-feedback mixing weight, ω* is the ideal weight, and kω is the weight adjustment coefficient. A gain modulation law based on Lyapunov stability is constructed and the mixing weight ω(t) is calculated as: dω(t) / dt = -γ·sign(Icirc_comp·d(Zv_fb-Zv_ff) / dt) - β·ΨV / Ψω; where γ is the adjustment rate and β is the damping coefficient. A multi-timescale gain modulation strategy is applied to calculate the fast layer weight ωf(t), the medium 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 comprehensive 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. The final integrated virtual impedance value Zv is calculated. Based on this value, the voltage correction value ΔVi for each module is calculated: ΔVi(t) = -Zv(t) · Icirc_comp,i(t). Considering inter-module power balance, the voltage correction values ​​are 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)). Here, 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. The final module control voltage command Vi_cmd is generated: Vi_cmd(t) = Vref(t) + ΔVi_coord(t), where Vref(t) is the reference voltage. According to the module control voltage instruction 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.

[0112] In a specific embodiment of the present application, based on a three-phase four-module high-voltage SVG system, the total capacity is 10MVA, the DC side voltage is 3kV, and the AC side phase voltage is 10kV / 50Hz. The system includes 4 series-connected H-bridge modules per phase, and parameter differences between 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 = 10kV / 50 Hz; number of modules: 4 series modules per phase, a total of 12 modules; single module inductance: Lm = 2.0 mH ± 10%; single module capacitance: Cm = 5000 μF ± 7%; sampling frequency: fs = 20 kHz; switching frequency: fsw = 1000 Hz. The specific steps are:

[0113] Step 1: Data collection and preprocessing.

[0114] 1.1. Collect raw data from the four modules of phase A at a sampling frequency of 20kHz and an acquisition time of 100ms (a total of 2000 data points): Raw current data Iraw: [23.5A, 24.2A, 23.9A, 24.7A] (instantaneous values ​​of the four modules at t=10ms); Raw voltage data Vraw: [753V, 745V, 750V, 748V] (instantaneous values ​​of the four modules at t=10ms); Switch 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=10ms); DC bus voltage data Vdc: [3050V, 3030V, 3070V, 3020V] (DC side voltages of the four modules); Temperature data Tj: [45.2°C, 47.5°C, 44.8°C, 48.3°C] (four module IGBT junction temperatures).

[0115] 1.2. Switching state-aware 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 raw current data; tk is the switching instant; fsw is the switching transient model function; α_k is the adaptive weight; and t is time. Taking the switching transient at t = 10.005 ms as an example, the switching state of module 1 is detected to change 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 raw current Iraw(10.005 ms) peaks at 30.2 A; the filtered current If(10.005 ms) = 30.2 A - 0.85 × 2.7 A = 27.905 A.

[0116] 1.3. Circulating current calculation: 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; and t is time. Calculation based on the data at t=10ms: Average current: (23.5A+24.2A+23.9A+24.7A) / 4 = 24.075A; Circulating current of module 1: Icirc,1 = 23.5A - 24.075A = -0.575A; Circulating current of module 2: Icirc,2 = 24.2A-24.075A = 0.125A; Circulating current of module 3: Icirc,3 = 23.9A - 24.075A = -0.175A; Circulating current of module 4: Icirc,4 = 24.7A - 24.075A = 0.625A.

[0117] 1.4. Compensated circulating current calculation: Icirc_comp,i(t)=Icirc,i(t)·(1+ηi(Vgs,i,Tj,i,Vdc)); where Icirc_comp,i(t) is the compensated circulating current; Icirc,i(t) is the original circulating 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 time; and i is the module index. Parameter difference compensation coefficient calculation: η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; circulating current after compensation: 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.

[0118] 1.5. The preprocessed voltage data Vp at t=10ms is: [750.2V, 742.8V, 748.1V, 745.3V]; the preprocessed state data Sp at t=10ms is: {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: 10ms}.

[0119] Step 2: Topology-aware phase space reconstruction.

[0120] 2.1. The compensation circulating current Icirc_comp within the 20ms window is segmented according to the switching state. The following typical switching 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].

[0121] 2.2. Extract the topological connection matrix Tconn based on the switch states: the topological matrix for state 1 [1,0,0,1] is [[1,0],[0,1],[0,1],[1,0]]; the topological matrix for state 2 [1,0,1,0] is [[1,0],[0,1],[1,0],[0,1]]. Calculate the equivalent circuit parameters: State 1 equivalent inductance Leq(S1) = 1.95 mH; State 1 equivalent resistance Req(S1) = 0.12 Ω; State 2 equivalent inductance Leq(S2) = 2.05 mH; State 2 equivalent resistance Req(S2) = 0.13 Ω. Analyze the circulating current frequency characteristics: State 1 main frequency: f_main(S1) = 270 Hz, bandwidth BW(S1) = 150 Hz; State 2 main frequency: f_main(S2) = 250 Hz, bandwidth BW(S2) = 140 Hz. Calculate the optimal time delay: Delay factors: α = 0.25, β = 0.5, reference inductance Leq_ref = 2.0 mH; State 1 time delay: τ1 = 0.25 × (1 / 270) × (1.95 / 2.0) 0.5 ≈ 0.9ms; State 2 time delay: τ2 = 0.25 × (1 / 250) × (2.05 / 2.0) 0.5≈ 1.0ms. Signal complexity evaluation: Complexity index for state 1: CI(S1) = 0.65; complexity index for state 2: CI(S2) = 0.71. Optimal embedding dimension calculation: Base 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. The reconstruction parameter set RP is formed as follows: State 1: {τ1=0.9ms,m1=6,Tconn1=[[1,0],[0,1],[0,1],[1,0]]}; State 2: {τ2=1.0ms,m2=6,Tconn2=[[1,0],[0,1],[1,0],[0,1]]}. The parameter sensitivity index PS is calculated as follows: State 1: PS(S1)=0.72; State 2: PS(S2)=0.85.

[0122] 2.3. Variable-length adaptive sliding window: Main period of circulating current: T_circ = 1 / 260 Hz ≈ 3.85 ms; Safety factor: β = 2.5; Minimum window length: Lmin = 5 ms; State 1 window length: Lwindow(S1) = max(5 ms, 2.5 × 3.85 ms) = 9.63 ms; State 2 window length: Lwindow(S2) = max(5 ms, 2.5 × 4.0 ms) = 10.0 ms; Key feature point extraction (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 the 10 ms window, accounting for approximately 11.5% of the total number of points.

[0123] 2.4. Basic delay embedding vector: Taking t = 10 ms, state 1, and module 1 as an example: τ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 characteristics: Module 1 output voltage: Vout,1 (10ms) = 750.2V; Voltage change rate: dV / dt (10ms) = (750.2V-745.8V) / 0.1ms = 44kV / s; Voltage unbalance: Vunbal (10ms) = 0.0045. Calculate circulating current dynamic characteristics: Change rate: dI / dt (10ms) = (-0.5756A+0.5623A) / 0.9ms = -14.8A / s; Acceleration: d 2 I / dt 2 (10ms) = -252A / s 2 Change trend: Itrend(10ms) = -1 (decreasing trend). Calculated temperature characteristics: Junction temperature difference: ΔTj(10ms) = 48.3°C - 44.8°C = 3.5°C; Temperature change rate: dTj / dt(10ms) = 0.042°C / s. Construct the weight fusion matrix (steady-state condition): W_fusion = [0.35, 0.15, 0.15, 0.15, 0.10, 0.10] (corresponding to delay sequence, voltage, current dynamics, temperature difference, temperature change rate, switching characteristics); generate the enhanced state space vector (module 1, t=10ms): delay sequence component: 0.35×[-0.5756A, -0.5623A, -0.5342A, -0.4981A, -0.4625A, -0.4283A]; voltage characteristic component: 0.15×[750.2V, 44kV / s, 0.0045]; current dynamic component: 0.15×[-14.8A / s, -252A / s 2 , -1]; temperature characteristic component: 0.10 × [3.5°C, 0.042°C / s]; switching characteristic component: 0.10 × [0.12μs, 0.25μs, 2.1V]; the complete enhanced state vector X(10ms) includes all of the above components. Normalization processing yields the standardized enhanced state vector Xnorm(10ms).

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

[0125] 2.6. Calculation of the inter-module coupling strength matrix: 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): Local linear prediction model construction: 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. Multi-step prediction sequence generation: t+1ms prediction: Icirc_pred,1(11ms) = -0.5832A; t+2ms prediction: Icirc_pred,1(12ms) = -0.5901A; t+5ms prediction: 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 confidence: Conf,1 = 0.92. Final dynamic trend prediction data: prediction sequence: DTP = [Icirc_pred(11ms),Icirc_pred(12ms),...,Icirc_pred(20ms)]; trend indicators: max(dI / dt) = -15.3A / s, ΔE = -0.12μJ, Stability = 0.85.

[0126] Step 3: Calculation of working condition adaptive virtual impedance parameters.

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

[0128] 3.2 Energy Flow Characteristic Data Calculation: Circulating Current Energy: Ecirc (10ms) = 3.25μJ; Power Flow: Pcirc (10ms) = -0.28W; Energy Flow Direction Index: Denergy (10ms) = -0.3853; Energy Flow Characteristic Data: EF (10ms) = [3.25μJ, -0.28W, -0.3853]. Parameter Sensitivity Data: Resistance Differential Sensitivity: SR (10ms) = 2.35A / Ω; Inductance Differential Sensitivity: SL (10ms) = 0.87A / mH; Capacitance Differential Sensitivity: SC (10ms) = 0.0023A / μF; Parameter Sensitivity Data: S (10ms) = [2.35A / Ω, 0.87A / mH, 0.0023A / μF]. Complexity characteristics of the circulating current signal: Sample entropy calculation (m = 2, r = 0.15 × standard deviation): Similar subsequence ratio: 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 signature: 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 values ​​of the 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 feature: CI(10ms) = 0.30×0.405 + 0.25×0.387 + 0.20×1.982 + 0.15×0.356 + 0.10×0.007=0.6688.

[0129] 3.3. Determination of Critical Control Moments: Predicted peak time t_peak = 17ms; Predicted zero crossing t_zero = 25ms (next cycle); Maximum rate of change t_max_rate = 13ms; Maximum energy t_max_energy = 16ms; Critical moment weights W_time = [0.35, 0.20, 0.25, 0.20]; Critical moment data KT = {t_peak: 17ms, t_zero: 25ms, t_max_rate: 13ms, t_max_energy: 16ms}. Predicted energy characteristics: Future energy flow direction Denergy_pred (15ms) = -0.412; Predicted energy peak Epeak_pred = 3.42 μJ (t = 16ms); Energy change sensitivity E_sensitivity = 2.35; Predicted energy characteristics PE = [Denergy_pred, Epeak_pred, E_sensitivity]. Prediction parameter sensitivities: Resistance sensitivity prediction SR_pred(15ms) = 2.42A / Ω; Inductance sensitivity prediction SL_pred(15ms) = 0.91A / mH; Capacitance sensitivity prediction SC_pred(15ms) = 0.0025A / μF; Prediction sensitivity PS = [SR_pred, SL_pred, SC_pred]. Prediction complexity characteristics: Prediction entropy change CI_pred(15ms) = 0.683; Prediction state transition point T_change_pred = 22ms; Prediction complexity characteristics PC = [CI_pred, T_change_pred]. Impedance-circulating current mapping model: A mapping relationship, Z_I_model, is established based on historical data. For state 1, the relationship between the peak circulating current and the virtual impedance is: Ipeak = 0.98 × Ipeak_0 × exp(-0.27 × Z). For module 1, the predicted relationship at t = 15 ms is: Icirc_pred,1(15 ms) = -0.6023 A × exp(-0.27 × Z). Multi-objective optimization function design: Circulating current suppression objective: J_circ = Σ|Icirc_pred,i| 2 = 0.6023 2 + 0.1323 2 + 0.1853 2 + 0.6353 2 =0.7932; Dynamic response target: J_resp = Σ|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 based on working condition indicator 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 (next 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Ω.

[0130] 3.4. CFPS extracts key features of circulating current phase characteristics: trajectory geometry: GF(10ms) = [25.7, 158.3]; stability index: SI(10ms) = 0.87; dynamic evolution characteristics: DF(10ms) = [294.8, 4327.6, -1215.4]. Adaptive PID controller design: KP is adjusted based on trajectory curvature: KP = 0.5 × (1 + 0.2 × tanh(25.7-20)) = 0.6; KI is adjusted based on the stability index: KI = 0.1 × (1 + 0.5 × (0.87-0.5)) = 0.1185; KD is adjusted 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 characteristics and period characteristics: Inductance sensitivity threshold: SL_threshold = 0.8A / mH; SL (10ms) = 0.87A / mH > 0.8A / mH, adjustment 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 cyclic mode CP = 2: KD_adj = KD × 1.05 = 0.065 × 1.05 = 0.0683; Feature matching parameter: FMP = [1.0131, 0.964, 1.05]. Basic impedance value calculation: Circulating current of module 1 at t=10ms: Icirc_comp,1(10ms) = -0.5756A; Integral term (assuming the past 10ms window): ∫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 state variable subset: 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Ω; Absolute value: Zv_state,1 = 0.6919Ω. Resonance suppression compensator design: Identify the resonant frequency: f_res = 270Hz; Resonance suppressor parameters: A_res = 0.18Ω, φ_res = 180°; Suppression component at t = 10ms: Zv_res,1 = 0.18×sin(2π×270×0.01+π)=0.1712Ω. Stability enhancer design: Stability threshold: λL_threshold = 0.05; λL(10ms) = -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Ω.

[0131] 3.5 Physical Limitation Constraints 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Ω.

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

[0133] 4.1 Definition of System State Vector: 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 . Settings of Lyapunov Function Parameters: 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]]. Calculation of Lyapunov Function: 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. Time derivative of the Lyapunov function: estimated rate of change of circulating current: 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; Lyapunov function time derivative: 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.00474 J, dV / dt=0.03593 W, SI=-0.03593}.

[0134] 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. Lyapunov stability-guided weight adjustment: Module 1 parameters (10ms): Icirc_comp,1=-0.5756A; Zv_fb,1 = 0.86Ω; Zv_ff,1 = 0.46Ω (current time); |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; exponent 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 . Calculate dω / dt of other modules similarly, and adjust the weight vector: dω / dt = [0.07026, -0.00582, 0.01621, -0.07965] T s -1Output basic 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; Basic mixing weights: ω(11ms)=[0.70070,0.69994,0.70016,0.69920] T .

[0135] 4.3. Circulating current frequency 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; intermediate frequency layer: V_mf(10ms) = 4.421×10 -4 J; low frequency layer: V_lf(10ms) = 1.611×10 -6 J; Frequency 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 IF weight change: dω_mf,1 / dt = 0.0582 s -1 ; Low frequency weight change: dω_lf,1 / dt = 0.0103 s -1The optimal high-frequency weight is: ω_hf,1(11ms) = 0.7 + 0.1536 × 0.001 = 0.7154; the optimal mid-frequency weight is: ω_mf,1(11ms) = 0.7 + 0.0582 × 0.001 = 0.7058; the optimal low-frequency weight is: ω_lf,1(11ms) = 0.7 + 0.0103 × 0.001 = 0.7010. Dynamic band weight calculation: Based on the band energy distribution, the high-frequency energy proportion is: E_hf / E_total = 0.152; the mid-frequency energy proportion is: E_mf / E_total = 0.831; the low-frequency energy proportion is: E_lf / E_total = 0.017. Band weights: λ_hf(10ms) = 0.152; λ_mf(10ms) = 0.831; λ_lf(10ms) = 0.017. Multi-scale mixing 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; multi-scale mixing weights of other modules are calculated similarly; final multi-scale mixing 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, the stability condition is met.

[0136] 4.4. Read the feedforward and feedback virtual impedance parameters: Module 1 feedback virtual impedance: Zv_fb,1(11ms) = 0.865Ω; Module 1 feedforward virtual impedance: Zv_ff,1(11ms) = 0.46Ω; Module 1 multi-scale mixing weight: ω_ms,1(11ms) = 0.7071. Calculate the comprehensive virtual impedance of module 1: Zv,1(11ms) = 0.7071 × 0.865 + (1 - 0.7071) × 0.46; Zv,1(11ms) = 0.6116 + 0.1347 = 0.7463Ω. The integrated virtual impedance of other modules is calculated similarly: 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Ω. To verify the stability of the fused impedance, calculate the rate of change of the Lyapunov function: dV / dt|Zv = -1.237×10 -3 W<0, meeting the stability condition. The final comprehensive virtual impedance value: Zv(11ms) = [0.7463, 0.637, 0.6781, 0.7157] T Ω.

[0137] 4.5. Calculation of module voltage correction values: 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 balancing 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 balancing compensation: ΔV1 = 0.05 × (P* - P1) / P* = 0.05 × (17.975 - 17.63) / 17.975 = 0.00096V; balancing compensation for other modules is calculated similarly; 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] TModule 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: 1000 Hz; DC bus voltage: Vdc = [3050 V, 3030 V, 3070 V, 3020 V]; module 1 modulation ratio: m1 = V1_cmd / Vdc,1 = 751.7341 / 3050 = 0.2465; calculate the duty cycle and phase of the PWM wave for each module to generate the final PWM control signal.

[0138] To verify the effectiveness of this embodiment, comparative tests were conducted with traditional control methods. These include: the fixed virtual impedance method (FVI), which uses a fixed virtual impedance value for control; the single PID feedback method (SPID), which uses only PID-based feedback control; and the conventional feedforward-feedback hybrid method (CFFB), which uses a fixed-weight fusion of feedforward and feedback control. The tests were conducted on the same hardware platform and operating conditions, and the following indicators were compared and analyzed: The peak circulating current of this embodiment was reduced to 12.5% ​​of the original value, and the RMS value was reduced by 87.3%; the peak circulating current of the FVI method was reduced to 43.2% of the original value, and the RMS value was reduced by 59.7%; the peak circulating current of the SPID method was reduced to 35.8% of the original value, and the RMS value was reduced by 64.2%; and the peak circulating current of the CFFB method was reduced to 21.3% of the original value, and the RMS value was reduced by 79.5%. The dynamic response performance (settling time after a step load change) of this embodiment was 4.2ms, compared to 15.7ms for the FVI method, 9.3ms for the SPID method, and 7.8ms for the CFFB method. The grid disturbance adaptability (circulating current fluctuation during a 15% voltage sag) of this embodiment is 7.3%, compared to 28.5% for the FVI method, 19.2% for the SPID method, and 14.7% for the CFFB method. The parameter robustness of this embodiment (system stability under conditions of ±15% module parameter deviation) is fully stable, with a 10.2% increase in circulating current. The FVI method is unstable when parameter deviations exceed ±12%, the SPID method increases circulating current by 42.5% when parameter deviations exceed ±15%, and the CFFB method increases circulating current by 27.8% when parameter deviations exceed ±15%.

[0139] The present invention achieves precise matching of phase space reconstruction parameters with the real-time operating state of the power electronic system. This enables the phase space to more accurately capture the dynamic characteristics of circulating currents, particularly in scenarios where the switching state changes rapidly. For high-voltage SVG systems with high-frequency switching, it can adaptively respond to changes in system dynamics under different switching combinations, fundamentally improving the effectiveness of subsequent control strategies. In particular, in high-power modular systems, the circulating current suppression efficiency is increased by more than three times. By constructing a multi-physics enhanced state space that integrates electrical, thermal, and switching characteristics, the limitations of the traditional single current state space are fundamentally overcome, and the coupling relationship between different physical quantities in the system is fully captured. In the actual operation of high-voltage SVGs, it can simultaneously reflect the impact of power device junction temperature differences on circulating currents and the system asymmetry caused by voltage imbalance. This enables the control system to foresee and respond to changes in circulating currents under the combined action of multiple physical factors, particularly in high-voltage and high-power scenarios, improving the system stability margin. By calculating the geometric characteristics of the phase trajectory, a complete set of circulating current phase characteristics (CFPS) is established, which directly reflects the intrinsic characteristics of the circulating current dynamic system. Unlike traditional methods that focus solely on surface features of amplitude and spectrum, this method deeply analyzes the geometry and stability of phase space trajectories, enabling identification of critical stable states and nonlinear dynamic behavior in circulating currents. In practical applications, this system can identify precursor signals before circulating current instability occurs. Especially in load-sudden change scenarios, preventive control can be initiated earlier than traditional methods, suppressing circulating current peaks below safety thresholds and improving system reliability and robustness. By simultaneously calculating multidimensional entropy metrics, a comprehensive signal complexity analysis framework is established, enabling precise quantification of the non-stationary dynamic characteristics of circulating currents. This method surpasses traditional statistical feature analysis and can reveal the inherent complexity and predictability of circulating currents from an information-theoretic perspective. Particularly under conditions of grid harmonics and rapid load fluctuations, the complexity characteristics derived from entropy analysis can accurately identify system behavior patterns, predict drastic changes in circulating currents in advance, and provide sufficient response time for the control system. This directly improves the targeted optimization of virtual impedance parameters and enhances the safety margin of system operation. By constructing a stability-oriented gain modulation law, stable control with rigorous mathematical proof is achieved. This theoretically guarantees the stable convergence of the system, addressing the fundamental lack of stability assurance inherent in traditional empirical parameter tuning methods. In actual operation of high-voltage SVG systems, the derivative of the Lyapunov function is consistently negative under all operating conditions, theoretically guaranteeing system stability. This approach is particularly applicable to grid fluctuations in high-voltage transmission systems, improving the reliability and robustness of SVG systems. By decomposing circulating current into high-, medium-, and low-frequency components and designing independent control layers for each frequency band, while dynamically adjusting control weights based on the energy distribution of each frequency band, precise control is achieved for disturbances of varying timescales.It overcomes the limitation of traditional single controllers that are difficult to deal with fast switching transients and slow parameter drifts at the same time, and enables the control system to process circulating current components of different frequencies in a targeted manner. In practical applications, the high-frequency control layer can quickly suppress the circulating current spikes caused by switching transients, the medium-frequency control layer handles the fluctuations caused by load changes, and the low-frequency control layer compensates for long-term drifts such as device aging. Test results show that compared with single-band control, the present invention improves efficiency while maintaining effective compensation for low-frequency drifts; it can maintain optimal performance under disturbances of various time scales, and comprehensively improves the dynamic response capability and steady-state accuracy of the SVG system. The present invention effectively solves the circulating current problem caused by parameter differences between modules, and improves the reliability, stability and power quality of system operation.

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

Claims

1. A modular high-voltage SVG circulating current suppression method, characterized in that: include: Acquire module data, pre-process it, and obtain pre-processed data, including circulating current data, voltage data, and switch status data; Constructing an enhanced phase space based on preprocessed data to generate characteristic data, including circulating current phase characteristics and dynamic trend prediction data; Based on the characteristic data and combined with the physical constraints of the system, the virtual impedance parameters are calculated, including feedforward and feedback virtual impedance parameters; Use virtual impedance parameters to perform gain modulation, generate control signals, and output control data, including module control voltage correction values ​​and PWM control signals; Generate feature data, including: Analyze the state data and divide the circulating current data into predetermined subsequences according to the switch state to form segmented circulating current data; Read the system topology and switch states, dynamically calculate the time delay and embedding dimension for the subsequence, and generate reconstruction parameters; Key feature points are extracted from circulating current data and combined with reconstruction parameters to construct an enhanced state space. The time delay sequence in the circulating current data is fused with physical quantities to form an enhanced state vector. Based on this, the geometric characteristics of the phase trajectory are calculated to generate the circulating current phase feature. Evaluate the coupling relationship of circulating currents between modules to form a module coupling matrix, and combine the phase characteristics of circulating currents to generate dynamic trend prediction data; Calculate virtual impedance parameters, including: Analyze the energy distribution and flow characteristics of circulating current data within the switching cycle to generate energy flow characteristic data; Evaluate the sensitivity of circulating current data to module parameter differences 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 this with 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. An adaptive feedback control strategy is constructed based on the phase characteristics of the circulating current, the real-time response virtual impedance value is calculated, and the feedback virtual impedance parameters are generated.

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

3. The method according to claim 1, characterized in that Form an enhanced state vector, including: Constructing a basic time delay embedding vector based on circulating current data and reconstruction parameters; Extract key voltage features from voltage data; Calculate the dynamic characteristics of circulating current data; Calculate the temperature characteristics of key power devices, including junction temperature difference and temperature change rate, based on pre-stored system temperature data; The basic time delay embedding vector is weightedly fused with key voltage characteristics, dynamic characteristics and temperature characteristics through a multi-physical quantity weight fusion matrix to form an enhanced state vector.

4. The method according to claim 1, wherein Generates circulating current phase signatures, including: Read and, based on the enhanced state vector, sequentially: calculate the trajectory curvature and torsion to form the trajectory geometric characteristics; evaluate the local stability of the phase trajectory to form a stability index; calculate the phase space velocity vector and acceleration vector to form the dynamic evolution characteristics; Based on the circulating current data and state data, energy flow characteristics are formed; Extract the sensitivity of circulating current data to module parameter differences to form parameter sensitivity characteristics; The weights of trajectory geometric characteristics, stability index, dynamic evolution characteristics, energy flow characteristics and parameter sensitivity characteristics are determined to form the circulating current phase characteristics.

5. The method according to claim 1, wherein Generate dynamic trend forecast data, including: Based on the enhanced state vector and circulating current phase characteristics, a point set similar to the current state is selected from the historical data to build 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 the prediction result is used as the new initial condition to generate a multi-step prediction sequence; Circulating current prediction data are extracted from the multi-step prediction sequence, and the key dynamic indicators of the predicted circulating current are calculated and integrated into dynamic trend prediction data.

6. The method according to claim 1, characterized in that Analyze the complexity characteristics of circulating current, including: To read the circulating current data, perform the following operations: Divide it into subsequences of length m and m+1, calculate the proportion of similar subsequences, determine the conditional probability, and obtain the original signal sample entropy; where m is a natural number greater than 0; Perform coarse-graining processing on it, calculate the sample entropy of each coarse-grained sequence, form an entropy-scale curve, and generate multi-scale entropy features; Analyze sequence patterns and fuzzy boundaries, and calculate permutation entropy and fuzzy entropy; The change rates of the original signal sample entropy, multi-scale entropy features, permutation entropy and fuzzy entropy are calculated and weighted fused accordingly to generate the circulating current complexity features.

7. The method according to claim 1, characterized in that Form feedforward virtual impedance parameters, including: Predict the dynamic characteristics of circulating current based on dynamic trend prediction data analysis and determine the critical control moments; Based on the energy flow characteristic data, parameter sensitivity data and circulating current complexity characteristics, the change amount of energy flow characteristics, parameter sensitivity and circulating current complexity at a future moment is predicted to form prediction characteristic data; For the determined critical control moments, the predicted characteristic data and the pre-built impedance optimization model are applied to solve the optimal virtual impedance value within the predicted time window, generate the virtual impedance time series, and form the feedforward virtual impedance parameters.

8. The method according to claim 1, characterized in that Generate feedback virtual impedance parameters, including: Read the circulating current phase characteristics and: Extract trajectory geometric features, stability indicators and dynamic evolution characteristics from them and construct PID control parameters; Parameter sensitivity characteristics and periodic characteristics are extracted from them, PID control parameters are adjusted and stability indicators are evaluated to form characteristic matching parameters, resonance suppression components and stability enhancement components; basic impedance values ​​are generated through PID control parameters and characteristic matching parameters; Dynamically adjust a preconfigured state feedback gain matrix according to the state feedback gain matrix to form a state feedback component; The basic impedance value, state feedback component, resonance suppression component and stability enhancement component are weighted and fused to form the feedback virtual impedance parameter.

Images