Wind storage combined operation optimization method based on energy storage life and frequency modulation performance

Abstract

The application relates to the technical field of energy storage control, and specifically discloses a wind storage combined operation optimization method based on energy storage life and frequency modulation performance, real-time collection of power grid automatic power generation control instructions and preprocessing to obtain a standardized frequency modulation instruction sequence; the instruction sequence is decomposed into a plurality of intrinsic mode function components through variational mode decomposition; Hilbert transformation is performed on each component to obtain an instantaneous frequency and an instantaneous amplitude, the instantaneous frequency is matched with a previously established frequency band-life damage mapping relationship to obtain an instantaneous damage coefficient, and a comprehensive damage coefficient of each component is obtained through weighted average taking the instantaneous amplitude as a weight; the damage change rate of each component is calculated, the comprehensive damage coefficient and the damage change rate are input into a two-dimensional decision table to determine an amplitude attenuation ratio, and the components exceeding a threshold value are compressed point by point; all the components are linearly superimposed to reconstruct a power instruction, and the power instruction is respectively issued to the energy storage and the wind power system after dynamic filtering for execution; the application realizes dynamic balance of the frequency modulation performance and the energy storage life.

Description

Technical Field

[0001] This invention relates to the field of energy storage control technology, and more specifically to an optimization method for wind and energy storage joint operation based on energy storage lifetime and frequency regulation performance. Background Technology

[0002] As the scale of wind power grid connection continues to expand, the volatility and intermittency of wind power output are having an increasingly significant impact on grid frequency stability. To meet grid frequency regulation requirements, wind farms typically configure energy storage systems to form a wind-storage integrated operation system, utilizing the rapid response characteristics of energy storage to participate in primary and secondary frequency regulation services. In actual operation, the energy storage system needs to frequently adjust its charging and discharging according to the grid's automatic generation control commands, and its frequency regulation performance directly affects the performance evaluation and profitability of the wind-storage integrated system.

[0003] Existing wind and energy storage joint operation control technologies generally treat grid frequency regulation commands as a whole time-domain signal for immediate response. They lack a deep understanding of the coupling relationship between different frequency band components within the command and the electrochemical aging mechanism of the energy storage battery. This forces the energy storage system to continue operating in a "high-risk frequency band" that is severely mismatched with its electrochemical impedance spectral characteristics in high-frequency, small-amplitude micro-circulation scenarios. This leads to directional cumulative fatigue damage of charge transfer resistance and diffusion resistance. This frequency-selective, covert aging process is severely out of step with the low marginal returns of the frequency regulation market. Ultimately, the core lifespan resources of the energy storage system are silently overdrawn and consumed while it appears to be operating in compliance with grid requirements. Traditional time-domain response-based control architectures cannot identify and block this "chronic mismatch" process from the frequency domain dimension. Summary of the Invention

[0004] The purpose of this invention is to provide an optimization method for the joint operation of wind and energy storage based on energy storage lifetime and frequency regulation performance, so as to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance includes the following steps: S1: Real-time acquisition of automatic power generation control command signals from the power grid, preprocessing of the command signals to obtain a standardized frequency regulation command sequence; S2: The standardized frequency modulation command sequence is decomposed into multiple eigenmode function components with different center frequencies through variational mode decomposition. S3: Match each intrinsic mode function component with the frequency band-lifetime damage mapping relationship established in advance based on electrochemical impedance spectroscopy test to obtain the lifetime damage coefficient corresponding to each component, and attach the lifetime damage coefficient to the corresponding component to form an intrinsic mode function component set with damage coefficient. S4: Filter the intrinsic mode function component set with damage coefficient, attenuate the amplitude of the component whose lifetime damage coefficient exceeds the preset threshold, and retain the component whose lifetime damage coefficient does not exceed the preset threshold to obtain the adjusted intrinsic mode function component set. S5: Linearly superimpose the adjusted intrinsic mode function component sets to reconstruct the optimized joint operating power command, and then send the joint operating power command to the energy storage system and wind power system for execution.

[0006] As a further aspect of the present invention: S2 specifically includes: A Fourier transform is performed on the standardized frequency modulation command sequence to obtain its spectral distribution. The number of mode decompositions for variational mode decomposition is determined based on the number of peaks in the spectrum whose energy amplitude exceeds a preset energy threshold. The initial center frequency of the variational mode decomposition is initialized with the center frequency corresponding to each peak. The variational optimization problem is solved iteratively using the alternating direction multiplier method. In each iteration, the center frequency and bandwidth of each intrinsic mode function component are updated, and the L2 norm of the difference between the components obtained in two adjacent iterations is calculated. The iteration is terminated when the L2 norm is less than the convergence threshold, and multiple intrinsic mode function components with different center frequencies are output.

[0007] As a further aspect of the present invention: the iterative solution of the variational optimization problem using the alternating direction multiplier method specifically includes: Initialize the Lagrange multipliers and the quadratic penalty factor; substitute the Lagrange multipliers and the quadratic penalty factor into the frequency domain update equation, and sequentially update each eigenmode function component in the frequency domain to obtain the updated eigenmode function components; update each center frequency according to the power spectrum centroid of the updated eigenmode function components to obtain the updated center frequency; update the Lagrange multipliers according to the updated eigenmode function components, the updated center frequency, and the Lagrange multipliers from the previous iteration to obtain the updated Lagrange multipliers; repeat the update process until the relative change of each eigenmode function component is less than the preset convergence threshold, and output the final eigenmode function components and their corresponding center frequencies.

[0008] As a further aspect of the present invention: S3 specifically includes: A Hilbert transform is performed on each intrinsic mode function component to obtain the instantaneous frequency sequence of each component. Each frequency value in the instantaneous frequency sequence is matched with a pre-constructed frequency band-lifetime damage mapping table to obtain the instantaneous damage coefficient corresponding to each frequency point. A weighted average is performed on all instantaneous damage coefficients within the same intrinsic mode function component, using the amplitude of the corresponding frequency point in the instantaneous frequency sequence as the weight, to calculate the comprehensive damage coefficient of the corresponding component. The comprehensive damage coefficient is assigned to the corresponding intrinsic mode function component as an additional attribute to form an intrinsic mode function component set with damage coefficients.

[0009] As a further aspect of the present invention: the calculation process of the instantaneous damage coefficient is as follows: Equivalent circuit fitting was performed on the impedance spectrum data obtained from electrochemical impedance spectroscopy to extract the characteristic curves of charge transfer resistance and diffusion resistance as a function of frequency. Based on the inflection point frequency of the sudden increase in internal resistance value in the characteristic curve, the frequency range was divided into a low-frequency diffusion region, a mid-frequency charge transfer region, and a high-frequency ohmic region, and a benchmark damage coefficient was set for each frequency range. Each frequency value in the instantaneous frequency sequence was sequentially assigned to its corresponding frequency range, and the benchmark damage coefficient corresponding to the frequency range was extracted. Based on the specific position of the frequency value within its range, linear interpolation was performed using the rate of change of internal resistance corresponding to the two endpoints of the frequency range to calculate the instantaneous damage coefficient corresponding to the frequency value.

