Adaptive coordination method for dynamic droop parameters of heterogeneous battery modules

Abstract

The present invention provides a method for adaptively coordinating dynamic droop parameters of heterogeneous battery modules, including: real-time identification of transfer function parameters of each heterogeneous module, using virtual pole supplementation technology to solve the pole alignment problem of different order systems; decomposing the multi-objective optimization problem into three independent sub-problems of low-frequency band SOC balance, mid-frequency band power regulation and high-frequency band transient suppression based on modal decoupling theory; realizing multi-band parallel gradient projection optimization through pre-calculated constrained feasible domain and versioned synchronization mechanism, effectively avoiding thread competition; adaptively adjusting fusion weights according to the convergence quality of each frequency band to generate coordinated droop parameters. The present invention fundamentally solves the problem of mismatched frequency domain responses of heterogeneous battery modules, achieves a technological breakthrough from static pole configuration to dynamic adaptive adjustment, and improves the stability and efficiency of heterogeneous energy storage systems.

Description

Technical Field

[0001] The present invention belongs to the field of control systems, and in particular to a method for adaptively coordinating dynamic droop parameters of heterogeneous battery modules. Background Art

[0002] With the rapid development of new energy technologies, heterogeneous battery energy storage systems have attracted widespread attention due to their ability to fully utilize the performance advantages of different types of batteries. In practical applications, battery modules with different chemical systems, such as lithium iron phosphate batteries and ternary lithium batteries, often need to be operated in parallel to achieve the complementary advantages of high energy density and long cycle life. However, due to the significant differences in electrochemical properties, internal resistance characteristics, and dynamic response characteristics, heterogeneous battery modules face severe challenges such as uneven power distribution, large circulating current losses, and poor system stability when operated in parallel. Droop control, as a distributed control strategy that does not require communication, has important application value in heterogeneous battery systems. However, how to achieve dynamic optimization of droop parameters so that battery modules with different characteristics can operate in coordination has become a key technical bottleneck restricting the large-scale application of heterogeneous energy storage systems.

[0003] At present, research on droop control of battery energy storage systems mainly focuses on homogeneous battery systems or static parameter configuration. Existing technologies usually adopt a droop parameter design method based on a steady-state model, determine a fixed droop coefficient through offline calculation, and keep it unchanged during system operation. Some studies have considered the impact of battery SOC differences on droop parameters and proposed a piecewise linear adjustment strategy based on SOC. In terms of heterogeneous battery systems, existing methods mostly adopt a unified control parameter design criterion, which does not fully consider the dynamic characteristics differences of different types of batteries. Some studies have attempted to optimize droop parameters through time-domain simulation, but the computational complexity is high and it is difficult to implement online applications. In terms of multi-module coordinated control, existing technologies mainly use centralized optimization algorithms to obtain the control parameters of each module by solving large-scale optimization problems, but there are problems with heavy computational burden and poor scalability.

[0004] However, existing technologies have key technical defects when dealing with the dynamic coordination of heterogeneous battery modules. First, in terms of frequency domain characteristic matching, existing methods fail to effectively solve the problem of coupling the response characteristics of batteries of different chemical systems in multiple frequency bands. When the system simultaneously has multi-time scale control requirements such as SOC balance, power regulation and transient suppression, unified parameter optimization often leads to performance deterioration in certain frequency bands and even causes cross-band oscillations. Secondly, in terms of adaptability to dynamic working conditions, the existing static pole configuration method cannot track the evolution of dynamic characteristics caused by factors such as battery aging, temperature changes and SOC drift, resulting in the preset control parameters gradually deviating from the optimal operating point and the continuous deterioration of system performance. Summary of the Invention

[0005] The purpose of the invention is to provide a method for adaptively coordinating dynamic droop parameters of heterogeneous battery modules, in order to solve at least one technical problem existing in the prior art.

[0006] The technical solution is an adaptive coordination method for dynamic droop parameters of heterogeneous battery modules, including:

[0007] Acquire real-time dynamic response data of heterogeneous battery modules and obtain a standardized module state matrix through normalization. Based on this, extract transfer function parameters and calculate the pole distribution matrix.

[0008] Based on the pole distribution matrix, the frequency domain characteristics are hierarchically decoupled through the decoupling transformation matrix to obtain a hierarchical optimization objective function set; based on this, multi-band synchronous solution is performed to obtain the frequency band matching parameter matrix;

[0009] Based on the frequency band matching parameter matrix, the droop parameter adjustment instruction is obtained through weighted fusion processing.

[0010] Beneficial effect: The present invention fundamentally solves the problem of mismatched frequency domain responses of heterogeneous battery modules, achieves a technological breakthrough from static pole configuration to dynamic adaptive adjustment, and improves the stability and efficiency of heterogeneous energy storage systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 A flowchart of the steps of a method for adaptively coordinating dynamic droop parameters of a heterogeneous battery module provided in an embodiment of the present application.

[0012] Figure 2 A flowchart of the steps for obtaining a hierarchical optimization objective function set provided in an embodiment of the present application.

[0013] Figure 3 A flowchart of the steps for obtaining a pole distribution matrix provided in an embodiment of the present application.

[0014] Figure 4 A flowchart of the steps for obtaining a frequency band matching parameter matrix provided in an embodiment of the present application. DETAILED DESCRIPTION

[0015] 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.

[0016] 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.

[0017] The research found that when the transfer functions of heterogeneous modules have inconsistent orders, traditional methods lack effective pole alignment and coordination mechanisms, which can easily lead to undesirable resonant peaks at specific frequencies. Furthermore, existing centralized optimization algorithms experience exponential computational complexity when dealing with large-scale heterogeneous systems and lack effective parallelization mechanisms, making them difficult to meet the demands of online real-time optimization.

[0018] like Figure 1 As shown in FIG, a dynamic droop parameter adaptive coordination method for heterogeneous battery modules is proposed, which includes the following steps:

[0019] S1. Obtain real-time dynamic response data of heterogeneous battery modules and obtain a standardized module state matrix through normalization processing;

[0020] S2. Based on the standardized module state matrix, the transfer function parameters are extracted and the pole distribution matrix is ​​calculated;

[0021] S3. Based on the pole distribution matrix, the frequency domain characteristics are hierarchically decoupled through the decoupling transformation matrix, and the optimization objective is decomposed into three independent sub-problems of low frequency band, medium frequency band and high frequency band, and a hierarchical optimization objective function set is obtained;

[0022] S4, performing multi-band synchronous solution based on the hierarchical optimization objective function set to obtain the band matching parameter matrix;

[0023] S5. Based on the frequency band matching parameter matrix, a droop parameter adjustment instruction is obtained through weighted fusion processing.

[0024] Specifically, real-time dynamic response data, including real-time voltage and current data, SOC (remaining battery charge) values, and temperature parameters, are processed through multi-channel synchronous sampling to obtain a standardized module state matrix. Based on a hierarchical optimization objective function set, a parallel gradient projection algorithm is used to perform a multi-band synchronous solution to obtain a band matching parameter matrix.

[0025] According to one aspect of the present application, obtaining a standardized module state matrix includes:

[0026] S11. Automatically identify the chemical type of each battery module based on real-time dynamic response data, and distinguish between lithium iron phosphate modules and ternary lithium modules;

[0027] S12. According to the chemical type identification results, the matching normalization coefficients of the lithium iron phosphate module and the ternary lithium module are respectively calculated;

[0028] S13. Based on the matched normalization coefficients, the module data of different chemical types are subjected to differentiated normalization processing to obtain a normalized module state matrix.

[0029] Specifically, the raw voltage, current, and temperature signals are read and sampled using a 16-bit ADC via a hardware synchronization trigger mechanism to generate a synchronized raw data set. This synchronized raw data set is then processed using a module type identification and parameter normalization algorithm to generate a standardized module state matrix. Different chemical types of modules use their own normalization coefficients: the lithium iron phosphate module's normalization coefficient is α_LFP, and the ternary lithium module's normalization coefficient is α_NCM.

[0030] like Figure 3 As shown, according to one aspect of the present application, obtaining the pole distribution matrix includes:

[0031] S21. Based on the standardized module state matrix, identify the chemical type of each heterogeneous battery module and set the corresponding characteristic coefficient; set the characteristic coefficient γ_LFP = 0.8 for the lithium iron phosphate module and set the characteristic coefficient γ_NCM = 1.2 for the ternary lithium module;

[0032] S22. Calculate heterogeneous property weighted prediction errors based on the chemical type and property coefficients, where modules of different chemical types use corresponding property coefficients for error weighting processing;

[0033] S23, based on the weighted prediction error of heterogeneous characteristics, dynamically adjust the forgetting factor value through the adaptive forgetting factor adjustment mechanism, calculate the adaptive forgetting factor, and obtain a dynamic forgetting factor sequence;

[0034] S24. Based on the dynamic forgetting factor sequence, perform parameter iteration to obtain a dynamic transfer function parameter set and calculate the pole distribution matrix accordingly.

[0035] According to one aspect of the present application, before extracting the transfer function parameters, the normalized module state matrix is ​​subjected to sliding window data preprocessing, specifically:

[0036] Based on the standardized module state matrix, the real-time change rate dY / dt of the signal is calculated to identify the severity of the dynamic change of the signal;

[0037] Based on the real-time change rate, the sliding window length is dynamically adjusted according to the relationship between the window lengths;

[0038] Based on the dynamically adjusted sliding window length, the standardized module state matrix is ​​segmented to obtain time domain response data segments.

[0039] Specifically, the standardized module state matrix is ​​read, and the data is segmented by the adaptive window length algorithm to obtain the time domain response data segment. The window length is dynamically adjusted according to the signal change rate: L_window=L_base×(1+β×|dY / dt|), where L_base is the basic window length of 100ms and β is the adjustment coefficient. The time domain response data segment is read, and the parameters are identified by the recursive least squares algorithm to obtain the dynamic transfer function parameter set. The forgetting factor adaptive adjustment mechanism Ξ(k)=Ξ0+(1-Ξ0)×exp(-σ×|e(k)|) is adopted, where Ξ0 is the initial forgetting factor, e(k) is the prediction error, and σ is the sensitivity coefficient, so that the algorithm can quickly adapt to the parameter changes of heterogeneous modules. The dynamic transfer function parameter set is read, and the poles of the transfer function of each module are calculated by the characteristic root solution algorithm to obtain the pole distribution matrix. For the second-order system G(s)=K / (s 2 +2ζω n s+ω n 2 ), the pole is calculated as Ξ1,2=-ζω n ±ω n sqrt(ζ 2 -1), and record the real and imaginary information of the poles at the same time. Where G(s) is the transfer function, which shows how the input signal affects the output signal; K is the system gain, s is the Laplace variable in the complex domain, ζ is the damping ratio, ω n is the natural frequency.