[0010] As a further aspect of the present invention: S4 specifically includes: Calculate the rate of change of the damage coefficient over time in each intrinsic mode function component with a damage coefficient to obtain the damage change rate of each component. Input the damage coefficient and damage change rate of each component into a preset two-dimensional decision table. The two-dimensional decision table is divided into multiple decision regions with the damage coefficient as the horizontal axis and the damage change rate as the vertical axis. Each decision region corresponds to a different amplitude attenuation ratio. According to the landing region of each component in the two-dimensional decision table, read the corresponding amplitude attenuation ratio. Compress the instantaneous amplitude sequence of each component point by point according to the amplitude attenuation ratio to obtain the attenuated intrinsic mode function components. Components with damage coefficients not exceeding a preset threshold are retained as intrinsic mode function components. All attenuated components and retained components are combined to form the adjusted intrinsic mode function component set.

[0011] As a further aspect of the present invention: the step of compressing the instantaneous amplitude sequence of each component point by point according to the amplitude attenuation ratio to obtain the attenuated intrinsic mode function components specifically includes: Hilbert transform is performed on each attenuated intrinsic mode function component to obtain the instantaneous phase sequence of each component. The instantaneous amplitude sequence of each component is divided into continuous half-wave segments according to the zero-crossing point of the instantaneous phase sequence. The peak value of the instantaneous amplitude in each half-wave segment is calculated, and the peak value is multiplied by the corresponding amplitude attenuation ratio to obtain the target peak value of the corresponding half-wave segment. With the target peak value as the upper limit, all instantaneous amplitudes in the half-wave segment are proportionally limited and compressed so that the maximum amplitude of the compressed half-wave segment is equal to the target peak value. All compressed half-wave segments are spliced ​​together in their original time sequence to form the attenuated intrinsic mode function components.

[0012] As a further aspect of the present invention: S5 specifically includes: The initial reconfiguration power command is obtained by summing all components in the adjusted intrinsic mode function component set time-by-time. The initial reconfiguration power command is then input into a finite impulse response low-pass filter to extract the low-frequency trend component as the wind power system execution command. The low-frequency trend component is subtracted from the initial reconfiguration power command to obtain the high-frequency compensation component as the energy storage system execution command. The current state of charge of the energy storage system and the current available regulating capacity of the wind power system are collected. The energy storage system execution command is then limited and corrected based on the difference between the energy storage system execution command and the high-frequency compensation component. The corrected energy storage system execution command and the wind power system execution command are then sent to the corresponding systems for execution.

[0013] As a further aspect of the present invention: the step of accumulating all components in the adjusted intrinsic mode function component set time-by-time to obtain the initial reconfiguration power command specifically includes: The rate of change of the current output power of the wind power system is obtained, and the rate of change is input into a preset cutoff frequency mapping function. The mapping function outputs the current cutoff frequency according to the rule that the cutoff frequency increases with the rate of change. The filter coefficient sequence of the finite impulse response low-pass filter is calculated based on the current cutoff frequency. The sum of the coefficients in the filter coefficient sequence is 1 and they are symmetrically distributed. The initial reconstructed power command is convolved with the filter coefficient sequence to obtain a smoothed power sequence after moving average. The smoothed power sequence is output as a low-frequency trend component, which is used as the wind power system execution command.

[0014] The beneficial effects of this invention are: (1) This invention decomposes the grid frequency regulation command into intrinsic mode function components of different frequency bands through variational mode decomposition, and quantifies the damage coefficient of each component by combining the frequency band-lifetime damage mapping relationship established by electrochemical impedance spectroscopy testing. On this basis, a two-dimensional decision table is introduced to compress the amplitude of components with high damage coefficients point by point based on half-wave segments, so that the energy storage system can effectively avoid the "high-risk frequency band" response that causes serious battery damage in high-frequency, small-amplitude fluctuation scenarios. Compared with traditional instantaneous response control, this method avoids the ineffective wear of energy storage batteries in low-yield frequency regulation tasks, reduces the cumulative damage of battery cycles, thereby extending the actual service life of the energy storage system and reducing the replacement and operation and maintenance costs of the wind-storage combined system throughout its entire life cycle.

[0015] (2) While delaying the degradation of energy storage lifespan, this invention reconstructs the frequency domain of the frequency regulation command in a refined manner, linearly superimposing the attenuated component with the retained component, and dynamically adjusting the cutoff frequency of the low-pass filter according to the real-time power change rate of the wind power system to extract the low-frequency trend component suitable for wind power tracking. The energy storage system then undertakes the high-frequency compensation component. This differentiated allocation strategy fully leverages the complementary characteristics of wind power and energy storage, while ensuring the rapid response capability and tracking accuracy of the combined system to the grid frequency regulation command. Furthermore, by combining the limiting correction of the energy storage state of charge and the available regulating capacity of wind power, it avoids command execution deviations caused by insufficient energy storage capacity. While meeting the grid frequency regulation assessment requirements, it achieves a dynamic balance between frequency regulation benefits and lifespan loss, thereby improving the overall economic efficiency of wind and energy storage joint operation. Attached Figure Description

[0016] The invention will now be further described with reference to the accompanying drawings.

[0017] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

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

[0019] Please see Figure 1 As shown, this invention is an optimization method for the joint operation of wind and energy storage based on energy storage lifetime and frequency regulation performance, including the following steps: S1: Real-time acquisition of automatic power generation control command signals from the power grid, preprocessing of the command signals to obtain a standardized frequency regulation command sequence; S2: The standardized frequency modulation command sequence is decomposed into multiple eigenmode function components with different center frequencies through variational mode decomposition. S3: Match each intrinsic mode function component with the frequency band-lifetime damage mapping relationship established in advance based on electrochemical impedance spectroscopy test to obtain the lifetime damage coefficient corresponding to each component, and attach the lifetime damage coefficient to the corresponding component to form an intrinsic mode function component set with damage coefficient. S4: Filter the intrinsic mode function component set with damage coefficient, attenuate the amplitude of the component whose lifetime damage coefficient exceeds the preset threshold, and retain the component whose lifetime damage coefficient does not exceed the preset threshold to obtain the adjusted intrinsic mode function component set. S5: Linearly superimpose the adjusted intrinsic mode function component sets to reconstruct the optimized joint operating power command, and then send the joint operating power command to the energy storage system and wind power system for execution.

[0020] In S1, the automatic generation control command signal of the power grid is acquired in real time, and the command signal is preprocessed to obtain a standardized frequency regulation command sequence, specifically including: First, the automatic generation control command signals issued by the power grid dispatching terminal are acquired in real time through a remote communication device. This remote communication device is connected to the dispatching data network and uses the IEC104 communication protocol to receive command data at a rate of 1 point per second, and stores the received command signals in the buffer as floating-point numbers.

[0021] Secondly, outlier removal is performed on the original instruction signals in the buffer. The arithmetic mean and standard deviation of all instruction signals within the previous 60 seconds are calculated. The absolute value of the difference between the current instruction signal and the arithmetic mean is compared with three times the standard deviation. If the absolute value is greater than three times the standard deviation, the current instruction signal is determined to be an outlier, and the arithmetic mean is used to replace the outlier to obtain the corrected instruction signal.

[0022] Next, the corrected command signal is normalized. All corrected command signals stored in the historical database for the past 30 days are traversed, and the maximum and minimum values ​​are extracted. For the corrected command signal at the current moment, the difference between it and the minimum value is calculated, and then the difference between the maximum and the minimum value is calculated. The ratio of these two differences is used as the standardized command value for the current moment. The range of this standardized command value is constrained to between -1 and 1.

[0023] Finally, the 60 standardized instruction values ​​obtained within 60 consecutive seconds are arranged in chronological order to form a standardized frequency modulation instruction sequence containing 60 data points, and this sequence is output to memory for use in the next step.

[0024] In S2, the standardized frequency modulation command sequence is decomposed into multiple eigenmode function components with different center frequencies through variational mode decomposition, specifically including: First, a Fourier transform is performed on the standardized frequency modulation (FM) command sequence. Using the 60-point FM command sequence as input, a Fast Fourier Transform (FFT) algorithm is applied to obtain the real and imaginary parts of the sequence in the frequency domain. For each frequency point, the sum of the squares of the real and imaginary parts is calculated, and the square root of this sum is taken to obtain the energy amplitude corresponding to that frequency point. All frequency points are arranged in ascending order of frequency, and together with their corresponding energy amplitudes, constitute the spectral distribution of the FM command sequence.

[0025] Secondly, determine the number of mode decompositions in the variational mode decomposition. Traverse all frequency points in the spectral distribution, comparing the energy amplitude of each frequency point with the energy amplitudes of its two preceding and following frequency points. If the energy amplitude of a frequency point is greater than both its preceding and following frequencies, it is identified as a peak, and its frequency and energy amplitude are recorded. A preset energy threshold is set to 5% of the maximum energy amplitude. Count the number of peaks whose energy amplitude exceeds this preset energy threshold, and use this number as the number of mode decompositions in the variational mode decomposition.

[0026] Next, the center frequencies of each intrinsic mode function component are initialized. The frequency values ​​corresponding to each selected peak are arranged in ascending order and used as the initial center frequencies of the 1st to Nth intrinsic mode function components, where N is the number of mode decompositions.

[0027] Then, the variational optimization problem is solved iteratively using the alternating direction multiplier method. The Lagrange multipliers are initialized to a zero sequence of the same length as the normalized frequency modulation command sequence, and the quadratic penalty factor is initialized to 2000. The iteration loop is then entered, and the following update steps are performed in each iteration: The first step is to update the frequency domain of each intrinsic mode function (IMF) component sequentially. For the k-th IMF component, first calculate the Fourier transform of the normalized frequency modulation command sequence to obtain its frequency domain representation. Then, subtract the current frequency domain estimates of all components except the k-th component from this frequency domain representation, and add the current frequency domain value of the Lagrange multiplier divided by the quadratic penalty factor to obtain an intermediate variable. Divide the intermediate variable by 1, add 2, multiply by the quadratic penalty factor, and multiply by the square of the difference between the frequency and the current center frequency of the k-th component to obtain the updated frequency domain value of the k-th component. Perform the above operation sequentially for all k components to complete the frequency domain update of all IMF components.

[0028] The second step is to update each center frequency. For the k-th intrinsic mode function component, take the frequency domain value of this component obtained after the update in the first step, and calculate its power spectral density from 0 to positive infinity. Treat the power spectral density of this component as a weighting function of the frequency, and calculate the weighted average frequency of the power spectral density. Specifically, multiply the frequency by the integral of the power spectral density over the range of 0 to positive infinity, and divide by the integral of the power spectral density over the range of 0 to positive infinity. The resulting value is used as the updated center frequency of the k-th component. Perform the above operation sequentially for all k components to complete the update of all center frequencies.

[0029] The third step is to update the Lagrange multipliers. Take the frequency domain representation of the normalized frequency modulation command sequence, subtract the sum of the frequency domain values ​​of all intrinsic mode function components obtained after the first step update, and obtain the frequency domain representation of the residual. Add the current frequency domain value of the Lagrange multipliers to the quadratic penalty factor multiplied by the frequency domain representation of the residual to obtain the updated frequency domain representation of the Lagrange multipliers.

[0030] After the update is completed, the convergence condition of the iteration is determined. For each intrinsic mode function component, the L2 norm of the difference between its updated time-domain value and its time-domain value in the previous iteration is calculated, and then divided by the L2 norm of its time-domain value in the previous iteration to obtain the relative change of the component. When the relative changes of all components are less than the preset convergence threshold of 0.001, the iteration is terminated; otherwise, the frequency domain values, center frequencies, and Lagrange multipliers of each component updated in the current iteration are used as initial values ​​to enter the next iteration.

[0031] Finally, when the iteration terminates, the frequency domain values ​​of each eigenmode function component obtained from the final update are subjected to inverse Fourier transform to obtain the time domain sequence corresponding to each component. The time domain sequence corresponding to each component is output together with the final center frequency of each component as multiple eigenmode function components with different center frequencies.

[0032] In S3, each intrinsic mode function component is matched with a frequency band-lifetime damage mapping relationship pre-established based on electrochemical impedance spectroscopy testing to obtain the lifetime damage coefficient corresponding to each component. The lifetime damage coefficient is then appended to the corresponding component to form an intrinsic mode function component set with damage coefficients, specifically including: First, a Hilbert transform is performed on each intrinsic mode function (IMF) component to obtain the instantaneous frequency and amplitude sequences of each component. For each IMF component, it is treated as the real part signal, and its analytic signal is constructed using a Hilbert transform. Specifically, a Fourier transform is performed on the component, doubling the amplitude of the positive frequency portion and setting the negative frequency portion to zero, followed by an inverse Fourier transform to obtain the imaginary part signal. After constructing the analytic signal from the real and imaginary parts, the phase angle of the analytic signal at each moment is calculated. The rate of change of this phase angle with time is the instantaneous frequency at that moment. Simultaneously, the square root of the sum of the squares of the real and imaginary parts at each moment is calculated to obtain the instantaneous amplitude at that moment. The instantaneous frequencies at all moments are arranged in chronological order to obtain the instantaneous frequency sequence of the component; the instantaneous amplitudes at all moments are arranged in chronological order to obtain the instantaneous amplitude sequence of the component. The above operations are performed sequentially on all IMF components to obtain the instantaneous frequency and amplitude sequences corresponding to each component.