[0040] In this embodiment, the weight of the weighted prediction error of heterogeneous characteristics is determined as follows: LFP battery: internal resistance reference R0 = 12mΩ, capacity attenuation rate δ = 0.8% / year; NCM battery: internal resistance reference R0 = 8mΩ, capacity attenuation rate δ = 1.2% / year; γ_LFP = (R0_NCM / R0_LFP) × (δ_LFP / δ_NCM) = (8 / 12) × (0.8 / 1.2) = 0.67 × 0.67 ≈ 0.8; γ_NCM = (R0_LFP / R0_NCM) × (δ_NCM / δ_LFP) = (12 / 8) × (1.2 / 0.8) = 1.5 × 1.5 = 1.2; prediction error e_LFP = |SOC_measured - SOC_predicted| = 0.03; weighted error e_weighted_LFP = 0.03 × 0.8 = 0.024. Where SOC_measured is the actual measured battery state, and SOC_predicted is the predicted battery state.

[0041] According to one aspect of the present application, setting the characteristic coefficient includes:

[0042] Based on a standardized module state matrix, the internal resistance characteristics and capacity decay patterns of battery modules of different chemical types are analyzed. Among them, lithium iron phosphate modules have low internal resistance stability, while ternary lithium modules have high energy density and fast response characteristics.

[0043] Based on the differences in internal resistance characteristics and capacity attenuation laws, a mapping relationship between chemical type and characteristic coefficient is established, and different characteristic coefficients are set for battery modules of different chemical types;

[0044] By comparing the parameter identification accuracy under different characteristic coefficients, the rationality of the characteristic coefficient setting is verified to ensure that the identification performance differences between heterogeneous battery modules are controlled within the preset range.

[0045] Specifically, prediction error calculation and heterogeneous characteristic identification are performed. The time domain response data segment and the pre-stored heterogeneous module feature library are read, and the prediction error of each module is calculated using the heterogeneous characteristic matching algorithm to obtain the heterogeneous characteristic weighted error vector. The specific processing process is as follows: the chemical type identifier of each module is identified, the characteristic coefficient γ_LFP is set to 0.8 for the lithium iron phosphate module, and the characteristic coefficient γ_NCM is set to 1.2 for the ternary lithium module; the weighted prediction error e_weighted(k) is calculated as γ_type×|y_actual(k)-y_predict(k)|, where y_actual is the actual output, y_predict is the model predicted output, and γ_type is the characteristic coefficient of different types.

[0046] Perform dynamic adjustment of the adaptive forgetting factor. Read the weighted error vector of heterogeneous characteristics, calculate the forgetting factor through the error sensitivity analysis algorithm, and obtain the dynamic forgetting factor sequence. The core process is: establish the error sensitivity function S(k)=|e_weighted(k)| / (μ+Σ j=1 k |e_weighted(j)| / k), where k is the time step or sequence index; μ = 0.01 is the regularization coefficient; the adaptive forgetting factor Ξ_adaptive(k) = Ξ_min + (Ξ_max - Ξ_min) × exp(-α × S(k)) is calculated, where the minimum forgetting factor Ξ_min = 0.95, the maximum forgetting factor Ξ_max = 0.995, and α = 2.0 are adjustment parameters; when S(k) > 0.1, Ξ_adaptive(k) is forced to be set to Ξ_min to quickly adapt to parameter changes.

[0047] Perform recursive parameter updates and convergence checks. Read the dynamic forgetting factor sequence Ξ_adaptive and the parameter estimate θ_prev at the previous moment, perform parameter iteration using the recursive least squares update algorithm, and obtain the updated transfer function parameters θ_updated. Detailed processing flow: Calculate the gain matrix K(k)=P(k-1)φ(k) / [Ξ_adaptive(k)+φ T (k)P(k-1)φ(k)], where φ(k) is the regression vector; update parameter estimate θ_updated(k)=θ_prev(k-1)+K(k)[y(k)-φ T (k)θ_prev(k-1)]; Update the covariance matrix P(k)=[IK(k)φ T (k)]P(k-1) / Ξ_adaptive(k); perform convergence test ||θ_updated(k)-θ_prev(k-1)||2≤ε_conv, where the convergence threshold ε_conv=1e -6 , T is the transpose, I is the identity matrix, and y(k) is the observed data.

[0048] Format the transfer function parameters output. Read the updated transfer function parameters θ_updated, and obtain the dynamic transfer function parameter set G_params through parameter validity verification and standardized format conversion. The processing process includes: verifying the physical rationality of the parameters, requiring the denominator polynomial coefficients to be positive and meet the stability conditions; converting the parameters θ_updated into the standard transfer function format G(s)=(b_ns n +...+b_1s+b_0) / (a_ms m +...+a_1s+a_0); the output format is G_params = {[b_coeff], [a_coeff], module_id, timestamp}. Where b_n is the highest-order coefficient of the numerator polynomial, indicating the weight of the transfer function input affecting the output; a_m is the highest-order coefficient of the denominator polynomial, which determines the dynamic characteristics and stability of the system; b_coeff is the set of all coefficients of the numerator polynomial, used to describe the input response of the system; a_coeff is the set of all coefficients of the denominator polynomial, used to define the pole distribution and stability conditions of the system; module_id is the module identifier, and timestamp is the timestamp.

[0049] like Figure 2 As shown, according to one aspect of the present application, obtaining a hierarchical optimization objective function set includes:

[0050] S31. Based on the pole distribution matrix, determine the frequency band attribution characteristics of each pole, classify the poles into low frequency band, medium frequency band and high frequency band according to the natural frequency, and obtain a pole frequency band attribution table;

[0051] S32, based on the extreme frequency band attribution table, analyzing the coupling terms of the preset original unified objective function, and identifying the cross-band correlation matrix C_cross;

[0052] S33. Based on the cross-band correlation matrix, a decoupling transformation matrix T_decouple=I-α×C_cross is constructed, and the decoupling transformation is applied to decompose the original unified objective function into a predetermined number of independent sub-objective functions, where I is the identity matrix and α is the decoupling strength coefficient;

[0053] S34. Based on the sub-objective function, verify the decoupling effect, require that the non-diagonal elements of the correlation coefficient matrix after decoupling are less than 0.05, and construct a hierarchical optimization objective function set.

[0054] Specifically, the pole distribution matrix is ​​read, and layered processing is performed through frequency domain characteristic analysis to obtain a three-layer frequency domain feature set. The traditional full-frequency domain optimization is decomposed into three independent frequency bands: the low-frequency band F_low (0.1-1Hz) mainly affects the SOC balance, the mid-frequency band F_mid (1-5Hz) mainly affects the power regulation, and the high-frequency band F_high (5-10Hz) mainly affects the transient response. Each frequency band uses a different performance indicator weight. Read the three-layer frequency domain feature set, reconstruct the objective function through the decoupling mapping algorithm, and obtain a hierarchical optimization objective function set. The complex multi-objective optimization problem J_total=Σ i |Ξ i -Ξ ref | 2 Decompose into three independent sub-problems: J_low=Σ i ∈F_loww_low×|Ξ i -Ξ ref_low | 2 ; J_mid=Σ i ∈F_midw_mid×|Ξ i -Ξ ref_mid | 2 ; J_high = Σ i ∈F_highw_high×|Ξ i -Ξ ref_high | 2 ;Where w_low, w_mid, w_high are the frequency band weight coefficients. i is the i-th extreme point; ref is the reference pole; J_low is the low-frequency optimization objective function, ref_lowis the reference pole of the low frequency band, J_mid is the optimization objective function of the mid-frequency band, ref_mid is the reference pole of the mid-frequency band, J_high is the optimization objective function of the high-frequency band, ref_high The hierarchical optimization objective function set is read, and a hierarchical constraint set is obtained through physical constraint analysis and feasible domain mapping. Unlike traditional unified constraint processing, this embodiment sets dedicated constraints for each frequency band: the key constraint for the low-frequency band is steady-state accuracy Δ_ss ≤ 2%, the key constraint for the mid-frequency band is damping ratio ζ ≥ 0.6, and the key constraint for the high-frequency band is response time t_r ≤ 50ms.

[0055] According to one aspect of the present application, constructing a cross-band correlation matrix includes:

[0056] Based on the pole frequency band attribution table, the frequency domain coupling strength between poles in different frequency bands is calculated, and the correlation between frequency bands is determined by analyzing the pole distribution density and the frequency domain overlapping area.

[0057] The degree of correlation is quantified as the matrix element c_ij of the cross-band correlation matrix C_cross, where c_ij represents the coupling strength between the i-th frequency band and the j-th frequency band;

[0058] Apply the decoupling transformation matrix to transform the original correlation, calculate the decoupled correlation coefficient matrix, and verify whether the non-diagonal elements are less than 0.05. If the verification fails, adjust the decoupling strength coefficient α and recalculate.

[0059] Specifically, the pole distribution matrix Ξ_matrix is ​​read and classified using pole frequency calculation and frequency band partitioning algorithms to obtain the pole frequency band assignment table Ξ_classification. The classification process: Determine the frequency band assignment characteristics of each pole: ω_n = |Ξ| and damping frequency ω_d = |Im(Ξ)|; Establish a frequency band determination rule: When ω_d∈[0.628, 6.28]rad / s, it belongs to the low frequency band F_low; when ω_d∈[6.28, 31.4]rad / s, it belongs to the mid frequency band F_mid; when ω_d∈[31.4, 62.8]rad / s, it belongs to the high frequency band F_high; Classify the complex conjugate pole pairs as a whole to avoid pole pair separation; Output Ξ_classification = {Ξ_low[], Ξ_mid[], Ξ_high[], classification_matrix}. Among them, Ξ_low[] is the pole set of the low frequency band, Ξ_mid[] is the pole set of the mid frequency band, Ξ_high[] is the pole set of the high frequency band, and classification_matrix is ​​the pole classification matrix. Read the pole band attribution table and the pre-stored heterogeneous module performance index Performance_targets, calculate the weight coefficient through the band importance evaluation algorithm, and obtain the adaptive band weight matrix W_adaptive. Different from the traditional fixed weight algorithm, evaluate the pole distribution density of each band ρ_i=N_i / (ω_max_i-ω_min_i), where N_i is the number of poles in the i-th band; calculate the band performance deviation Δ_i=|Performance_current_i-Performance_target_i| / Performance_target_i; adaptive weight calculation w_i=(ρ_i×Δ_i) / Σ j (ρ_j×Δ_j), and normalized Σ iw_i=1; differentiated weight corrections are applied to heterogeneous modules: the low-frequency weight of lithium iron phosphate modules is increased by 20%, and the high-frequency weight of ternary lithium modules is increased by 15%. ω_max_i is the maximum frequency of the i-th frequency band, ω_min_i is the minimum frequency of the i-th frequency band, Performance_current_i is the actual performance indicator of the current frequency band, and Performance_target_i is the target performance indicator. The adaptive band weight matrix and extreme band assignment table are read, and a three-layer frequency domain feature set is obtained through feature vector organization and hierarchical construction. Construction process: Create the feature vector F_low={Ξ_low[], w_low, ζ_target_low, ω_target_low}, F_mid={Ξ_mid[], w_mid, ζ_target_mid, ω_target_mid}, F_high={Ξ_high[], w_high, ζ_target_high, ω_target_high} for each frequency band; set the target parameters of each frequency band: low-frequency band target damping ratio ζ_target_low=0.8, target frequency ω_target_low=2π×0.5rad / s, mid-frequency band target damping ratio ζ_target_mid=0.6, target frequency ω_target_mid=2π×3rad / s, high-frequency band target damping ratio ζ_target_high=0.4, target frequency ω_target_high=2π×8rad / s.

[0060] Read the three-layer frequency domain feature set, decompose and reconstruct it through the standard function decoupling algorithm, and obtain the decoupled objective function group {J_decoupled}. Core decoupling process: Analyze the original unified objective function J_total=Σ i |Ξ i -Ξ ref | 2The coupling terms in [C] are used to identify the cross-band correlation matrix C_cross = [c_ij], where c_ij represents the coupling strength between the i-th and j-th bands. A decoupling transformation matrix T_decouple = I - α × C_cross is established, where I is the identity matrix and α = 0.1 is the decoupling strength coefficient. The decoupling transformation is applied to obtain independent objective functions: the independent objective function for the low-frequency band is J_low_independent = (T_decouple × J_low_original), the independent objective function for the mid-frequency band is J_mid_independent = (T_decouple × J_mid_original), and the independent objective function for the high-frequency band is J_high_independent = (T_decouple × J_high_original). The decoupling effect is verified by calculating the post-decoupling correlation coefficient matrix, requiring that the off-diagonal elements are less than 0.05. J_low_original is the original objective function for the low-frequency band, J_mid_original is the original objective function for the mid-frequency band, and J_high_original is the original objective function for the high-frequency band. Read the decoupled objective function group and the adaptive frequency band weight matrix, and obtain the layered optimization objective function set J_layered through the parameter mapping algorithm and weight optimization allocation. Processing process: Establish a parameter mapping relationship, map the physical parameters to the optimization variable space, low-frequency band optimization variable x_low=[ζ_low, ω_n_low], mid-frequency band optimization variable x_mid=[ζ_mid, ω_n_mid], high-frequency band optimization variable x_high=[ζ_high, ω_n_high]; construct the final layered objective function J_layered={J_low(x_low), J_mid(x_mid), J_high(x_high)}, where J_low(x_low)=w_low×Σ i ∈F_low[(ζ i -ζ_target_low) 2 +(ω_n_i-ω_target_low) 2 ], and other frequency bands are constructed similarly; the dimension normalization coefficient of the objective function is set to ensure that the numerical range of the objective function of each frequency band is consistent.

[0061] The calculation of the decoupling transformation matrix C_cross is as follows: Calculate the cross-correlation of the frequency domain response, the response of each module at frequency ω: H1(jω)=0.85 / (jω+2.5) (LFP module); H3(jω)=1.2 / (jω) 2 +6.4(jω)+26.57)(NCM module); construct the cross-band correlation matrix C_cross; c_ij=|∫[H i(jω)×H j *(jω)]dω| / (||H i ||×||H j ||); where the integration range is the frequency interval of each frequency band; the actual calculation result C_cross=[[1.0, 0.15, 0.23, 0.18], [0.15, 1.0, 0.31, 0.28], [0.23, 0.31, 1.0, 0.42], [0.18, 0.28, 0.42, 1.0]]; determine the decoupling parameter α; minimize the condition number by α=0.6; T_decouple=I-0.6×C_cross.

[0062] like Figure 4 As shown, according to one aspect of the present application, obtaining a frequency band matching parameter matrix includes:

[0063] S41, discretizing the continuous constraint space into grid points, adaptively adjusting the grid density according to the constraint gradient, performing constraint checks on each grid point, and constructing a pre-calculated constraint feasible domain database and a spatial index tree;

[0064] S42, performing optimization calculations by parallel gradient projection based on the hierarchical optimization objective function set and the pre-computed constrained feasible domain database to obtain parallel optimization intermediate results; including: starting three parallel threads to process low-frequency band, mid-frequency band, and high-frequency band optimization respectively, each thread independently performing a gradient projection operation based on the pre-computed constrained feasible domain database;

[0065] S43. Based on the intermediate results of the parallel optimization, the priority negotiation mechanism parameters are adjusted through the thread synchronization mechanism, and the frequency band matching parameter matrix is ​​obtained by summarizing.

[0066] Specifically, the hierarchical constraint condition set is read, and the parallel computing configuration parameters are obtained through multi-threaded task allocation and resource pre-allocation. The optimization tasks of the three frequency bands are respectively assigned to independent computing threads to avoid mutual interference. The parallel computing configuration parameters and the hierarchical optimization objective function set are read, and the improved parallel gradient projection algorithm is used to solve them to obtain the sub-band pole configuration results. The pre-calculation constraint feasible region method is adopted to convert the constraint check from online calculation to offline table lookup, which significantly reduces the calculation time. The gradient projection formula is: Ξ(k+1)=P_Ω[Ξ(k)-α▽J(Ξ(k))], where P_Ω is the projection operator on the constraint set Ω and α is the step size parameter. The sub-band pole configuration results are read, and the band matching parameter matrix is ​​obtained through convergence judgment and cross-band consistency test. Set the convergence judgment condition: ||Ξ(k+1)-Ξ(k)||2≤ε, where ε=1e -4is the convergence threshold. When a frequency band fails to converge, the iteration continues using the progressive constraint relaxation strategy. Where Ξ(k) is the pole distribution matrix at the kth iteration, ▽ is the gradient, and J() is the objective function.

[0067] According to one aspect of the present application, constructing a spatial index tree includes:

[0068] Based on the hierarchical optimization objective function set, the constraint gradients of each region in the constraint space are calculated to identify the regions with drastic constraint changes and the regions with gentle constraint changes.

[0069] Based on the constraint gradient, mesh refinement is implemented in areas where the constraints change drastically, and mesh sparsification is performed in areas where the constraints change gently, forming an adaptive mesh structure.

[0070] Based on the adaptive grid structure, a spatial coordinate index is established for each grid point, and a spatial index tree that supports fast neighbor search is constructed.

[0071] According to one aspect of the present application, summarizing the frequency band matching parameter matrix includes:

[0072] Based on the intermediate results of parallel optimization, the minimum spacing between the poles of adjacent frequency bands is checked to identify boundary conflicts caused by the poles being too close (less than a threshold);

[0073] When a boundary conflict occurs, the priority of frequency band adjustment is determined according to the module importance. Frequency bands with a priority higher than the threshold remain unchanged, while frequency bands with a priority lower than the threshold undergo parameter configuration adjustment.

[0074] Multi-thread synchronous collection is achieved through the exchange of boundary information between threads, and the adjusted optimization results of each frequency band are unified and integrated to obtain the frequency band matching parameter matrix.

[0075] Specifically, the hierarchical constraint set and pre-stored physical parameter boundaries are read and constructed offline using the feasible domain pre-computation algorithm to obtain a pre-computed constrained feasible domain database. The pre-computation process: discretize the continuous constraint space into grid points, and the grid density is adaptively adjusted according to the constraint gradient, with high grid density in areas with large gradients. Constraint checks are performed on each grid point (ζ_grid, ω_grid): stability constraint Re(Ξ) < 0, performance constraint 0.2≤ζ≤1.0, frequency constraint 0.1≤ω_n≤10, response time constraint t_s≤2 / ζω_n≤100ms. A feasible domain lookup table LUT_feasible is constructed, with the storage format {grid_point, feasibility_flag, distance_to_boundary}. A spatial index tree structure is established for fast search, supporting nearest neighbor search with O(logn) complexity. Where grid_point is the grid point, feasibility_flag is the feasibility flag, and distance_to_boundary is the distance to the constraint boundary. The grid adaptation process is as follows: Initial grid: droop parameter range [0.1, 2.0] × [0.1, 2.0], density 20 × 20; Compute constraint gradient gradg(x, y) = [dg / dx, dg / dy]; For example, the gradient at a point gradg = [0.15, 0.23]; Gradient modulus ||gradg|| = 0.275; Adaptive adjustment criteria: if ||gradg|| > 0.2: Grid density increased by × 2; if ||gradg|| < 0.05: Grid density reduced by ÷ 2; Result: The region is resized to a 40 × 40 grid. Feasible region boundary detection: Constraint: K_p + K_i ≤ 1.5; In this embodiment, the feasible point set: {(0.8, 0.6), (1.0, 0.4), (0.5, 0.9), ...}.

[0076] Read the pre-calculated constraint feasible domain database and the hierarchical optimization objective function set, perform optimization calculations through the parallel gradient projection algorithm, and obtain the parallel optimization intermediate result x_parallel. Parallel processing flow: Start three parallel threads to handle low, medium, and high frequency band optimization respectively; each thread independently calculates the objective function gradient ▽J_i(x_i) using the finite difference method: ▽J_i≈[J_i(x_i+Δe_j)-J_i(x_i-Δe_j)] / (2Δ), where Δ=1e -6, e_j is a unit vector; perform the projection operation x_new = P_Ω[x_old - α▽J(x_old)] . Projection is achieved by searching the precomputed feasible region: searching LUT_feasible to find the nearest feasible point. Projection time complexity is O(logn). Threads exchange boundary information through shared memory to ensure solution consistency, with a synchronization period of every 10 iterations. x_new is the new optimization variable, and x_old is the optimization variable for the current iteration.

[0077] Read the intermediate results of parallel optimization and, through thread synchronization and result aggregation, obtain the sub-band pole configuration result Ξ_config. During the synchronization process, a synchronization barrier is set to ensure that all three threads complete the current iteration simultaneously. Cross-band boundary constraints are checked: adjacent band poles cannot be too close together, with a minimum spacing of Δω_min = 0.5 rad / s. When boundary conflicts arise, a negotiation mechanism is adopted: high-priority bands (ordered by module importance) remain unchanged, while low-priority bands adjust their parameters. The results of each thread, Ξ_config = {Ξ_low_final, Ξ_mid_final, Ξ_high_final, convergence_flags, iteration_counts}, are collected, and the convergence status and iteration count are recorded. The output format includes the final position of each pole, the corresponding physical parameter value, and the optimization quality indicator. Among them, Ξ_low_final is the final point configuration result of the low frequency band, Ξ_mid_final is the final point configuration result of the mid frequency band, Ξ_high_final is the final point configuration result of the high frequency band, convergence_flags is the convergence flag, and iteration_counts is the number of iterations.