[0033] Secondly, a frequency band-lifetime damage mapping table was pre-constructed. This mapping table was established through electrochemical impedance spectroscopy (EIS) testing and cyclic aging experiments on energy storage batteries. The specific process is as follows: Sample batteries of the same model as those used in the wind-storage combined system were selected. Under a constant temperature environment of 25 degrees Celsius, frequency scanning tests were performed on the batteries using an electrochemical workstation. The scanning frequency range was set from 0.01 Hz to 10000 Hz. A sinusoidal AC excitation signal with an amplitude of 5 mV was applied at each frequency point, and the real and imaginary parts of the battery impedance were recorded to obtain impedance spectrum data at different frequencies. The obtained impedance spectrum data were fitted with an equivalent circuit. The equivalent circuit consisted of an ohmic internal resistance, a charge transfer internal resistance, and a diffusion impedance element connected in series. The parameter values ​​of each element were obtained by fitting using the complex nonlinear least squares method, and the characteristic curves of the charge transfer internal resistance and diffusion internal resistance as a function of frequency were extracted. Based on the inflection point frequency at which the internal resistance value increases abruptly with frequency in the characteristic curve, the frequency range is divided into three intervals: the low-frequency diffusion region (frequency below 0.1 Hz), the mid-frequency charge transfer region (frequency between 0.1 Hz and 100 Hz), and the high-frequency ohmic region (frequency above 100 Hz). Subsequently, constant-frequency charge-discharge cycle aging tests were conducted on sample batteries at typical frequency points within these three frequency intervals. The capacity decay rate per cycle was measured at each frequency point, and the capacity decay rate was used as a measure of the degree of damage. The average capacity decay rate of all test points within each frequency interval was normalized and used as the baseline damage coefficient for that interval. Specifically, the baseline damage coefficient for the low-frequency diffusion region was set to 0.8, the baseline damage coefficient for the mid-frequency charge transfer region was set to 1.0, and the baseline damage coefficient for the high-frequency ohmic region was set to 0.5. The frequency interval divisions and corresponding baseline damage coefficients were stored in a lookup table to construct a frequency band-lifetime damage mapping table.

[0034] Next, each frequency value in the instantaneous frequency sequence of each intrinsic mode function component is sequentially matched with the frequency band-lifetime damage mapping table to obtain the instantaneous damage coefficient corresponding to each frequency point. The matching process is as follows: For a certain frequency value in the instantaneous frequency sequence... First, determine the frequency range to which it belongs. If If the frequency is less than 0.1 Hz, it belongs to the low-frequency diffusion region, and the lower limit frequency of this region is denoted as ______. The upper limit frequency is 0.01 Hz. The corresponding baseline damage coefficient is 0.1 Hz. The baseline damage factor for the low-frequency diffusion region is 0.8. The baseline damage coefficient for adjacent intervals (mid-frequency charge transfer region) is 1.0; if The frequency range between 0.1 Hz and 100 Hz falls within the mid-frequency charge transfer region, denoted as... It is 0.1 Hz. 100 Hz It is 1.0. The baseline damage coefficient for the high-frequency ohmic region is 0.5; if Frequency values ​​greater than 100 Hz fall within the high-frequency ohmic region. However, since the extrapolation of the high-frequency region is undefined, a baseline damage coefficient of 0.5 is directly used. For frequency points falling within the low-frequency diffusion region and the mid-frequency charge transfer region, a linear interpolation formula is used to calculate the instantaneous damage coefficient. The calculation formula is as follows: ; in, This is the current frequency value. and These are the lower and upper frequency limits of the frequency range in which the frequency falls, respectively. and These are the reference damage coefficients corresponding to the lower and upper frequency limits of the range, respectively. For frequency points located in the high-frequency ohmic region, we directly take... It equals 0.5. Perform the above operation on each frequency point in all instantaneous frequency sequences to obtain the instantaneous damage coefficient sequence corresponding to that component, whose length is the same as that of the instantaneous frequency sequence.

[0035] Then, a weighted average is performed on all instantaneous damage coefficients within the same intrinsic mode function component, using the amplitude at the corresponding moment in the instantaneous amplitude sequence as the weight, to calculate the comprehensive damage coefficient of that component. Let the instantaneous amplitude sequence of the k-th intrinsic mode function component be... , ,..., Where N is the sequence length, and the corresponding instantaneous damage coefficient sequence is: , ,..., Then the overall damage coefficient of this component. Calculated using the weighted average formula: ; In the formula, Indicates the sequence number of the time sampling point. The k-th eigenmode function component represents the component at the k-th eigenmode function. The instantaneous amplitude at each moment, The k-th eigenmode function component represents the component at the k-th eigenmode function. The instantaneous damage coefficient corresponding to each moment is calculated by multiplying the instantaneous amplitude at each moment by the corresponding instantaneous damage coefficient in the numerator and the sum of the instantaneous amplitudes at each moment in the denominator. The resulting ratio is the comprehensive damage coefficient of that component. This coefficient reflects the average degree of damage to battery life caused by that intrinsic mode function component throughout the entire time history.

[0036] Finally, the calculated comprehensive damage coefficient The corresponding k-th intrinsic mode function component is assigned an additional attribute, that is, a field is added to the data storage structure of this component to record its comprehensive damage coefficient. After performing the above operation on all intrinsic mode function components, a set of intrinsic mode function components labeled with damage coefficients is obtained. These components are then aggregated to form a set of intrinsic mode function components with damage coefficients for use in subsequent steps.

[0037] In S4, the intrinsic mode function component set with damage coefficients is filtered. Components with lifetime damage coefficients exceeding a preset threshold are attenuated in amplitude, while components with lifetime damage coefficients below the preset threshold are retained, resulting in an adjusted intrinsic mode function component set, specifically including: First, the rate of change of the damage coefficient over time is calculated for each intrinsic mode function component with a damage coefficient, thus obtaining the damage change rate of each component. For each intrinsic mode function component with a damage coefficient, its instantaneous damage coefficient sequence is obtained, which contains the damage coefficient values ​​at N time points from the start time to the end time. Starting from the second time point, the difference between the damage coefficient at the current time and the damage coefficient at the previous time is calculated sequentially, and then divided by the sampling time interval to obtain the instantaneous change rate at the current time. The sampling time interval is 0.0167 seconds, corresponding to 60 sampling points per second. The arithmetic mean of the instantaneous change rates at all times is calculated, and this average value is taken as the damage change rate of that component. This damage change rate reflects the trend of the component's contribution to the aggravation or mitigation of battery life damage; a positive value indicates an aggravating trend of damage, and a negative value indicates a mitigating trend of damage.