[0078] In another embodiment of the present application, the versioned synchronization mechanism is as follows: initial state: 4 parallel threads process 4 frequency band optimizations respectively; version number management thread_version = [v1.0, v1.0, v1.0, v1.0]; global_version = 1.0; CAS operation of parameter exchange: thread 1 completes low frequency band optimization and prepares to update boundary parameters: if (compare_and_swap(global_version, 1.0, 1.1)) { / / Successfully obtained write permission update_boundary_params(freq_low_boundary, new_value); broadcast_version_update(1.1);}else { / / Version conflict, rollback and recalculate rollback_to_version(1.0);}; Conflict detection: When thread 2 and thread 3 try to update at the same time: thread 2: expected version 1.0 → actual version 1.1 (conflict); processing: thread 2 rolls back and recalculates based on version 1.1; result: parameter consistency is ensured.

[0079] According to one aspect of the present application, after obtaining the droop parameter adjustment instruction, a physical feasibility test is also included:

[0080] Based on the droop parameter adjustment instructions, multi-dimensional physical constraint checks are performed, including droop coefficient value range checks, response time limit checks, and system stability margin checks;

[0081] Verify the parameter adaptability of each module based on the different chemical types of heterogeneous battery modules, and apply corresponding constraint boundary conditions for lithium iron phosphate modules and ternary lithium modules;

[0082] The parameters that have passed the physical constraint check and adaptability verification are packaged in an output format, including the specific parameter values ​​of each module, execution timestamp and integrity check code.

[0083] Specifically, the frequency band matching parameter matrix is ​​read and cross-band fusion is performed through an adaptive weight allocation algorithm to obtain a unified droop parameter matrix R_unified. The weighted fusion formula is: R_unified = w_low × R_low + w_mid × R_mid + w_high × R_high, and the weight coefficient is dynamically adjusted according to the convergence quality of each frequency band: w_i = (1 / ε_i) / Σ j(1 / ε_j), where ε_i is the convergence error for each frequency band. The unified droop parameter matrix is ​​read and, through the physical limitations of the heterogeneous modules and safety boundary constraints, the achievable droop parameter set R_feasible is obtained. Verification conditions include: droop coefficient range 0.01 ≤ R_droop ≤ 1.0, response time limit t_response ≤ 100ms, and stability margin PM ≥ 45°. The achievable droop parameter set is read and, through module address mapping and instruction format encapsulation, the final droop parameter adjustment instruction R_droop_final is obtained. The output format is: R_droop_final = [R_module1, R_module2, ..., R_moduleN, timestamp, checksum], which contains the specific droop parameter values, timestamp, and checksum for each heterogeneous module. R_low, R_mid, and R_high represent the droop parameter matrices for the low, mid, and high frequency bands, respectively.

[0084] According to one aspect of the present application, the weighted fusion process includes calculating the fusion weight, specifically:

[0085] Based on the frequency band matching parameter matrix, the convergence error, standard deviation of the solution and parameter stability of each frequency band are calculated respectively, and a multi-dimensional convergence quality evaluation index is established;

[0086] Based on the multi-dimensional convergence quality evaluation index, the comprehensive quality evaluation value of each frequency band is calculated and the weight coefficient of each frequency band in the fusion process is determined;

[0087] Based on the weight coefficient, a protective weight reduction is implemented for the frequency bands whose comprehensive quality evaluation values ​​are lower than the threshold, and an adjusted fusion weight coefficient matrix is ​​obtained.

[0088] Specifically, read the frequency band pole configuration results, and obtain the fusion weight coefficient matrix W_fusion through the convergence quality analysis algorithm and dynamic weight calculation. Quality assessment process: Calculate the convergence quality index Q_i=1 / (1+ε_i+σ_i) of each frequency band, where ε_i is the convergence error and σ_i is the standard deviation of the solution; evaluate parameter stability S_i=exp(-|Ξ_final_i-Ξ_target_i| / Ξ_target_i); Comprehensive quality evaluation M_i=η×Q_i+(1-η)×S_i, where η=0.6 is the trade-off coefficient; Dynamic weight calculation w_fusion_i=M_i / Σ jM_j, ensure that the weight sum is 1; for frequency bands with too low quality (M_i<0.3), a protective weight reduction is used to prevent unstable solutions from affecting the overall performance. Read the fusion weight coefficient matrix and the sub-band pole configuration results, and obtain the unified droop parameter matrix R_unified through weighted fusion algorithm and consistency verification. Fusion calculation process: Execute weighted fusion R_unified=Σ i w_fusion_i × R_i, where R_i is the droop parameter for each frequency band. Parameter consistency check: Check the physical rationality of the fused parameters, requiring the droop coefficient to be 0.01 ≤ R_droop ≤ 1.0. Verify system stability: Calculate the eigenvalues ​​of the fused closed-loop system to ensure that all poles are located in the left half plane. Perform performance check: Overshoot σ ≤ 20%, Settling time t_s ≤ 100 ms, and Steady-State Error e_ss ≤ 2%. If verification fails, adjust the fusion weights and recalculate, iterating up to three times. The output format is R_unified = {R_droop_matrix, stability_margin, performance_indices, fusion_quality}. Where Ξ_final_i is the final point configuration for the i-th frequency band, Ξ_target_i is the target pole configuration for the i-th frequency band, R_droop_matrix is ​​the unified droop parameter matrix, stability_margin is the stability margin, performance_indices is the set of performance indicators, and fusion_quality is the fusion quality indicator.

[0089] The dynamic weight fusion process in quality assessment is as follows: The optimization result error of each frequency band is: ε1=||x_k+1-x_k|| / ||x_k||=0.023 (ultra-low frequency band); ε2=0.031 (low frequency band); ε3=0.015 (mid-frequency band); ε4=0.042 (high frequency band); S_i=1-|Δparameter_i| / parameter_i; S1=1-|0.05| / 1.2=0.958; S2=1-|0.08| / 1.5=0.947; S3=1-|0.03| / 0.8=0.963; S4=1-|0.12| / 2.1=0 .943; Mi = (1-εi) × Si × ηi; trade-off coefficient η = [0.3, 0.3, 0.25, 0.15]; M1 = (1-0.023) × 0.958 × 0.3 = 0.281; M2 = (1-0.031) × 0.947 × 0.3 = 0.275; M3 = (1-0.015) × 0.963 × 0.25 = 0.237; M4 = (1-0.042) × 0.943 × 0.15 = 0.136; dynamic weight update wi = Mi / Σ(Mi); final weight w = [0.303, 0.296, 0.255, 0.146]. Where parameter_i is the reference parameter value for the current frequency band, and Δparameter_i is the change in the parameter value.

[0090] According to another aspect of the present application, the process of obtaining the pole distribution matrix may also be:

[0091] Based on the dynamic transfer function parameter set, the system order of each heterogeneous battery module is identified through the transfer function order analysis algorithm, and the module order distribution vector and heterogeneous type identification are obtained;

[0092] Based on the module order distribution vector and dynamic transfer function parameter set, the pole number of modules of different orders is aligned through the virtual pole generation algorithm, in which the low-order modules supplement the virtual poles to the highest order to obtain the aligned pole set;

[0093] Based on the aligned pole set and heterogeneous type identification, a pole importance evaluation algorithm is used to assign a first weight of 1.0 to actual poles and a second weight of 0.01 to virtual poles, constructing a standardized pole distribution matrix and a module pole mapping table. The second weight is smaller than the first weight.

[0094] Specifically, read the dynamic transfer function parameter set. Through the heterogeneous module pole unification processing algorithm, obtain the standardized pole distribution matrix and the module pole mapping table: Read the dynamic transfer function parameter set. Through the transfer function order analysis algorithm, identify the system order of each module, and obtain the module order distribution vector and the heterogeneous type identifier. The processing process is as follows: Analyze the denominator polynomial order n_i of each module in the dynamic transfer function parameter set G_params; Identify the module type: first-order system (n = 1), second-order system (n = 2), high-order system (n > 2); Output the module order distribution vector order_vector = [n_1, n_2,..., n_N], and the heterogeneous type identifier hetero_flags = [type_1, type_2,..., type_N]. Read the module order distribution vector and the dynamic transfer function parameter set, and perform pole number alignment through the virtual pole generation algorithm to obtain the aligned pole set Ξ_aligned. The processing process is as follows: Determine the maximum order: n_max = max(order_vector); For the module with order n_i < n_max, generate (n_max - n_i) virtual poles; The virtual pole position is set to: Ξ_virtual = -100 ± j×0 (far from the working frequency band); Construct the alignment matrix: Ξ_aligned[i, :] = [Ξ_real_i, Ξ_virtual_i]. Read the aligned pole set and the heterogeneous type identifier, and generate a weighted pole matrix through the pole importance evaluation algorithm to obtain the standardized pole distribution matrix Ξ_matrix_unified and the module pole mapping table Ξ_mapping. The processing process is as follows: Assign a weight w_real = 1.0 to the actual pole, and a weight w_virtual = 0.01 to the virtual pole; Construct the mapping table Ξ_mapping to record the source module and the authenticity identifier of each pole; Output Ξ_matrix_unified = [Ξ_aligned, weight_matrix], with a dimension of N×(2×n_max); Output Ξ_mapping = {module_id, pole_index, is_real, weight}. Where Ξ_real_i is the actual pole of the i-th module, Ξ_virtual_i is the virtual pole of the i-th module, module_id is the module identifier, pole_index is the pole index, weight_matrix is the weight matrix, is_real is the authenticity identifier, and weight is the pole weight.

[0095] According to another aspect of the present application, determine the frequency band attribution characteristics of each pole, or calculate the frequency domain classification parameters of each pole, including:

[0096] Based on the standardized pole distribution matrix and the module pole mapping table, the effective frequency of each pole is calculated using a weighted frequency extraction algorithm, where virtual poles are excluded through weight filtering to obtain a weighted pole frequency matrix;

[0097] Based on the weighted pole frequency matrix, the frequency band boundaries are dynamically determined by the K-means clustering algorithm, requiring that the interval between adjacent frequency band boundaries is greater than 0.5 decades, and the dynamic frequency band boundary vector is obtained;

[0098] Based on the weighted pole frequency matrix, dynamic frequency band boundary vector and module pole mapping table, the poles are assigned to the corresponding frequency bands through the heterogeneous perception classification algorithm, maintaining the frequency band consistency of the conjugate pole pairs, and obtaining the pole frequency band attribution table containing the correspondence between the poles and the modules.

[0099] In other words, the frequency domain characteristic parameters calculated based on the position of the pole in the complex plane are used to determine the frequency band to which the pole belongs in the multi-band control strategy.

[0100] Specifically, the standardized pole distribution matrix and module pole mapping table are read and layered using a weighted frequency domain analysis algorithm to obtain a three-layer frequency domain feature set. The standardized pole distribution matrix and module pole mapping table are read and the effective pole frequencies are calculated using a weighted frequency extraction algorithm to obtain a weighted pole frequency matrix ω_weighted. The processing is as follows: extract the frequency of each pole: ω_i = |Im(Ξ_i)|; apply importance weights: ω_effective_i = ω_i × weight_i; filter virtual poles: if weight_i < 0.1, then ω_effective_i = NaN; output ω_weighted = [ω_effective, module_id, original_index]. The weighted pole frequency matrix is ​​read and the frequency band boundaries are dynamically determined using a cluster analysis algorithm to obtain the dynamic frequency band boundary vector freq_boundaries. The processing procedure is as follows: K-means clustering (K=4) is performed on the effective pole frequencies; frequency band boundaries are determined based on the cluster centers to avoid pole segmentation caused by fixed boundaries; the rationality of the boundaries is verified: the interval between adjacent frequency band boundaries is greater than 0.5 decade; and the output is freq_boundaries = [ω_ultra_low, ω_low, ω_mid, ω_high]. The weighted pole frequency matrix, dynamic frequency band boundary vector, and module pole mapping table are read, and frequency domain features are generated using the heterogeneous perception classification algorithm to obtain a three-layer frequency domain feature set F_layers. The processing procedure is as follows: poles are assigned to corresponding frequency bands based on freq_boundaries; conjugate pole pairs of the same module are kept in the same frequency band; the pole distribution of different module types in each frequency band is recorded; and the output is F_layers = {F_ultralow[poles, modules], F_low[poles, modules], F_mid[poles, modules], F_high[poles, modules]}. Where weight_i is the weight value of each pole, ω_effective_i is the weighted effective pole frequency, and NaN is used to mark invalid data. For example, when weight_i < 0.1, the pole is filtered out and assigned a value of NaN. original_index is the index position of the pole in the original data. ω_ultra_low is the boundary frequency of the ultra-low frequency band. F_ultralow[poles, modules] is the pole of the ultra-low frequency band and the corresponding module information.

[0101] In another embodiment of the present application, the boundaries are determined from the cluster centers, as follows: based on the effective pole frequency data of the aforementioned embodiments: ω_effective1=0rad / s (real part dominant pole) of module 1; ω_effective2=0rad / s (real part dominant pole) of module 2; ω_effective3=4.1rad / s (complex conjugate pole pair) of module 3; ω_effective4=4.0rad / s (complex conjugate pole pair) of module 4; the input data set after filtering the virtual poles: ω_input=[0, 0, 4.1, 4.0]rad / s. Execute K-means clustering (K=4) calculation process: Initial cluster centers are randomly distributed: [0.01, 1.0, 4.0, 10.0] rad / s; First iteration: assign data point [0, 0] to cluster 1, [4.1, 4.0] to cluster 3; Update cluster centers: Cluster 1 center = (0+0) / 2=0, Cluster 3 center = (4.1+4.0) / 2=4.05; Final cluster centers after convergence: [0, 0.8, 4.05, 15.3 ]rad / s; The mapping relationship between cluster centers and boundaries is as follows: The geometric mean is used to determine the boundaries of adjacent frequency bands: The upper bound of the ultra-low frequency band is sqrt(0×0.8)=0 (special treatment is 0.1rad / s); the upper bound of the low frequency band is sqrt(0.8×4.05)=1.8rad / s; the upper bound of the mid-frequency band is sqrt(4.05×15.3)=7.9rad / s; the lower bound of the high frequency band is 15.3rad / s (open interval); the boundary interval verification is log 10 (1.8 / 0.1)=1.26>0.5decade;log 10 (7.9 / 1.8)=0.64>0.5decade;log 10 (15.3 / 7.9)=0.29<0.5decade, adjusted to log 10(20.0 / 7.9)=0.40; the final output freq_boundaries=[0.1,1.8,7.9,20.0]rad / s. Determine which frequency band the pole belongs to, as follows: Based on the determined frequency band boundaries [0.1,1.8,7.9,20.0]rad / s, perform pole allocation: Frequency band allocation rule: when ω_effective≤0.1, it belongs to F_ultralow; when 0.1<ω_effective≤1.8, it belongs to F_low; when 1.8<ω_effective≤7.9, it belongs to F_mid; when ω_effective>7.9, it belongs to F_high; Actual allocation process: Module 1 pole: ω_effective1=0≤0.1→ belongs to F_ultralow; Module 2 pole: ω_effective2=0≤0.1→ belongs to F _ultralow; Module 3 pole: ω_effective3=4.1, satisfying 1.8<4.1≤7.9→belongs to F_mid; Module 4 pole: ω_effective4=4.0, satisfying 1.8<4.0≤7.9→belongs to F_mid; If the pole frequency is exactly on the boundary, the processing process is as follows: Set the boundary tolerance δ=0.01rad / s, and when |ω_effective-boundary|≤δ, it is considered to be on the boundary: Boundary processing rules: If the poles come from the conjugate pole pair of the same module, they are kept in the same frequency band first; if the poles come from different modules, the belonging is determined according to the module importance ranking; if the importance is the same, the proximity principle is used to assign them to the frequency band with a closer numerical distance. For example: Assume that the pole frequency of module 5 is ω_effective5 = 1.81 rad / s ≈ 1.8 (boundary value); calculate the distance |1.81-1.8| = 0.01≤δ, which is determined to be a boundary condition; compare the frequency band distances, the distance to the F_low center vs. the distance to the F_mid center; F_low center estimate: (0.1+1.8) / 2=0.95, distance = |1.81-0.95|=0.86; F_mid center estimate: (1.8+7.9) / 2=4.85, distance = |1.81-4.85|=3.04; the conclusion is that it belongs to F_low (the distance is closer). If there is no actual frequency data, the process is as follows: When some modules have only virtual poles (weight < 0.1) resulting in no valid frequency data: Default frequency band boundary settings: Theoretical boundaries based on the physical characteristics of the battery system: F_ultralow: [0.001, 0.1] rad / s (SOC equilibrium time scale); F_low: [0.1, 1.0] rad / s (thermal dynamic time scale); F_mid: [1.0, 10.0] rad / s (electrochemical dynamic time scale); F_high: [10.0, 100.0] rad / s (electrical dynamic time scale).For example, when only 1-2 valid frequencies remain after filtering, the theoretical boundary is used; when the number of valid frequencies is ≥3, K-means is used to dynamically determine the boundary; in mixed mode, the theoretical boundary is used as a constraint and the K-means result is used as the optimization target. Based on the above allocation results, the frequency domain feature set is constructed: F_layers={F_ultralow[poles=[-2.5,-2.3], modules=[LFP_1, LFP_2]], F_low[poles=[], modules=[]], F_mid[poles=[-3.2±j4.1,-3.0±j4.0], modules=[NCM_3, NCM_4]], F_high[poles=[], modules=[]]}. Module pole mapping table update: Ξ_mapping={{module_id=1, pole_index=1, frequency_band="ultralow", weight=1.0}, {module_id=2, pole_index=1, frequency_band="ultralow", weight=1.0}, {module_id=3, pole_index=1, frequency_band="mi d", weight=1.0}, {module_id=3, pole_index=2, frequency_band="mid", weight=1.0}, {module_id=4, pole_index=1, frequency_band="mid", weight=1.0}, {module_id=4, pole_index=2, frequency_band="mid", weight=1.0}}.

[0102] According to another aspect of the present application, the frequency band matching parameter matrix may also be obtained as follows:

[0103] Based on the hierarchical optimization objective function set, a shared memory region and a mutex matrix are created through a thread-safe initialization algorithm, and a shared memory manager is established to store frequency band boundary information.

[0104] Based on the shared memory manager, mutex matrix and the optimized state of each thread, the boundary parameters are exchanged through a versioned synchronization algorithm, in which a read-write lock mechanism and version numbers are used to prevent data competition and obtain the synchronization boundary parameters.

[0105] Based on the synchronization boundary parameters and thread priority table, the boundary constraint conflict is resolved through the conflict negotiation algorithm. The high-priority thread parameters remain unchanged, and the low-priority thread execution is rolled back. The coordinated boundary parameters and thread synchronization status table are obtained.

[0106] Based on the coordinated boundary parameters, thread synchronization state table and local optimization results of each thread, the global solution is updated using the Compare-And-Swap operation through a lock-free aggregation algorithm to obtain the frequency band matching parameter matrix.

[0107] Specifically, the pre-calculated constraint feasible domain database and the hierarchical optimization objective function set are read, and the optimization calculation is performed through the parallel gradient projection algorithm with synchronization mechanism to obtain the parallel optimization intermediate results and thread synchronization state table. The hierarchical optimization objective function set is read, and the shared memory structure is established through the thread-safe initialization algorithm to obtain the shared memory manager shared_mem_manager and the mutex matrix mutex_matrix. Processing process: Create a shared memory area: boundary_info[3][3] to store frequency band boundary information; initialize the read-write lock: create an independent mutex_matrix[i][j] for each boundary; set the atomic operation flag: atomic_flags[3] to indicate the update status of each thread; establish a memory barrier: ensure cross-thread memory visibility; output shared_mem_manager={boundary_info, version_numbers, timestamp}. Read the shared memory manager shared_mem_manager, the mutex matrix mutex_matrix and the current thread optimization state thread_state, and coordinate the boundaries through the versioned synchronization algorithm to obtain the synchronization boundary parameter boundary_synchronized. Processing process:

[0108] foreachthread_iinparallel:

[0109] / / Read adjacent band boundaries (using read lock);

[0110] acquire_read_lock(mutex_matrix[i][i-1]);

[0111] local_boundary_low=boundary_info[i-1][i];

[0112] version_low=version_numbers[i-1][i];

[0113] release_read_lock(mutex_matrix[i][i-1]);

[0114] / / Calculate local updates;

[0115] new_boundary = compute_local_update(thread_state[i]);

[0116] / / Write the update (using write lock + version check);

[0117] acquire_write_lock(mutex_matrix[i][i]);

[0118] if (version_numbers[i][i] == local_version):

[0119] boundary_info[i][i] = new_boundary;

[0120] version_numbers[i][i]++;

[0121] atomic_flags[i] = UPDATED;

[0122] else:

[0123] / / Version conflict, need to recompute;

[0124] conflict_flag = true;

[0125] release_write_lock(mutex_matrix[i][i]).

[0126] Read the synchronized boundary parameter boundary_synchronized, atomic operation flag atomic_flags, and thread priority table thread_priority. Solve the race condition through a conflict negotiation algorithm to obtain the coordinated boundary parameter boundary_coordinated and thread synchronization status table sync_status. Processing process: Detect conflicts: Triggered when |boundary[i] - boundary[i - 1]| < threshold; Priority arbitration: Decide which value to retain based on thread_priority and convergence quality; Rollback mechanism: Low-priority threads roll back to the previous stable state; Global synchronization point: Perform a global barrier synchronization every 10 iterations; Output sync_status = {conflict_count, rollback_count, sync_timestamp}.

[0127] The coordinated boundary parameter boundary_coordinated, the thread synchronization status table sync_status, and each thread's local optimization results local_results are read. A lock-free aggregation algorithm is used to generate the final result, resulting in the parallel optimization intermediate result x_parallel. Processing: Use Compare-And-Swap (CAS) operations to update the global optimal solution. Each thread writes the results to a separate buffer to avoid write conflicts. The main thread polls atomic_flags and waits for all threads to complete. Aggregate the results: x_parallel = aggregate(local_results, boundary_coordinated). Record the parallel efficiency metrics: speedup_ratio and conflict_ratio.