[0038] Secondly, a two-dimensional decision table is pre-constructed. This table uses the damage coefficient as the x-axis and the damage change rate as the y-axis, dividing the two-dimensional plane into multiple decision regions, each corresponding to a specific amplitude attenuation ratio. The specific construction process is as follows: the damage coefficient, ranging from 0 to 2, is uniformly divided into 10 intervals, each with a width of 0.2; the damage change rate, ranging from -0.05 to +0.05 per second, is uniformly divided into 10 intervals, each with a width of 0.01 per second. This forms 100 rectangular decision regions, each 10 x 10. Through offline simulation experiments, different amplitude attenuation ratios are applied to the intrinsic mode function components falling into each decision region. The optimal amplitude attenuation ratio for each decision region is determined with the goal of achieving the optimal comprehensive weighted index of energy storage lifetime extension and frequency regulation performance loss. Each decision region and its corresponding optimal amplitude attenuation ratio are stored in a two-dimensional lookup table, forming the pre-constructed two-dimensional decision table. In this decision table, the region with the larger the damage coefficient and the greater the rate of damage change corresponds to a higher amplitude attenuation ratio, with values ​​ranging from 0.2 to 0.9.

[0039] Next, based on the landing area of ​​each component in the two-dimensional decision table, the corresponding amplitude attenuation ratio is read. For each intrinsic mode function component with an impairment coefficient, its comprehensive impairment coefficient is used as the abscissa value, and its impairment change rate is used as the ordinate value. The specific rectangular area where these two coordinate values ​​fall in the two-dimensional decision table is located. The amplitude attenuation ratio corresponding to this rectangular area is found and used as the amplitude attenuation ratio of that component. If the comprehensive impairment coefficient of the component is less than the preset impairment coefficient threshold of 0.6, the component is determined to be a low-impairment component, and no amplitude attenuation processing is performed. It is directly marked as a retained intrinsic mode function component. If the comprehensive impairment coefficient is greater than or equal to 0.6, the component is determined to be a high-impairment component, and amplitude attenuation is required. The attenuation ratio obtained from the lookup table is used as the basis for subsequent compression processing.

[0040] Then, the instantaneous amplitude sequences of each high-damage component are compressed point-by-point according to the amplitude attenuation ratio to obtain the attenuated intrinsic mode function components. For each intrinsic mode function component that needs amplitude attenuation, a Hilbert transform is first performed to obtain the instantaneous phase sequence of that component. The specific method of the Hilbert transform is the same as in step S3: a Fourier transform is performed on the component, doubling the amplitude of the positive frequency part and setting the negative frequency part to zero, followed by an inverse Fourier transform to obtain the imaginary part signal of the component. After constructing an analytic signal from the real part signal and the imaginary part signal, the phase angle of the analytic signal at each moment is calculated to obtain the instantaneous phase sequence.

[0041] Next, the instantaneous amplitude sequence of each component is divided into continuous half-wave segments based on the zero-crossing points of the instantaneous phase sequence. Traversing the instantaneous phase sequence, starting from the initial moment, each time a phase value crosses zero from a negative value to a positive value, or from a positive value to a negative value, that moment is recorded as a zero-crossing point. The time interval between two adjacent zero-crossing points constitutes a half-wave segment. Segmenting the instantaneous amplitude sequence according to these zero-crossing points yields several continuous half-wave segments, each corresponding to the amplitude envelope of a positive or negative half-wave.

[0042] Then, the peak value of the instantaneous amplitude within each half-wave segment is calculated, which is the maximum value among all instantaneous amplitudes within that segment. This peak value is multiplied by the amplitude attenuation ratio corresponding to that component to obtain the target peak value for that half-wave segment. Using this target peak value as the upper limit, all instantaneous amplitudes within the half-wave segment are proportionally limited and compressed. The specific compression method is as follows: for the instantaneous amplitude at each moment within the half-wave segment, its ratio to the original peak value is calculated, and then multiplied by the target peak value to obtain the compressed instantaneous amplitude value. After this processing, the maximum instantaneous amplitude within the half-wave segment is exactly equal to the target peak value, while the relative proportional relationship of the amplitudes at different moments within the segment remains unchanged.

[0043] Finally, all compressed half-wave segments are sequentially spliced ​​together in their original time order to form a complete instantaneous amplitude sequence. This sequence is then combined with the instantaneous phase information of the original component to reconstruct the attenuated intrinsic mode function (IMF) component time-domain signal. The above compression operation is performed on all high-impairment components requiring attenuation to obtain the corresponding attenuated IMF components. All attenuated IMF components are then combined with the previously marked low-impairment IMF components to form an adjusted IMF component set for use in subsequent steps.

[0044] In S5, the adjusted intrinsic mode function component sets are linearly superimposed to reconstruct the optimized joint operating power command, which is then issued to the energy storage system and wind power system for execution. Specifically, this includes: First, the initial reconfiguration power command is obtained by accumulating all components in the adjusted intrinsic mode function (IMF) component set time-by-time. The specific process is as follows: The adjusted IMF component set is obtained, containing M components, each a time-domain sequence with 60 data points. Starting from time 1, the values ​​of all M components at that time are summed to obtain the accumulated value for time 1. The same accumulation operation is performed sequentially from time 2 to time 60, resulting in accumulated values ​​for all 60 times. These 60 accumulated values ​​are arranged in chronological order to form the initial reconfiguration power command sequence, which has a length of 60 and uses the same unit as the original AGC command, namely megawatts.

[0045] Secondly, a finite impulse response (FIR) low-pass filter is constructed and low-frequency trend components are extracted. The cutoff frequency of this filter is dynamically adjusted according to the rate of change of the current output power of the wind power system. The specific implementation steps are as follows: The first step is to obtain the rate of change of the current output power of the wind power system. The active power output value of the wind turbines at the current moment is collected in real time through the wind farm's remote control terminal, while the active power output value of the previous moment is read from the historical database. The sampling interval is 1 second. The rate of change of power at the current moment is calculated using the formula: current power value minus previous power value, then divided by the sampling interval of 1 second, yielding the rate of change of power in megawatts per second. The absolute value of this rate of change is used as the input for subsequent calculations.

[0046] The second step involves inputting the power change rate into a preset cutoff frequency mapping function to obtain the current cutoff frequency. This mapping function is pre-established offline. Specifically, the maximum trackable frequency of the wind power system under different power change rates is determined through simulation experiments. The power change rate range from 0 to 10 MW / s is divided into 20 intervals, each with a width of 0.5 MW / s. A corresponding cutoff frequency value is set for each interval, with the cutoff frequency ranging from 0.01 Hz to 0.5 Hz. The above correspondence is stored in the form of a lookup table. During actual operation, based on the absolute value of the current power change rate, the interval to which it belongs is looked up, and the corresponding cutoff frequency value for that interval is read as the current filter's cutoff frequency. This cutoff frequency controls the filter's passband width; the larger the power change rate, the higher the cutoff frequency, allowing more rapidly changing low-frequency components to pass through.