[0128] In another embodiment of the present application, the normalization process can also include: reading the synchronized raw data set and SOC value, processing them through the SOC adaptive normalization algorithm, and obtaining the normalized module state matrix X_state. Processing process: SOC segmented normalization coefficient design: low SOC range (0-20%): γ_low = 1.2 (enhanced sensitivity); medium SOC range (20-80%): γ_mid = 1.0 (standard response); high SOC range (80-100%): γ_high = 0.8 (reduced sensitivity). Differentiated SOC processing of heterogeneous modules: lithium iron phosphate module: α_LFP(SOC)=α_LFP_base×(1+0.2×(0.5-SOC)); ternary lithium module: α_NCM(SOC)=α_NCM_base×(1+0.3×(SOC-0.5)); standardization formula: X_state[i]=(X_raw[i]-μ[i]) / (σ[i]×α_type(SOC[i])×γ_SOC_range). Among them, α_LFP_base is the basic adjustment coefficient of the lithium iron phosphate module; α_LFP(SOC) is the adjustment coefficient of the lithium iron phosphate module; SOC is the state of charge of the battery; α_NCM(SOC) is the adjustment coefficient of the ternary lithium module, and α_NCM_base is the basic adjustment coefficient of the ternary lithium module; X_raw[i] is the value of the i-th sample in the original data set, μ[i] is the mean of the i-th sample, σ[i] is the standard deviation of the i-th sample, α_type is the corresponding adjustment coefficient selected according to the module type; γ_SOC_range is the SOC segment normalization coefficient.

[0129] In another embodiment of the present application, the physical feasibility check can also include: adding parameter constraints related to the SOC operating range: SOC sensitive range protection: when any module's SOC is less than 0.2 or greater than 0.8; lowering the droop coefficient upper limit to R_droop_max = 0.5 (normally 1.0); and relaxing the response time requirement: t_response_max = 150ms (normally 100ms). SOC balancing rate constraints: calculating the maximum allowable balancing current based on the current SOC distribution; and verifying that droop parameters do not result in excessive balancing current.

[0130] In another embodiment of the present application, the process of obtaining the extreme frequency band attribution table can also be: based on the theoretical analysis of the multi-time scale dynamic characteristics of the battery system, a three-layer frequency domain division strategy is established: extremely low frequency band (0.001-0.1Hz, 0.00628-0.628rad / s): corresponding to the SOC balance process, time constant τ_SOC=100-1000s; low frequency band (0.1-1Hz, 0.628-6.28rad / s): corresponding to the thermal dynamic process, time constant τ_thermal=1-10s; medium frequency band (1-10Hz, 6.28-62.8rad / s): corresponding to the electrochemical dynamic process, time constant τ_electrochemical=0.1-1s; high frequency band (10-100Hz, 62.8-628rad / s): corresponding to the electrical dynamic process, time constant τ_electrical=0.01-0.1s. Considering the importance of SOC balance, the frequency bands are redefined as: F_ultralow (0.001-0.1Hz): SOC balance control; F_low (0.1-1Hz): thermal management coordination; F_mid (1-10Hz): power distribution adjustment; F_high (10-100Hz): transient response suppression.

[0131] Pole frequency domain assignment and classification: Read the pole distribution matrix Ξ_matrix and process it using a physical time-scale-based pole classification algorithm to obtain a four-layer pole frequency band assignment table Ξ_classification. Processing: Calculate the time constant corresponding to each pole: τ_i = -1 / Re(Ξ_i). Assign frequency bands based on the time constants: When τ_i > 100s, assign to the ultra-low frequency band F_ultralow (SOC equilibrium); when 10s < τ_i ≤ 100s, assign to the low frequency band F_low (thermal dynamics); when 0.1s < τ_i ≤ 10s, assign to the mid-frequency band F_mid (electrochemical dynamics); when τ_i ≤ 0.1s, assign to the high frequency band F_high (electrical dynamics). Output: Ξ_classification = {Ξ_ultralow[], Ξ_low[], Ξ_mid[], Ξ_high[], time_constant_map}.

[0132] In another embodiment of the present application, the process of obtaining a dynamic forgetting factor sequence is as follows: reading the heterogeneous feature weighted error vector e_weighted, calculating the forgetting factor using an improved error-sensitive adaptive algorithm, and obtaining a dynamic forgetting factor sequence Ξ_adaptive. The specific process is as follows: calculating the normalized error: e_norm(k)=|e_weighted(k)| / (ε+max(|e_weighted|)), where ε=0.001; constructing the inverse forgetting factor function: Ξ_adaptive(k)=Ξ_max-(Ξ_max-Ξ_min)×(1-exp(-β×e_norm(k))); where Ξ_min=0.90 (fast adaptation value for large errors); Ξ_max=0.995 (stable tracking value for small errors); and β=3.0 (error sensitivity coefficient). Sudden error detection and processing: Calculate the error change rate: Δe(k) = |e_weighted(k) - e_weighted(k-1)|; when Δe(k) > θ_jump (θ_jump = 0.5), force Ξ_adaptive(k) = Ξ_min; this enables rapid model reset to adapt to sudden operating conditions. Differentiated processing for heterogeneous modules: Lithium iron phosphate modules: Ξ_range_LFP = [0.92, 0.995] (more stable); ternary lithium modules: Ξ_range_NCM = [0.88, 0.99] (more sensitive). This design ensures: Ξ → Ξ_min (fast adaptation) for large errors, and Ξ → Ξ_max (stable tracking) for small errors.

[0133] In another embodiment of the present application, after constructing the decoupling transformation matrix, it also includes: obtaining the decoupled objective function group {J_decoupled} through the objective function decomposition algorithm based on the modal decoupling theory, specifically: coupling system state space description: considering the state space model of multi-band coupling: x*=Ax+Bu; y=Cx+Du; wherein the state matrix A of the system contains cross-band coupling terms, x* is the time derivative of the system state variable, x is the state vector of the system, y is the output vector, B is the input matrix, u is the input vector, C is the output matrix, and D is the direct transfer matrix. Modal analysis and decoupling: Calculate the eigenvalue decomposition of the system matrix A: A=VΞV -1 ; Modal coordinate transformation: z=V -1 x; in modal coordinates, the system decoupling is: z*=Ξz+V -1 Bu. Where Ξ is a diagonal matrix representing the eigenvalue matrix of the system matrix A; V is the eigenvector matrix of the system matrix A, representing the basis of the modal coordinate transformation; z* is the time derivative of the modal coordinate z. Frequency domain coupling strength quantification: The elements of the cross-band coupling strength matrix C_cross are defined as: c_ij=|∫H_i(jω)×H_j(jω)dω| / sqrt(∫|H_i(jω)| 2 dω×∫|H_j(jω)| 2 dω); where H_i(jω) is the frequency response function of the i-th frequency band. Nonlinear decoupling transformation design: Considering the nonlinear characteristics of the system, the extended decoupling transformation is adopted: T_decouple=I-α1×C_cross-α2×C_cross 2 +α3×diag(C_cross); where α1=0.1 (first-order coupling compensation coefficient); α2=0.01 (second-order coupling compensation coefficient); α3=0.05 (self-coupling enhancement coefficient). Decoupling effect verification: Calculate the residual coupling matrix after decoupling: C_residual=T_decouple T ×C_cross×T_decouple; verification condition: ||C_residual|| F <0.05×||C_cross|| F ; If not satisfied, optimize the α coefficient through iteration. Application of objective function decoupling: decompose the original coupling objective function J_total through modal coordinate transformation: J_decoupled_i=∑_kw_ik×||T_decouple[i,:]×(Ξ_k-Ξ_ref)|| 2 , ensuring the independence of optimization of each frequency band while retaining the necessary coordination constraints. F is the Frobenius norm, Ξ_k is the eigenvalue matrix of the kth mode, and Ξ_ref is the reference eigenvalue matrix.

[0134] In another embodiment of the present application, the numerical implementation of the decoupling transformation is specifically as follows: Coupling strength matrix calculation algorithm: fori=1:num_bands; forj=1:num_bands; / / Calculate the cross-correlation function of frequency bands i and j; correlation=xcorr(H_i, H_j); / / Normalize to obtain the coupling strength; c_ij=max(abs(correlation)) / sqrt(energy_ienergy_j); end; end. Iterative optimization of the decoupling coefficient: while(||C_residual||_F>tolerance); / / Gradient descent update; grad_α1=Ψ||C_residual||_F / Ψα1; α1=α1-η×grad_α1; / / Update the decoupling matrix; T_decouple=update_decouple_matrix(α1, α2, α3); iteration++; end. Stability guarantee: / / Ensure that the decoupling transformation does not cause system instability; eigenvalues ​​= eig(T_decouple); ifany(real(eigenvalues) < 0.5); / / Modify the decoupling strength to ensure stability; α1 = α10.8; end. Where Ψ is the partial derivative.

[0135] In another embodiment of the present application, the traditional forgetting factor λ_traditional = 0.95 (fixed value); the reverse forgetting factor calculation process is as follows: Calculate the prediction error variance: σ 2 _error=Var(y_measured-y_predicted); In the embodiment: σ 2 _error=0.032; calculate the reverse adjustment factor β_reverse=1-exp(-σ 2 _error / 0.1)=1-exp(-0.32)=0.274; calculate the adaptive forgetting factor λ_adaptive=λ_traditional×(1+β_reverse×γ_type); for the LFP module, λ_LFP=0.95×(1+0.274×0.8)=1.158; for the NCM module, λ_NCM=0.95×(1+0.274×1.2)=1.263; generate the dynamic forgetting factor sequence λ_sequence=[1.158, 1.263, 1.158, 1.263] (corresponding to 4 modules).

[0136] This embodiment solves the fundamental problem of inconsistent transfer function orders in heterogeneous battery modules. It ensures mathematical dimensionality uniformity while avoiding the interference of virtual poles on the actual system dynamics. It enables coordinated optimization of the first-order system of lithium iron phosphate batteries and the second-order system of ternary lithium batteries within a unified mathematical framework, eliminating the frequency band resonance problem caused by dimensionality mismatch in traditional methods, improving system stability margins, and effectively preventing high-frequency oscillations when heterogeneous modules are operated in parallel. By constructing a nonlinear decoupling transformation, the previously mutually influencing SOC balancing, power regulation, and transient suppression objectives are completely separated, allowing each frequency band to be optimized independently. This avoids the contradiction in traditional unified optimization where improved SOC balancing in the low-frequency band leads to worsened transient response in the high-frequency band. This shortens the SOC balancing time while reducing transient overshoot, achieving simultaneous improvements in multi-objective performance and fundamentally resolving the problem of inter-coupling of multi-timescale control objectives. It also reduces the computational complexity of online optimization using pre-calculated constrained feasible regions, achieving the engineering feasibility of real-time optimization. When a version conflict is detected, the low-priority thread automatically rolls back and recalculates based on the latest data, avoiding thread blocking caused by traditional mutex mechanisms. This fully utilizes the computing resources of multi-core processors, enabling real-time control of large-scale heterogeneous systems. When convergence quality in a frequency band degrades due to a sudden change in operating conditions, the system automatically reduces the weight of that band to the protection lower limit to prevent local failures from impacting global performance. This reduces parameter configuration deviations in the face of disturbances such as module aging and temperature fluctuations, improving long-term operational stability. The reverse forgetting factor mechanism enables rapid adaptation and stable tracking, enhancing the system's adaptability to changing operating conditions.

[0137] In a specific embodiment of the present application, an energy storage system is provided that includes four heterogeneous battery modules, wherein two lithium iron phosphate (LFP) modules and two ternary lithium (NCM) modules are operated in parallel.

[0138] Step 1: Obtain dynamic response data of heterogeneous modules.

[0139] 1.1. The system is configured with a 16-bit ADC and a 1kHz sampling rate. Data is collected using a hardware synchronization trigger mechanism: Module 1 (LFP): V_raw1 = 3.25V, I_raw1 = 48.5A, T_raw1 = 25.3°C; Module 2 (LFP): V_raw2 = 3.23V, I_raw2 = 47.8A, T_raw2 = 25.1°C; Module 3 (NCM): V_raw3 = 3.78V, I_raw3 = 52.3A, T_raw3 = 26.2°C; Module 4 (NCM): V_raw4 = 3.76V, I_raw4 = 51.9A, T_raw4 = 26.0°C. At the same time, the SOC value is read: SOC_vector = [65%, 63%, 68%, 67%]; the synchronized raw data set Data_sync = {V_raw[1:4], I_raw[1:4], T_raw[1:4], SOC_vector, timestamp = 1698234567}. V_raw is the raw voltage data, I_raw is the raw current data, and T_raw is the raw temperature data.

[0140] 1.2. Identify the module type: Modules 1 and 2 are LFP type, and modules 3 and 4 are NCM type. Calculate the SOC segment normalization coefficient. The SOC of module 3 = 68% belongs to the medium SOC range (20-80%), γ_mid = 1.0: α_NCM(SOC3) = α_NCM_base×(1+0.3×(SOC3-0.5)) = 1.2×(1+0.3×(0.68-0.5)) = 1.2×1.054 = 1.265; normalization processing: calculate the mean: μ3_V = 3.7V (based on historical data); calculate the standard deviation: σ3_V = 0.1V; normalization: X_state[3, V] = (3.78-3.7) / (0.1×1.265×1.0) = 0.632; fully output the normalized module state matrix X_state, dimension 4×3 (4 modules, 3 state quantities).

[0141] Step 2: Identify the transfer function characteristics of each module in real time.

[0142] 2.1. Calculate the signal change rate, using the module 3 voltage as an example: dV / dt = (3.78 - 3.75) / 0.001 = 30 V / s. Dynamically adjust the window length: L_window = L_base × (1 + β × |dV / dt|) = 100 ms × (1 + 0.01 × 30) = 130 ms. Divide X_state into 130 ms windows to obtain the time domain response data segments Y_segments.

[0143] 2.2. Identify the chemical type and set the characteristic coefficients: Modules 1 and 2 (LFP): γ_LFP = 0.8; Modules 3 and 4 (NCM): γ_NCM = 1.2; Calculate the weighted prediction error of module 3: e_weighted(k) = γ_NCM × |y_actual(k) - y_predict(k)| = 1.2 × |3.78 - 3.75| = 0.036. Calculate the normalized error: e_norm(k)=|e_weighted(k)| / (ε+max(|e_weighted|))=0.036 / (0.001+0.045)=0.783; apply the reverse forgetting factor formula: Ξ_adaptive(k)=Ξ_max-(Ξ_max-Ξ_min)×(1-exp(-β×e_norm(k))); where Ξ_max=0.99 is the stable tracking value for small errors; Ξ_min=0.88 is the fast adaptation value for large errors (NCM module); β=3.0 is the error sensitivity coefficient; calculation result: Ξ_adaptive(k)=0.99-(0.99-0.88)×(1-exp(-3.0×0.783))=0.891. Update the gain matrix (taking the second-order system as an example): K(k)=P(k-1)φ(k) / [Ξ_adaptive(k)+φ T (k)P(k-1)φ(k)]; where the regression vector φ(k)=[-y(k-1),-y(k-2),u(k-1),u(k-2)] T ; Update parameters: θ_updated(k)=θ_prev(k-1)+K(k)[y(k)-φ T (k)θ_prev(k-1)]; Convergence test: ||θ_updated(k)-θ_prev(k-1)||2=0.00008<ε_conv(1e -6 ). Convert the parameters to the standard transfer function format: Module 1 (LFP): G1(s)=0.85 / (s+2.5); Module 3 (NCM): G3(s)=1.2 / (s 2 +6.4s+26.57); Output G_params={[b_coeff], [a_coeff], module_id, timestamp}.

[0144] 2.3. Analyze the denominator order of the transfer function: Modules 1 and 2: n = 1 (first-order system); Modules 3 and 4: n = 2 (second-order system); Output order_vector = [1, 1, 2, 2], hetero_flags = [LFP_1st, LFP_1st, NCM_2nd, NCM_2nd]. Determine the maximum order: n_max = max(order_vector) = 2; Add virtual poles to the first-order system: Module 1 Real poles: Ξ 11 =-2.5 (weight w=1.0); module 1 virtual pole: Ξ 12 = -100 + j × 0 (weight w = 0.01); Construct the alignment matrix, module 1: Ξ_aligned[1, :] = [-2.5, -100], weight[1, :] = [1.0, 0.01]. Construct the normalized pole distribution matrix: Ξ_matrix_unified = [[-2.5, -100], [1.0, 0.01], module 1; [-2.3, -100], [1.0, 0.01], module 2; [-3.2 + j4.1, -3.2 - j4.1], [1.0, 1.0], module 3; [-3.0 + j4.0, -3.0 - j4.0], [1.0, 1.0]], module 4. Output module pole mapping table Ξ_mapping records the source and authenticity of each pole.

[0145] Step 3: Hierarchical decoupling and optimization of pole placement targets.

[0146] 3.1. Extracting Pole Frequencies and Applying Weights: Module 1: ω 11 =|Im(-2.5)|=0,ω_effective 11 =0×1.0=0; Module 3: ω 31 =|Im(-3.2+j4.1)|=4.1,ω_effective 31 =4.1×1.0=4.1; filter virtual poles with weights <0.1 and output ω_weighted. Perform K-means clustering (K=4) on the effective pole frequencies: cluster centers: [0.05, 0.8, 4.2, 15.3] rad / s; determine boundaries: freq_boundaries=[0.01, 0.3, 2.0, 10.0] rad / s. Assign poles according to boundaries: ultra-low frequency band F_ultralow: real dominant poles of modules 1 and 2; mid-frequency band F_mid: complex conjugate pole pairs of modules 3 and 4; output F_layers={F_ultralow[Ξ 11 ,Ξ 21 ],F_low[],F_mid[Ξ 31 ,Ξ 32,Ξ 41 ,Ξ 42 ], F_high[]}.

[0147] 3.2. Calculate the cross-band coupling strength, taking the ultra-low frequency and medium frequency as an example: c 13 =|∫H_ultralow(jω)×H_mid*(jω)dω| / sqrt(∫|H_ultralow(jω)| 2 dω×∫|H_mid(jω)| 2 dω); where H_ultralow(jω) is the frequency response of the ultra-low frequency band; H_mid(jω) is the frequency response of the mid-frequency band; the integral is calculated to get: c 13 =0.28. Construct the complete coupling matrix: C_cross=[[0, 0.15, 0.28, 0.08], [0.15, 0, 0.35, 0.12], [0.28, 0.35, 0, 0.25], [0.08, 0.12, 0.25, 0]]. Calculate the decoupling transformation matrix: T_decouple=I-α1×C_cross-α2×C_cross 2 +α3×diag(C_cross); where I is a 4×4 unit matrix; α1=0.1 is the first-order coupling compensation coefficient; α2=0.01 is the second-order coupling compensation coefficient; α3=0.05 is the self-coupling enhancement coefficient. Apply decoupling transformation: J_mid_independent=T_decouple[3,:]×J_mid_original=0.82×J_mid_original. Establish optimization variables: ultra-low frequency band: x_ultralow=[ζ_ultralow, ω_n_ultralow]=[0.9, 0.1]; mid-frequency band: x_mid=[ζ_mid, ω_n_mid]=[0.6, 4.0]; Construct hierarchical objective function: J_mid(x_mid)=w_mid×Σ i ∈F_mid[(ζ i -ζ_target_mid) 2 +(ω_n_i-ω_target_mid) 2 ]=0.35×[(0.52-0.6) 2 +(4.1-4.0) 2 ]=0.0028.

[0148] 3.3. Set frequency band-specific constraints: Very low frequency band: steady-state error Δ_ss ≤ 2%; Medium frequency band: damping ratio ζ ≥ 0.6; High frequency band: response time t_r ≤ 50 ms; Output layered constraint condition set C_layered.

[0149] Step 4: Solve the matching parameters of each frequency band in parallel.

[0150] 4.1. Allocate computing resources: Thread 0: Extremely low frequency band optimization; Thread 1: Low frequency band optimization; Thread 2: Medium frequency band optimization; Thread 3: High frequency band optimization; Set thread priority: thread_priority=[3, 2, 1, 0] (Extremely low frequency band has the highest priority).

[0151] 4.2. Discretize the continuous space, using the mid-frequency band as an example: ζ ranges from [0.2, 1.0] with a step size of 0.01; ω_n ranges from [2.0, 8.0] with a step size of 0.1. Generate 80 × 60 = 4800 grid points. For each grid point (ζ = 0.6, ω_n = 4.0), check the following constraints: stability: Re(Ξ) = -ζω_n = -2.4 < 0; response time: t_s = 2 / (ζω_n) = 0.83s < 100ms; construct a lookup table: LUT_feasible[(0.6, 4.0)] = {feasible = false, distance_to_boundary = 0.73s}. Create shared memory structure: boundary_info[4][4]=initialized to 0; mutex_matrix[4][4]=create 16 mutex locks; atomic_flags[4]=[IDLE, IDLE, IDLE, IDLE]; version_numbers[4][4]=initialized to 0. Thread 2 (mid-frequency band) executes: acquire read lock mutex_matrix[2][1]; read low-frequency boundary: boundary_info[1][2]=2.0rad / s, version=5; release read lock; calculate new boundary: new_boundary=2.2rad / s; acquire write lock mutex_matrix[2][2]; check version number: version_numbers[2][2]==5; update: boundary_info[2][2]=2.2, version_numbers[2][2]=6; set atomic_flags[2]=UPDATED. A boundary conflict was detected: |boundary[2]-boundary[1]|=|2.2-2.0|=0.2<0.5 (threshold). Priority arbitration was performed: thread 1 priority = 2, thread 2 priority = 1; thread 1's boundary value was retained, and thread 2 rolled back; output: sync_status = {conflict_count = 1, rollback_count = 1, sync_timestamp = 1698234568}. Thread optimization results: thread 0: Ξ_ultralow_config = [-0.1, -0.12]; thread 2: Ξ_mid_config = [-2.4+j3.8, -2.4-j3.8]; a CAS operation was used to update the global optimal solution, outputting x_parallel.