[0047] The third step is to calculate the filter coefficient sequence of the finite impulse response (FIR) low-pass filter based on the current cutoff frequency. The filter order is fixed at 20, meaning the filter coefficient sequence contains 21 coefficients. First, the normalized cutoff frequency is calculated based on the current cutoff frequency. The normalized cutoff frequency is equal to the current cutoff frequency divided by the Nyquist frequency, where the Nyquist frequency is half the sampling frequency (1 Hz), hence 0.5 Hz. Then, the filter coefficients are designed using the window function method: using a Hamming window as the window function, the unit impulse response of the ideal low-pass filter is calculated. The unit impulse response of the ideal low-pass filter is the Singer function centered at time 0, expressed as 2 multiplied by the normalized cutoff frequency multiplied by the Singer function. The independent variable of the Singer function is 2 multiplied by the normalized cutoff frequency multiplied by the integer index. For the 21 integer points from -10 to +10, the ideal unit impulse response value is calculated, and then multiplied by the corresponding Hamming window function value to obtain 21 unnormalized filter coefficients. Finally, these 21 coefficients are normalized by dividing each coefficient by the sum of all 21 coefficients, so that the sum of all normalized coefficients equals 1. Furthermore, since the window function is even-symmetric and the unit impulse response of the ideal low-pass filter is also even-symmetric, the final filter coefficient sequence exhibits a symmetrical distribution; that is, the 1st coefficient is equal to the 21st coefficient, the 2nd coefficient is equal to the 20th coefficient, and so on.

[0048] The fourth step involves convolving the initial reconstructed power command with the filter coefficient sequence to obtain a smoothed power sequence after moving average. The specific process of convolution is as follows: For the j-th time point in the initial reconstructed power command sequence, take this point and the 20 time points preceding it (padding with zeros if necessary) to form a segment of length 21. Multiply this segment with the filter coefficient sequence at corresponding positions and sum them to obtain the convolution result at the j-th time point. This operation is performed sequentially from time point 1 to time point 60, resulting in 60 convolution results, which, arranged in chronological order, constitute the smoothed power sequence. This smoothed power sequence filters out high-frequency fluctuations in the original command while retaining low-frequency trend components that match the response capability of the wind power system.

[0049] The smoothed power sequence is output as a low-frequency trend component, which serves as the execution command of the wind power system and is sent to the power control system of the wind turbine, whereby the wind turbine tracks and executes it.

[0050] Next, the energy storage system executes the command. The initial reconfiguration power command is subtracted from the aforementioned low-frequency trend component time-by-time, yielding the difference at each moment. This difference sequence is the high-frequency compensation component. This component contains the high-frequency fluctuations filtered out from the original command and needs to be handled by the energy storage system with a faster response speed. This high-frequency compensation component is used as the initial value for the energy storage system's execution command.

[0051] Then, the commands executed by the energy storage system are limited and corrected. The current state of charge of the energy storage system is collected in real time, that is, the percentage of the battery's remaining capacity relative to its rated capacity, with a value ranging from 0 to 100%. At the same time, the current available adjustable capacity of the wind power system is collected, that is, the power that the wind turbine can increase or decrease from its current output, in megawatts. Based on the sign and magnitude of the energy storage system's execution command, and considering the current state of charge (SBC) and available wind power regulation capacity, the following adjustments are made: If the energy storage system's execution command is positive (indicating a need for energy storage discharge) and the current SBC is below 20%, the energy storage discharge capacity is deemed insufficient, and the command value is reduced to 50% of its original value. This reduction is then added to the wind power system's execution command, provided the wind power system has sufficient upward regulation capacity. If the energy storage system's execution command is negative (indicating a need for energy storage charging) and the current SBC is above 80%, the energy storage charging capacity is deemed insufficient, and the command value is reduced to 50% of its original value. This reduction is then added to the wind power system's execution command, provided the wind power system has sufficient downward regulation capacity. If the absolute value of the energy storage system's execution command exceeds the rated power of the energy storage converter, it is limited to the rated power value. After these adjustments, the final energy storage system execution command is obtained.

[0052] Finally, the revised energy storage system execution command and the aforementioned wind power system execution command are sent to the energy storage converter and the wind turbine power controller respectively through their respective remote communication devices. The energy storage system and the wind power system then coordinate to perform joint operation power regulation according to the command values.

[0053] In S5, the adjusted intrinsic mode function component sets are linearly superimposed to reconstruct the optimized joint operating power command, which is then issued to the energy storage system and wind power system for execution. Specifically, this includes: First, the initial reconfiguration power command is obtained by accumulating all components in the adjusted intrinsic mode function (IMF) component set time-by-time. The specific process is as follows: The adjusted IMF component set is obtained, containing M components, each a time-domain sequence with 60 time-domain data points. Starting from time 1, the values ​​of all M components at that time are summed to obtain the accumulated value for time 1. The same accumulation operation is performed sequentially from time 2 to time 60, resulting in accumulated values ​​for all 60 time points. These 60 accumulated values ​​are arranged in chronological order to form the initial reconfiguration power command sequence. This sequence has a length of 60, and the unit is the same as the original automatic generation control command, in megawatts (MW).

[0054] Secondly, a finite impulse response (FIR) low-pass filter is constructed. The parameters of this filter are dynamically determined based on the real-time operating status of the wind power system. This filter is used to extract low-frequency trend components suitable for wind power system tracking from the initial reconfiguration power command. Specifically, this includes the following sub-steps: The first step is to obtain the rate of change of the current output power of the wind power system. The active power output value of the wind turbines at the current moment is collected in real time through the wind farm's remote control terminal, while the active power output value of the previous moment is read from the historical database. The sampling interval is 1 second. The rate of change of power at the current moment is calculated, specifically by subtracting the previous moment's power value from the current moment's power value, and then dividing by the sampling interval of 1 second, yielding the rate of change of power in megawatts per second. The absolute value of this rate of change is used as the input parameter for subsequent calculations.

[0055] The second step involves inputting the power change rate into a preset cutoff frequency mapping function to obtain the current cutoff frequency. This mapping function is pre-established offline. Specifically, it is constructed by determining the maximum trackable frequency of the wind power system under different power change rates through simulation experiments. The power change rate range from 0 to 10 MW / s is divided into 20 intervals, each with a width of 0.5 MW / s. A corresponding cutoff frequency value is set for each interval, ranging from 0.01 Hz to 0.5 Hz. This correspondence is stored in local memory as a lookup table. During actual operation, based on the absolute value of the current power change rate, the corresponding interval is looked up, and the cutoff frequency value for that interval is read as the current filter's cutoff frequency. This cutoff frequency controls the filter's passband width. Its setting logic is that the higher the power change rate, the higher the cutoff frequency, allowing more rapidly changing low-frequency components to pass, thus meeting the requirement of tracking a wider frequency band when the wind power system experiences severe fluctuations.