[0152] 4.3. Calculate the convergence error: ||Ξ_mid(k+1)-Ξ_mid(k)||2=0.00006<ε(1e-4); Output the frequency band matching parameter matrix K_freq={K_ultralow, K_low, K_mid, K_high}.

[0153] Step 5: Fusion outputs the coordinated droop parameters.

[0154] 5.1. Calculate the quality indicators of each frequency band: Ultra-low frequency band: ε1=0.0001, σ1=0.002, Q1=1 / (1+0.0001+0.002)=0.998; Mid-frequency band: ε3=0.0003, σ3=0.005, Q3=1 / (1+0.0003+0.005)=0.995; Parameter stability: S1=exp(-||Ξ_final1-Ξ_target1|| / Ξ_target1)=exp(-0.1 / 0.1)=0.368; Overall quality: M1=0.6×Q1+0.4×S1=0.6×0.998+0.4×0.368=0.746; Dynamic weight: w_fusion1=M1 2 / Σ j (M j 2 )=0.746 2 / (0.746 2 +0.682 2 +0.735 2 +0.512 2 )=0.304. Weighted fusion calculation and result verification, perform fusion: R_unified=Σ i w_fusion_i×R_i=0.304×0.02+0.254×0.01+0.289×0.015+0.153×0.005=0.0138; verify physical rationality: 0.01≤0.0138≤1.0; calculate the poles of the fusion system and verify stability.

[0155] 5.2. Verify the droop coefficient range: 0.01 ≤ R_droop ≤ 1.0; Verify the response time: t_response = 45ms ≤ 100ms; Verify the stability margin: PM = 52° ≥ 45°. Considering SOC protection, module 1's SOC = 65%, no adjustment is required. The output can achieve the droop parameter set R_feasible.

[0156] 5.3. Build output format: R_droop_final=[R_module1=0.0142, R_module2=0.0139, R_module3=0.0135, R_module4=0.0137, timestamp=1698234570, checksum=0xA5F3].

[0157] The present invention addresses the problem that existing methods fail to effectively address the coupling of response characteristics of batteries with different chemical systems in multiple frequency bands. Based on a multi-objective hierarchical decomposition method based on modal decoupling theory, the cross-band coupling strength matrix is ​​calculated and a nonlinear decoupling transformation is designed to completely separate the three mutually influencing control objectives of SOC balance, power regulation, and transient suppression into independent frequency bands for optimization. This fundamentally eliminates the contradiction in traditional unified optimization where improved low-frequency band performance leads to deteriorated high-frequency band performance, allowing each frequency band to independently configure optimal parameters based on its own characteristics, effectively avoiding the occurrence of cross-band oscillation. In response to the problem that existing static pole configuration methods are unable to track the evolution of dynamic characteristics caused by factors such as battery aging, temperature changes, and SOC drift, an inverse forgetting factor is adopted and an adaptive forgetting factor is designed to achieve the ideal characteristics of rapid adaptation in the case of large errors and stable tracking in the case of small errors. Combined with the differentiated forgetting factor interval settings of heterogeneous modules, the system can track the dynamic characteristic changes of different types of batteries in real time, ensuring that the control parameters always remain near the optimal operating point, avoiding continuous performance degradation. In order to address the difficulty in pole alignment and coordination caused by inconsistent transfer function orders of heterogeneous modules, a virtual pole supplementation technology is proposed to supplement low-order systems with virtual poles located far away from the working frequency band and assign them extremely low weights. While ensuring the uniformity of mathematical dimensions, the interference of virtual poles on the actual system dynamics is avoided, so that heterogeneous modules of different orders can be coordinated and controlled under a unified optimization framework, completely solving the problem of specific frequency resonance caused by dimensional mismatch in traditional methods. In response to the high computational complexity and lack of parallelization mechanism of existing centralized optimization algorithms, a parallel gradient projection algorithm for pre-calculating the constrained feasible domain is proposed, which converts online constraint projection into a table lookup operation, and adopts a versioned synchronization mechanism to achieve safe parallel optimization of multiple frequency bands. It reduces the optimization calculation time and thread conflict rate, making real-time online optimization of large-scale heterogeneous systems possible and meeting the stringent requirements of millisecond-level control cycle.

[0158] 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 method for adaptively coordinating dynamic droop parameters of heterogeneous battery modules, characterized in that: include: Acquire real-time dynamic response data of heterogeneous battery modules and obtain a standardized module state matrix through normalization processing; Based on this, the transfer function parameters are extracted and the pole distribution matrix is ​​calculated; Based on the pole distribution matrix, the frequency domain characteristics are hierarchically decoupled through the decoupling transformation matrix to obtain a hierarchical optimization objective function set; Based on this, multi-band synchronous solution is performed to obtain the frequency band matching parameter matrix; Based on the frequency band matching parameter matrix, the droop parameter adjustment instruction is obtained through weighted fusion processing; The hierarchical optimization objective function set includes: Based on the pole distribution matrix, the frequency band attribution characteristics of each pole are determined, and the poles are classified into low frequency band, medium frequency band and high frequency band accordingly to obtain the pole frequency band attribution table; Based on the extreme frequency band attribution table, the coupling terms of the preset original unified objective function are analyzed to identify the cross-band correlation matrix C_cross. Based on this, the decoupling transformation matrix T_decouple=I-α×C_cross is constructed to decompose the original unified objective function into a predetermined number of independent sub-objective functions, where I is the identity matrix and α is the decoupling strength coefficient. Based on the sub-objective functions, verify the decoupling effect and construct a hierarchical optimization objective function set; Obtaining the pole distribution matrix includes: Based on the standardized module state matrix, the chemical type of each heterogeneous battery module is identified and the corresponding characteristic coefficient is set. Based on this, the heterogeneous characteristic weighted prediction error is calculated, where the error is weighted using the corresponding characteristic coefficient for modules of different chemical types. Based on the weighted prediction error of heterogeneous characteristics, the adaptive forgetting factor is calculated to obtain a dynamic forgetting factor sequence; based on this, parameter iteration is performed to obtain a dynamic transfer function parameter set and the pole distribution matrix is ​​calculated based on this; The obtained frequency band matching parameter matrix includes: Discretize the continuous constraint space into grid points, adaptively adjust the grid density according to the constraint gradient, perform constraint checks on each grid point, and construct a pre-calculated constraint feasible domain database and spatial index tree. Based on the hierarchical optimization objective function set and the pre-calculated constraint feasible region database, the optimization calculation is performed through parallel gradient projection to obtain the parallel optimization intermediate results; Based on the intermediate results of parallel optimization, the priority negotiation mechanism parameters are adjusted through the thread synchronization mechanism, and the frequency band matching parameter matrix is ​​summarized; The weighted fusion process includes calculating the fusion weights, specifically: Based on the frequency band matching parameter matrix, the convergence error, standard deviation of the solution, and parameter stability of each frequency band are calculated respectively, and a multi-dimensional convergence quality evaluation index is established. Based on this, the comprehensive quality evaluation value of each frequency band is calculated and the weight coefficient of each frequency band in the fusion process is determined. Based on the weight coefficient, a protective weight reduction is implemented for the frequency bands whose comprehensive quality evaluation values ​​are lower than the threshold, and an adjusted fusion weight coefficient matrix is ​​obtained; The standardized module state matrix is ​​obtained including: Based on real-time dynamic response data, it automatically identifies the chemical type of each battery module and distinguishes between lithium iron phosphate modules and ternary lithium modules; According to the chemical type identification results, the matching normalization coefficients of lithium iron phosphate modules and ternary lithium modules are respectively; Based on the matching normalization coefficients, the module data of different chemical types are subjected to differentiated normalization processing to obtain the standardized module state matrix.

2. The method according to claim 1, characterized in that Get the dynamic transfer function parameter set and calculate the pole distribution matrix based on it, including: Based on the dynamic transfer function parameter set, the system order of each heterogeneous battery module is identified, and the module order distribution vector and heterogeneous type identification are obtained; Based on the module order distribution vector, the pole numbers of modules of different orders are aligned, where the low-order modules are supplemented with virtual poles to the highest order to obtain an aligned pole set; Based on the aligned pole set and the heterogeneous type identifier, the actual poles are assigned a first weight and the virtual poles are assigned a second weight to form a standardized pole distribution matrix and a module pole mapping table.

3. The method according to claim 2, characterized in that Determine the frequency band characteristics of each extreme point, including: Based on the standardized pole distribution matrix and the module pole mapping table, the effective frequency of each pole is calculated, where the virtual poles are excluded through weight filtering to obtain a weighted pole frequency matrix; Based on the weighted pole frequency matrix, the frequency band boundaries are dynamically determined by the K-means clustering algorithm, requiring that the interval between adjacent frequency band boundaries is greater than 0.5 decades, and the dynamic frequency band boundary vector is obtained; Based on the weighted pole frequency matrix, dynamic frequency band boundary vector and module pole mapping table, the poles are assigned to the corresponding frequency bands, the frequency band consistency of the conjugate pole pairs is maintained, and a pole frequency band attribution table containing the correspondence between the poles and the modules is obtained.

4. The method according to claim 1, wherein Before extracting the transfer function parameters, the sliding window data preprocessing of the standardized module state matrix is ​​also included, specifically: Based on the standardized module state matrix, the real-time rate of change of the signal dY / dt is calculated. Based on this, the sliding window length is dynamically adjusted according to the relationship of window length L_window = L_base × (1 + β × |dY / dt|), where L_base is the base window length, β is the adjustment coefficient, Y is the output variable of the system, d is the differential operator, and t is time. Based on the dynamically adjusted sliding window length, the standardized module state matrix is ​​segmented to obtain time domain response data segments.

5. The method according to claim 1, wherein The frequency band matching parameter matrix obtained by summarizing includes: Based on the intermediate results of parallel optimization, the minimum spacing between the poles of adjacent frequency bands is checked to identify boundary conflicts caused by the pole spacing being less than the threshold. When a boundary conflict occurs, the priority of frequency band adjustment is determined according to the module importance. Frequency bands with a priority above the threshold remain unchanged, while frequency bands with a priority below the threshold undergo parameter configuration adjustment. Multi-thread synchronous collection is achieved through the exchange of boundary information between threads, and the adjusted optimization results of each frequency band are unified and integrated to obtain the frequency band matching parameter matrix.

Images