[0056] The third step is to calculate the filter coefficient sequence of the finite impulse response (FIR) low-pass filter based on the current cutoff frequency. The filter order is fixed at 20, meaning the filter coefficient sequence contains 21 coefficients. First, the normalized cutoff frequency is calculated based on the current cutoff frequency, which is equal to the current cutoff frequency divided by the Nyquist frequency. The Nyquist frequency is half the sampling frequency, which is 1 Hz, so the Nyquist frequency is 0.5 Hz. Then, the filter coefficients are designed using the window function method: using a Hamming window as the window function, the unit impulse response of the ideal low-pass filter is calculated. For 21 integer points from -10 to +10, the ideal unit impulse response value is calculated, and then multiplied by the corresponding Hamming window function value to obtain 21 unnormalized filter coefficients. Finally, these 21 coefficients are normalized by dividing each coefficient by the sum of all 21 coefficients, so that the sum of all normalized coefficients equals 1. Meanwhile, since the Hamming window is even-symmetric and the unit impulse response of the ideal low-pass filter is also even-symmetric, the final filter coefficient sequence is symmetrically distributed, that is, the first coefficient is equal to the 21st coefficient, the second coefficient is equal to the 20th coefficient, and so on.

[0057] The fourth step involves convolving the initial reconstructed power command with the filter coefficient sequence to obtain a smoothed power sequence after moving average. The specific process of convolution is as follows: For the j-th time point in the initial reconstructed power command sequence, this point and the 20 time points preceding it form a segment of length 21. If the number of points before the current time point is less than 20, it is padded with zeros to a length of 21. This segment is multiplied at corresponding positions by the filter coefficient sequence and accumulated to obtain the convolution result at the j-th time point. This operation is performed sequentially from time point 1 to time point 60, resulting in 60 convolution results, which, arranged in chronological order, constitute the smoothed power sequence. This smoothed power sequence filters out high-frequency fluctuations in the original command while retaining low-frequency trend components that match the response capability of the wind power system.

[0058] The smoothed power sequence is output as a low-frequency trend component, which serves as the execution command of the wind power system. This command is sent to the power controller of the wind turbine through a remote communication device, and the wind turbine tracks and executes it.

[0059] Next, the energy storage system executes the command. The initial reconfiguration power command is subtracted from the aforementioned low-frequency trend component time-by-time, yielding the difference at each moment. This difference sequence is the high-frequency compensation component. This component contains the high-frequency fluctuations filtered out from the original command and needs to be handled by the energy storage system with a faster response speed. This high-frequency compensation component is used as the initial value for the energy storage system's execution command.

[0060] Then, the commands executed by the energy storage system are limited and corrected. The current state of charge of the energy storage system is collected in real time, that is, the percentage of the battery's remaining capacity relative to its rated capacity, with a value ranging from 0 to 100%. At the same time, the current available adjustable capacity of the wind power system is collected, that is, the power that the wind turbine can increase or decrease from its current output, in megawatts. Based on the sign and magnitude of the energy storage system's execution command, and considering the current state of charge (SBC) and available wind power regulation capacity, the following corrections are made: If the energy storage system's execution command is positive and the current SBC is below 20%, the energy storage discharge capacity is deemed insufficient, and the command value is reduced to 50% of its original value. This reduction is then added to the wind power system's execution command, but the accumulated wind power system execution command must not exceed its upper regulation capacity limit. If the energy storage system's execution command is negative and the current SBC is above 80%, the energy storage charging capacity is deemed insufficient, and the command value is reduced to 50% of its original value. This reduction is then added to the wind power system's execution command, but the accumulated wind power system execution command must not fall below its lower regulation capacity limit. If the absolute value of the energy storage system's execution command exceeds the rated power of the energy storage converter, it is limited to the rated power value. After these corrections, the final energy storage system execution command is obtained.

[0061] Finally, the revised energy storage system execution command and the aforementioned wind power system execution command are sent to the energy storage converter and the wind turbine power controller respectively through their respective remote communication devices. The energy storage system and the wind power system then coordinate to perform joint operation power regulation according to the command values.

[0062] The working principle of this invention is as follows: Real-time acquisition and preprocessing of automatic power generation control command signals from the power grid yields a standardized frequency modulation command sequence. This sequence is then decomposed into multiple intrinsic mode function (IMF) components with different center frequencies using variational mode decomposition. A Hilbert transform is performed on each IMF component to obtain instantaneous frequency and amplitude sequences. The instantaneous frequency sequence is matched with a pre-established frequency-lifetime damage mapping table based on electrochemical impedance spectroscopy (EIS) testing to obtain instantaneous damage coefficients. A weighted average of these coefficients, weighted by instantaneous amplitude, yields the comprehensive damage coefficient for each component, forming an IMF component set with damage coefficients. The rate of change of each component's damage coefficient over time is calculated to obtain the damage change rate. The comprehensive damage coefficient is then... The amplitude attenuation ratio is determined by inputting the damage rate and the number of damage components into a preset two-dimensional decision table. For components whose comprehensive damage coefficient exceeds the preset threshold, the instantaneous amplitude is compressed point by point based on the half-wave segment according to the amplitude attenuation ratio to obtain the adjusted intrinsic mode function component set. All adjusted components are accumulated at each time step to obtain the initial reconfiguration power command. The cutoff frequency of the finite impulse response low-pass filter is dynamically adjusted according to the current power change rate of the wind power system, and the low-frequency trend component is extracted as the wind power system execution command. The initial reconfiguration power command is subtracted from the low-frequency trend component to obtain the high-frequency compensation component as the energy storage system execution command. The energy storage system execution command is then amplitude-limited and corrected by combining the energy storage system's state of charge and the wind power system's available regulating capacity before being issued for execution.

[0063] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. An optimization method for wind-storage joint operation based on energy storage lifetime and frequency regulation performance, characterized in that, Includes the following steps: S1: Real-time acquisition of automatic power generation control command signals from the power grid, preprocessing of the command signals to obtain a standardized frequency regulation command sequence; S2: The standardized frequency modulation command sequence is decomposed into multiple eigenmode function components with different center frequencies through variational mode decomposition. S3: Match each intrinsic mode function component with the frequency band-lifetime damage mapping relationship established in advance based on electrochemical impedance spectroscopy test to obtain the lifetime damage coefficient corresponding to each component, and attach the lifetime damage coefficient to the corresponding component to form an intrinsic mode function component set with damage coefficient. S4: Filter the intrinsic mode function component set with damage coefficient, attenuate the amplitude of the component whose lifetime damage coefficient exceeds the preset threshold, and retain the component whose lifetime damage coefficient does not exceed the preset threshold to obtain the adjusted intrinsic mode function component set. S5: Linearly superimpose the adjusted intrinsic mode function component sets to reconstruct the optimized joint operating power command, and then send the joint operating power command to the energy storage system and wind power system for execution.

2. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 1, characterized in that, S2 specifically includes: A Fourier transform is performed on the standardized frequency modulation command sequence to obtain its spectral distribution. The number of mode decompositions for variational mode decomposition is determined based on the number of peaks in the spectrum whose energy amplitude exceeds a preset energy threshold. The initial center frequency of the variational mode decomposition is initialized with the center frequency corresponding to each peak. The variational optimization problem is solved iteratively using the alternating direction multiplier method. In each iteration, the center frequency and bandwidth of each intrinsic mode function component are updated, and the L2 norm of the difference between the components obtained in two adjacent iterations is calculated. The iteration is terminated when the L2 norm is less than the convergence threshold, and multiple intrinsic mode function components with different center frequencies are output.

3. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 2, characterized in that, The iterative solution of the variational optimization problem using the alternating direction multiplier method specifically includes: Initialize the Lagrange multipliers and the quadratic penalty factor; substitute the Lagrange multipliers and the quadratic penalty factor into the frequency domain update equation, and sequentially update each eigenmode function component in the frequency domain to obtain the updated eigenmode function components; update each center frequency according to the power spectrum centroid of the updated eigenmode function components to obtain the updated center frequency; update the Lagrange multipliers according to the updated eigenmode function components, the updated center frequency, and the Lagrange multipliers from the previous iteration to obtain the updated Lagrange multipliers; repeat the update process until the relative change of each eigenmode function component is less than the preset convergence threshold, and output the final eigenmode function components and their corresponding center frequencies.

4. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 1, characterized in that, S3 specifically includes: A Hilbert transform is performed on each intrinsic mode function component to obtain the instantaneous frequency sequence of each component. Each frequency value in the instantaneous frequency sequence is matched with a pre-constructed frequency band-lifetime damage mapping table to obtain the instantaneous damage coefficient corresponding to each frequency point. A weighted average is performed on all instantaneous damage coefficients within the same intrinsic mode function component, using the amplitude of the corresponding frequency point in the instantaneous frequency sequence as the weight, to calculate the comprehensive damage coefficient of the corresponding component. The comprehensive damage coefficient is assigned to the corresponding intrinsic mode function component as an additional attribute to form an intrinsic mode function component set with damage coefficients.

5. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 4, characterized in that, The calculation process for the instantaneous damage coefficient is as follows: Equivalent circuit fitting was performed on the impedance spectrum data obtained from electrochemical impedance spectroscopy to extract the characteristic curves of charge transfer resistance and diffusion resistance as a function of frequency. Based on the inflection point frequency of the sudden increase in internal resistance value in the characteristic curve, the frequency range was divided into a low-frequency diffusion region, a mid-frequency charge transfer region, and a high-frequency ohmic region, and a benchmark damage coefficient was set for each frequency range. Each frequency value in the instantaneous frequency sequence was sequentially assigned to its corresponding frequency range, and the benchmark damage coefficient corresponding to the frequency range was extracted. Based on the specific position of the frequency value within its range, linear interpolation was performed using the rate of change of internal resistance corresponding to the two endpoints of the frequency range to calculate the instantaneous damage coefficient corresponding to the frequency value.

6. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 1, characterized in that, S4 specifically includes: Calculate the rate of change of the damage coefficient over time in each intrinsic mode function component with a damage coefficient to obtain the damage change rate of each component. Input the damage coefficient and damage change rate of each component into a preset two-dimensional decision table. The two-dimensional decision table is divided into multiple decision regions with the damage coefficient as the horizontal axis and the damage change rate as the vertical axis. Each decision region corresponds to a different amplitude attenuation ratio. According to the landing region of each component in the two-dimensional decision table, read the corresponding amplitude attenuation ratio. Compress the instantaneous amplitude sequence of each component point by point according to the amplitude attenuation ratio to obtain the attenuated intrinsic mode function components. Components with damage coefficients not exceeding a preset threshold are retained as intrinsic mode function components. All attenuated components and retained components are combined to form the adjusted intrinsic mode function component set.

7. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 6, characterized in that, The step of compressing the instantaneous amplitude sequence of each component point by point according to the amplitude attenuation ratio to obtain the attenuated intrinsic mode function components specifically includes: Hilbert transform is performed on each attenuated intrinsic mode function component to obtain the instantaneous phase sequence of each component. The instantaneous amplitude sequence of each component is divided into continuous half-wave segments according to the zero-crossing point of the instantaneous phase sequence. The peak value of the instantaneous amplitude in each half-wave segment is calculated, and the peak value is multiplied by the corresponding amplitude attenuation ratio to obtain the target peak value of the corresponding half-wave segment. With the target peak value as the upper limit, all instantaneous amplitudes in the half-wave segment are proportionally limited and compressed so that the maximum amplitude of the compressed half-wave segment is equal to the target peak value. All compressed half-wave segments are spliced ​​together in their original time sequence to form the attenuated intrinsic mode function components.

8. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 1, characterized in that, S5 specifically includes: The initial reconfiguration power command is obtained by summing all components in the adjusted intrinsic mode function component set time-by-time. The initial reconfiguration power command is then input into a finite impulse response low-pass filter to extract the low-frequency trend component as the wind power system execution command. The low-frequency trend component is subtracted from the initial reconfiguration power command to obtain the high-frequency compensation component as the energy storage system execution command. The current state of charge of the energy storage system and the current available regulating capacity of the wind power system are collected. The energy storage system execution command is then limited and corrected based on the difference between the energy storage system execution command and the high-frequency compensation component. The corrected energy storage system execution command and the wind power system execution command are then sent to the corresponding systems for execution.

9. The wind-storage joint operation optimization method based on energy storage lifetime and frequency regulation performance according to claim 8, characterized in that, The step of accumulating all components in the adjusted intrinsic mode function component set time-by-time to obtain the initial reconfiguration power command specifically includes: The rate of change of the current output power of the wind power system is obtained, and the rate of change is input into a preset cutoff frequency mapping function. The mapping function outputs the current cutoff frequency according to the rule that the cutoff frequency increases with the rate of change. The filter coefficient sequence of the finite impulse response low-pass filter is calculated based on the current cutoff frequency. The sum of the coefficients in the filter coefficient sequence is 1 and they are symmetrically distributed. The initial reconstructed power command is convolved with the filter coefficient sequence to obtain a smoothed power sequence after moving average. The smoothed power sequence is output as a low-frequency trend component, which is used as the wind power system execution command.

Images