Multi-mode hybrid energy storage power distribution method in wind and light storage system

Abstract

The invention provides a multi-mode hybrid energy storage power distribution method in a wind and light storage system. The multi-mode hybrid energy storage power distribution method is a hybrid energy storage power distribution method based on improved crown porcupine algorithm-adaptive noise complete set empirical mode decomposition. The method comprises the following steps: firstly, carrying out Savitzky-Golay filtering on wind-solar composite power to obtain multi-mode HES reference power PHES, and laying a foundation for subsequent power distribution; secondly, an ICPO algorithm is established, and falling into a local minimum value is avoided while the convergence speed and the search quality of the algorithm are improved; an ICPO-ICEEMDAN method is provided, so that the decomposition effect on the PHES is improved; and finally, realizing intrinsic mode component reconstruction of the multi-mode HES according to the relative entropy, and performing power distribution on the multi-mode HES while reducing pseudo modes. According to the method, Cubic chaotic mapping and population reduction technologies are introduced to form an ICPO algorithm to optimize parameters of ICEEMDAN, IMF reconstruction of multi-modal HES is realized on the basis of relative entropy, examples show that the proposed power distribution method can reduce pseudo modals, power distribution conforming to energy storage characteristics is realized, the problem of power distribution when the number of HES modals is larger is solved, and the power distribution efficiency is improved. And compared with other heuristic algorithms, the advantages of the ICPO algorithm are verified.

Description

Technical Field

[0001] The present invention relates to the technical field of energy storage power distribution, and in particular to a multi-mode hybrid energy storage power distribution method in a wind-solar storage system. Background Art

[0002] With the continuous growth of the global economy and the continuous increase in population, energy demand is showing a rapid growth trend. Traditional fossil energy such as coal, oil and natural gas are not only limited in resources but also produce a large amount of greenhouse gases during use, which has a serious impact on the environment. Therefore, the global demand for renewable energy has risen sharply, and the proportion of new energy installed capacity, mainly wind and solar, has continued to increase. However, both wind and solar energy have problems of intermittency and randomness, resulting in large fluctuations in output power, which causes certain interference to the stable operation of the system. It is increasingly important to configure energy storage and form a wind, solar and storage joint operation system. Energy storage has multiple modes, such as power-type flywheels and energy-type lithium batteries. Hybrid Energy Storage (HES) can give full play to the complementary characteristics of different energy storage technologies, balance the volatility of wind and solar power generation, improve the overall efficiency and stability of the system, and achieve more efficient and flexible energy storage and management.

[0003] Since hybrid energy storage systems have the advantages of smoothing power fluctuations, enhancing grid resilience, and improving grid recovery capabilities, their position in the power system is becoming increasingly important, and their development trend is also diversified and rapid. Reasonable power allocation of hybrid energy storage systems can optimize the charging and discharging process of different energy storage units, reduce energy loss and waste, and improve the utilization efficiency of energy storage resources. At the same time, it helps to reduce the construction and operation costs of hybrid energy storage systems, and improve the redundancy of power systems to ensure continuous power supply to power systems. At present, most power allocation strategies are based on a given "1 energy type + 1 power type" two-mode HES solution to study power allocation. In actual applications, HES is often affected by factors such as the type of energy storage technology, the complexity of system design, and the specific needs of the application scenario. The number of modes is not limited to two modes and is not fixed. Therefore, studying the power allocation method of multi-modal HES has become a new hotspot. Summary of the invention

[0004] In view of the deficiencies in the prior art, the present invention aims to provide a multi-modal hybrid energy storage power allocation method in a wind-solar-storage system to solve the power allocation problem when the number of HES modes is greater. It is a hybrid energy storage power allocation method based on the Improved Crested Porcupine Optimizer (ICPO)-Adaptive Noise Complete Ensemble Empirical Mode Decomposition (ICEEMDAN). First, the wind-solar composite power is filtered by Savitzky-Golay (SG) to obtain the multi-modal HES reference power P HES , laying the foundation for subsequent power allocation; then the ICPO algorithm is established to improve the convergence speed and search quality of the algorithm while avoiding falling into the local minimum; and the ICPO-ICEEMDAN method is proposed to improve the P HES decomposition effect; finally, the intrinsic mode component (IMF) of the multi-modal HES is reconstructed based on the relative entropy, so as to achieve the purpose of power allocation for the multi-modal HES while reducing the pseudo-mode.

[0005] In order to solve the above technical problems, the technical solution provided by the present invention is:

[0006] A multi-mode hybrid energy storage power distribution method in a wind-solar storage system comprises the following steps:

[0007] S1. Perform Savitzky-Golay filtering on the wind-solar composite power to obtain the multi-modal HES reference power P HES , laying the foundation for subsequent power allocation;

[0008] S2. Set the maximum number of iterations t, population size N, number of noise additions Nstd and the optimization range of amplitude weight Ne, introduce Cubic chaotic mapping as a random source, introduce the number of cycles to reduce the population, establish the ICPO algorithm, and perform Cubic chaos initialization of the ICPO population;

[0009] S3. Under the ICPO algorithm established in S2, the noise addition process and iteration mechanism are improved, and the ICPO-ICEEMDAN method is obtained by adaptively adjusting the noise level added to the signal. ICPO-ICEEMDAN is used to decompose the P obtained in S1. HES, the fitness function is iterated with the minimum composite entropy of the IMF component as the fitness function. When the number of iterations is less than t, the population position is updated and the iteration continues. When the number of iterations is greater than or equal to t, the optimal parameter combination is obtained, and the optimal IMF component is output based on the parameter combination.

[0010] S4. Based on the optimal IMF components obtained in S3, the high and low frequency areas are divided according to the relative entropy, and the areas are divided according to the rated number of cycles. P power-type and q energy-type energy storages are responsible for the IMF components corresponding to the high and low frequency areas respectively, realizing the power distribution of the n-mode HES.

[0011] Preferably, at time t, wind power, photovoltaic power, composite power, HES output power and direct grid-connected power are represented by P wind (t), P pv (t), P fuhe (t), P HES (t), P grid (t); said S1 comprises the following steps,

[0012] S1.1 Establish the relationship between composite output and output power

[0013] P fuhe The calculation formula of (t) is as follows:

[0014] P fuhe (t) = P wind (t)+P pv (t) (1)

[0015] S1.2 Calculate the direct grid-connected power P grid (t)

[0016] SG filtering uses a sliding window method to perform high-order polynomial fitting in each window to smooth the signal and directly connect the grid power P grid The calculation formula of (t) is as follows, and its fluctuation degree meets the preset grid connection standard;

[0017]

[0018] In the formula, c j is the weighting coefficient, which is related to the number of polynomial least squares fitting; w is the window size; P ref is the limit value of fluctuation within Δt time;

[0019] S1.3 Calculation of HES output power

[0020] P HES The expression of (t) is as follows:

[0021] P HES (t) = P fuhe(t)-P grid (t) (3).

[0022] Preferably, said S2 comprises the following steps:

[0023] S2.1 Classification of CP defense strategies

[0024] As the predator gets closer to the center of the crested porcupine, the CP's defense strategy gradually upgrades, from the least aggressive to the most aggressive, which are divided into Zone A - flapping back spines, Zone B - making noise, Zone C - secreting odors, and Zone D - backstab attacks;

[0025] S2.2 Improve CPO algorithm and establish ICPO algorithm

[0026] Cubic chaotic mapping is introduced as a random source, and cyclic population reduction technology is introduced to establish ICPO to improve the convergence speed and search quality of the algorithm while avoiding falling into local minima. Specifically:

[0027] S2.2.1 Chaotic map initialization

[0028] CPO starts the search process from the initial set of individuals, and the calculation formula is as follows:

[0029] x i =L+r×(UL),i=1,2...,N (4)

[0030] x i+1 =ρx i ×(1-x i 2 ) (5)

[0031] Where N represents the population size; x i is the i-th candidate solution in the search space; L and U are the lower and upper limits of the search range respectively; r is a random value in [0,1]; ρ is the Cubic mapping control parameter;

[0032] The cyclic population reduction technology means that not all CPs activate the defense mechanism. Only the threatened CPs will defend. The calculation formula is as follows:

[0033]

[0034] In the formula, T is an integer variable that determines the number of cycles, T max is the maximum value; t is the number of iterations; % represents the remainder operation; N min is the minimum number of individuals in the population. As t increases, the size of the new population N' gradually decreases;

[0035] S2.2.2 Exploration stage

[0036] Area A: After the CP finds the predator, it flaps its spines to deter it. If the predator still moves towards the CP, the distance between them decreases, and the predator is encouraged to explore area A. If the predator chooses to leave, the distance between them is maximized to explore whether there is a solution in the unknown area. The calculation formula is as follows:

[0037]

[0038] In the formula, x t CP is the solution at iteration t, y t i is the vector between the current CP and a randomly selected CP from the population, used to represent the position of the predator at time t, τ1 is a normally distributed random number, τ2 is a random value in [0,1], and r0 is a random value in [1,N];

[0039] Area B: CP makes noise to deter predators. When a predator approaches it, the sound becomes louder, otherwise it becomes smaller. The calculation formula is as follows:

[0040]

[0041] Where τ3 is a random value in [0,1], r1 and r2 are random numbers in [1,N];

[0042] When |U1|=0, the predator stops moving, the distance remains unchanged, and the sound is weakened; when |U1|=1, the predator is still nearby and the sound is strengthened; otherwise, the sound is normal;

[0043] S2.2.3 Development stage

[0044] Area C: CP secretes odor to prevent predators from approaching. The calculation formula is as follows:

[0045]

[0046] Where r3 is a random number in [1,N], δ is the search direction control parameter, and γ t is the defense factor, S t i is the odor diffusion factor;

[0047] When |U1|=0, the predator stops moving, the distance remains unchanged, and the predator stops emitting odor; when |U1|=1, the predator is still nearby and emits a stronger odor; otherwise, the predator maintains a safe distance from the CP and emits a general odor;

[0048] Area D: CP uses backstabs to attack predators, analogous to a one-dimensional inelastic collision, and the calculation formula is as follows:

[0049]

[0050] In the formula, x t CP is the optimal solution, namely CP; x t i represents the position of the i-th individual at iteration t, i.e., the predator; β is the convergence speed factor, τ4 is a random value in [0,1]; F i t To affect the attack power of the CP of the i-th predator, the relevant calculation formula is as follows:

[0051]

[0052] In the formula, is the mass of the ith individual at iteration t, f is the objective function, is the final velocity of the ith individual at time t+1, and is allocated based on the random solution selected from the current population, is the initial velocity of the ith individual at iteration t, τ6 is a random vector in [0,1], ε is a minimum value to avoid being divided by 0, and r4 is a random number in [1,N].

[0053] Preferably, S3 is specifically:

[0054] S3.1 Decomposition of P using ICEEMDAN HES

[0055] S3.1.1 Use empirical mode decomposition to decompose Gaussian white noise and obtain the IMF component of the noise l i(i＝1,2,.,N) , add the original signal x(t) to each component after processing, and get the extended signal x i (t):

[0056] x i (t) = x(t) + a0l i (15)

[0057] Where, α0=βstd(x(t)) / std(l1) is the standard deviation of white noise, β is a constant, and std represents the standard deviation operation;

[0058] S3.1.2 Define the first residual component res1:

[0059]

[0060] In the formula, M[x i (t)] represents x i (t) minus the first component IMF1 obtained by EMD decomposition. At this step, the first component IMF1 of the original signal x(t) is as follows:

[0061] IMF1=x(t)-res1 (17)

[0062] S3.1.3 Calculate the second residual component res2 and derive IMF2:

[0063]

[0064] IMF2=res1-res2 (19)

[0065] Where, α1 = βstd(r1);

[0066] S3.1.4 Use the same method to calculate the kth residual and IMF:

[0067]

[0068] IMF k =res k-1 -res k (twenty one)

[0069] In the formula, α k-1 =βstd(r k-1 ),k>1;

[0070] S3.1.5 Iterate until the maximum number of iterations is met;

[0071] S3.2 constructs composite entropy as an evaluation index. The composite entropy is composed of the envelope entropy E p and mutual information entropy MI;

[0072] S3.2.1 Calculation of envelope entropy E p

[0073] Envelope entropy can reflect whether the frequency characteristics of each IMF are obvious. When the IMF characteristic information is less, the envelope entropy value is larger, otherwise, it is smaller. p It can be expressed as:

[0074]

[0075] Where a(j) is the envelope signal obtained after Hilbert modulation of signal x(j); p j is the normalized form of a(j);

[0076] S3.2.2 Calculation of mutual information entropy MI

[0077] Mutual information in information theory indicates the degree of correlation between two events and is not easily affected by external factors. The larger the mutual information entropy value MI, the stronger the correlation between the two events. As far as IMF is concerned, the richer the original signal feature information contained is. The calculation formula is as follows:

[0078] MI(X,Y)=H(Y)-H(Y|X) (24)

[0079] In the formula, X and Y are different events, H(Y) is the entropy of Y, and X is P HES , Y is IMF, H(Y|X) is the conditional entropy of Y when X is known;

[0080] S3.2.3 Establishing a composite fitness function

[0081] The expression of the composite fitness function is as follows:

[0082]

[0083] ICEEMDAN improves the stability and robustness of the results by repeatedly adding noise and performing decomposition;

[0084] S3.2.4 Select Nstd and Ne as optimization objects, use compound entropy as fitness function, and use ICPO to optimize the parameter combination [Nstd, Ne] to improve ICEEMDAN's P HES Decomposition effect.

[0085] Preferably, S4 is specifically:

[0086] P HES After ICPO-ICEEMDAN decomposition, the spectrum is obtained by fast Fourier transform;

[0087] S4.1 Calculation of relative entropy D

[0088] Considering the relative entropy between the IMF spectra, the relative entropy between the IMF spectra is calculated to characterize the degree of spectrum differentiation. The IMF is divided into n combinations. The relative entropy D is calculated as follows:

[0089]

[0090] In the formula, x is the amplitude and N is the number of time points;

[0091] S4.2 Setting power distribution point

[0092] The maximum value of relative entropy D is taken as the high frequency component P H and low frequency component P L The dividing point K is obtained from the dividing point K, and after the high-frequency area and the low-frequency area are obtained, the p power types and q energy types in the HES selection results are further allocated;

[0093] For the high-frequency area, the relative entropies are arranged in descending order, and the first p-1 points with larger relative entropy values ​​are taken as power points, and the high-frequency area is divided into p segments. The same is true for the low-frequency area. The energy storage with a larger rated cycle number is selected as the bearer of the power in the higher frequency band. The power point setting of M IMFs can be summarized as follows: the maximum relative entropy value is used as the dividing point, the left side contains K IMF powers, which is the high-frequency area and is borne by the power-type energy storage; the right side contains MK IMF powers, which is the low-frequency area and is borne by the energy-type energy storage.

[0094] The decomposed residual res is placed in the low-frequency area and borne by the energy type with a larger number of cycles;

[0095] S4.3 Reconstruct the power of each mode of energy storage

[0096] The calculation formula of each mode energy storage power is as follows:

[0097]

[0098] Where x, y, x', y' are the relative entropy division points.

[0099] The beneficial effects of the present invention are:

[0100] The present invention studies the multimodal power distribution of HES in wind-solar-storage system, establishes a hybrid energy storage power distribution strategy based on ICPO-ICEEMDAN method, and optimizes the parameters of ICEEMDAN by introducing Cubic chaotic mapping and population reduction technology to form ICPO algorithm, and realizes IMF reconstruction of multimodal HES based on relative entropy. The example shows that the proposed power distribution method can reduce pseudo-modes and realize power distribution that meets the characteristics of energy storage. The advantages of ICPO algorithm are verified by comparing with other heuristic algorithms. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0102] Figure 1 Flow chart of the multi-mode hybrid energy storage power allocation method in the wind-solar storage system of the present invention

[0103] Figure 2 CP defense mechanism diagram

[0104] Figure 3 Set up diagram for power split point

[0105] Figure 4 The wind and solar power output diagram for the embodiment

[0106] Figure 5HES reference power diagram for example

[0107] Figure 6 This is the power decomposition result diagram of the example ICPO-ICEEMDAN

[0108] Figure 7 This is the power output result diagram of the example "lead-carbon battery + flow battery, supercapacitor"

[0109] Figure 8 Convergence curve comparison diagram for the embodiment

[0110] Figure 9a The power decomposition result diagram of ICEEMDAN is used in the embodiment. Figure 9b The power decomposition result diagram using VMD

[0111] Fig.10 The VMD power output result diagram is shown in the example DETAILED DESCRIPTION

[0112] The preferred embodiments of the present invention are described below in conjunction with the accompanying drawings. It should be understood that the following embodiments are provided only for the purpose of illustration and are not intended to limit the scope of the present invention. Those skilled in the art may make various modifications and substitutions to the present invention without departing from the purpose and spirit of the present invention.

[0113] The present invention provides a multi-modal hybrid energy storage power allocation method in a wind-solar-storage system to solve the power allocation problem when the number of HES modes is larger. It is a hybrid energy storage power allocation method based on the Improved Crested Porcupine Optimizer (ICPO)-Adaptive Noise Complete Ensemble Empirical Mode Decomposition (ICEEMDAN). First, the wind-solar composite power is filtered by Savitzky-Golay (SG) to obtain the multi-modal HES reference power P HES , laying the foundation for subsequent power allocation; then the ICPO algorithm is established to improve the convergence speed and search quality of the algorithm while avoiding falling into the local minimum; and the ICPO-ICEEMDAN method is proposed to improve the P HES The decomposition effect is obtained; finally, the intrinsic mode function (IMF) of the multi-modal HES is reconstructed based on relative entropy, and the purpose of power allocation for the multi-modal HES is achieved while reducing pseudo-modes. The specific steps are as follows:

[0114] 1. Calculation of multimodal HES reference power

[0115] (1) At time t, wind power, photovoltaic power, composite power, HES output power and grid connection point power are expressed as P wind (t), P pv (t), P fuhe (t), P HES (t), P grid (t), P fuhe The calculation formula of (t) is as follows:

[0116] P fuhe (t) = P wind (t)+P pv (t) (1)

[0117] (2) Calculation of direct grid-connected power P grid (t)

[0118] SG filtering uses a sliding window method to perform high-order polynomial fitting in each window to achieve signal smoothing. The calculation formula for direct grid-connected power is as follows: It is the direct grid-connected power obtained by SG filtering, and its fluctuation degree meets the pre-set grid-connected standard.

[0119]

[0120] In the formula, c j is the weighting coefficient, which is related to the number of polynomial least squares fitting; w is the window size; P ref It is the limit value of fluctuation within Δt time.

[0121] (3) Calculate HES output power, P HES The expression of (t) is as follows:

[0122] P HES (t) = P fuhe (t)-P grid (t) (3)

[0123] 2. Establish the Improved Crested Porcupine Optimization Algorithm (ICPO)

[0124] (1) Classification of CP defense strategies

[0125] Crested Porcupine Optimizer (CPO) is a new meta-heuristic algorithm. As the predator gets closer to the center of the Crested Porcupine (CP), the CP's defense strategy gradually upgrades. From the least aggressive to the most aggressive, they are Area A - flapping backstabs, Area B - making noise, Area C - secreting odors, and Area D - backstab attacks. Figure 2As shown in Figure 2. CPO uses a random generation method to initialize the population, which will lead to uneven distribution and easy to fall into local extreme values. However, chaotic mapping has the characteristics of randomness and pervasiveness.

[0126] (2) Improve the CPO algorithm and establish the ICPO algorithm

[0127] 1) Methods for improving CPO

[0128] The present invention introduces Cubic chaotic mapping as a random source and cyclic population reduction technology to establish ICPO, so as to improve the convergence speed and search quality of the algorithm while avoiding falling into a local minimum.

[0129] 2) ICPO algorithm process

[0130] ① Chaotic map initialization

[0131] CPO starts the search process from the initial set of individuals, and the calculation formula is as follows:

[0132] x i =L+r×(UL),i=1,2...,N (4)

[0133] x i+1 =ρx i ×(1-x i 2 ) (5)

[0134] Where N represents the population size; x i is the i-th candidate solution in the search space; L and U are the lower and upper limits of the search range respectively; r is a random value in [0,1]; ρ is the Cubic mapping control parameter.

[0135] The cyclic population reduction technology means that not all CPs activate the defense mechanism. Only the threatened CPs will defend. The calculation formula is as follows:

[0136]

[0137] In the formula, T is an integer variable that determines the number of cycles, T max is the maximum value; t is the number of iterations; % represents the remainder operation; N min is the minimum number of individuals in the population. As t increases, the size of the new population N' gradually decreases.

[0138] ② Exploration stage

[0139] Area A: After the CP finds the predator, it flaps its spines to deter it. If the predator still moves towards the CP, the distance between them decreases. This choice encourages exploration of area A. If the predator chooses to leave, the distance between them is maximized to explore whether there is a solution in the unknown area. The calculation formula is as follows:

[0140]

[0141] In the formula, x t CP is the solution at iteration t, y t i is the vector between the current CP and a randomly selected CP from the population, used to represent the position of the predator at time t, τ1 is a normally distributed random number, τ2 is a random value in [0,1], and r0 is a random value in [1,N].

[0142] Area B: CP makes noise to deter predators. When a predator approaches it, the sound becomes louder, otherwise it becomes smaller. The calculation formula is as follows:

[0143]

[0144] Where τ3 is a random value in [0,1], and r1 and r2 are random numbers in [1,N]. When |U1|=0, the predator stops moving, the distance remains unchanged, and the sound is weakened; when |U1|=1, the predator is still nearby and the sound is strengthened; otherwise, it emits a normal sound.

[0145] ③ Development stage

[0146] Area C: CP secretes odor to prevent predators from approaching. The calculation formula is as follows:

[0147]

[0148] Where r3 is a random number in [1,N], δ is the search direction control parameter, and γ t is the defense factor, S i t is the odor diffusion factor. When |U1|=0, the predator stops moving, the distance remains unchanged, and the predator stops emitting odor; when |U1|=1, the predator is still nearby and emits a stronger odor; otherwise, the predator maintains a safe distance from the CP and emits a general odor.

[0149] Area D: CP uses backstabs to attack predators, analogous to a one-dimensional inelastic collision, and the calculation formula is as follows:

[0150]

[0151] In the formula, x t CP is the optimal solution, namely CP; x t i represents the position of the i-th individual at iteration t, i.e., the predator; β is the convergence speed factor, τ4 is a random value in [0,1]; F i tTo affect the attack power of the CP of the i-th predator, the relevant calculation formula is as follows:

[0152]

[0153] In the formula, is the mass of the ith individual at iteration t, f is the objective function, is the final velocity of the ith individual at time t+1, and is allocated based on the random solution selected from the current population, is the initial velocity of the ith individual at iteration t, τ6 is a random vector in [0,1], ε is a minimum value to avoid being divided by 0, and r4 is a random number in [1,N].

[0154] 3. Propose ICPO-ICEEMDAN method

[0155] (1) Principle of ICEEMDAN

[0156] ICEEMDAN ensures that all information of the original signal can be retained after decomposition, improves the correctness of decomposition, and has been widely used in signal processing. Different from the traditional method of adding ordinary Gaussian white noise, ICEEMDAN improves the noise addition process and iterative mechanism, and improves the decomposition quality and reduces the modal aliasing phenomenon by adaptively adjusting the noise level added to the signal.

[0157] Decomposition of P using ICEEMDAN HES The specific decomposition process is as follows:

[0158] 1) Use Empirical Mode Decomposition (EMD) to decompose Gaussian white noise and obtain the IMF component of the noise l i(i＝1,2,.,N) , add the original signal x(t) to each component after processing, and get the extended signal x i (t):

[0159] x i (t) = x(t) + a0l i (15)

[0160] Where: α0=βstd(x(t)) / std(l1) is the standard deviation of white noise, β is a constant, and std represents the standard deviation operation.

[0161] 2) Define the first residual component res1:

[0162]

[0163] Where: M[x i (t)] represents x i(t) minus the first component IMF1 obtained by EMD decomposition. At this step, the first component IMF1 of the original signal x(t) is as follows:

[0164] IMF1=x(t)-res1 (17)

[0165] 3) Calculate the second residual component res2 and derive IMF2.

[0166]

[0167] IMF2=res1-res2 (19)

[0168] Where: α1 = βstd(r1).

[0169] 4) Similarly, use the same method to calculate the kth residual and IMF:

[0170]

[0171] IM F k =res k-1 -res k (twenty one)

[0172] Where: α k-1 =βstd(r k-1 ),k>1.

[0173] 5) Iterate until the maximum number of iterations is met.

[0174] (2) Constructing composite entropy as an evaluation index

[0175] The composite entropy is given by the envelope entropy E p and mutual information entropy MI.

[0176] 1) Envelope entropy E p

[0177] Envelope entropy can reflect whether the frequency characteristics of each IMF are obvious. When the IMF feature information is less, the envelope entropy value is larger, otherwise, it is smaller. p It can be expressed as:

[0178]

[0179] Where a(j) is the envelope signal obtained after Hilbert modulation of signal x(j); p j is the normalized form of a(j).

[0180] 2) Mutual Information Entropy MI

[0181] Mutual information is used in information theory to indicate the degree of correlation between two events and is not easily affected by external factors. The larger the mutual information entropy value MI, the stronger the correlation between the two events. In terms of IMF, the richer the original signal feature information contained is. The calculation formula is as follows:

[0182] MI(X,Y)=H(Y)-H(Y|X) (24)

[0183] In the formula, X and Y are different events, H(Y) is the entropy of Y, and X is P HES , Y is the IMF, and H(Y|X) is the conditional entropy of Y when X is known.

[0184] (3) Establishing a composite fitness function

[0185] The expression of the composite fitness function is as follows:

[0186]

[0187] ICEEMDAN improves the stability and robustness of the results by repeatedly adding noise and performing decomposition. Among them, the two parameters of the number of noise additions Nstd and the amplitude weight Ne have a significant impact on the decomposition effect. A larger Nstd can enhance the processing ability of complex and IMF modes in the signal, but it will also increase the computational cost and affect the accuracy and stability of the decomposition; a larger Ne will reduce the accuracy of signal decomposition, and a smaller Ne will cause aliasing in the decomposition results.

[0188] The present invention selects Nstd and Ne as optimization objects, takes composite entropy as fitness function, and uses ICPO to optimize the parameter combination [Nstd, Ne] to realize ICEEMDAN for P HES Improved decomposition effect.

[0189] 4. Perform multi-modal HES power allocation

[0190] P HES After ICPO-ICEEMDAN decomposition, the spectrum can be obtained through Fast Fourier Transform (FFT).

[0191] (1) Relative entropy D

[0192] Considering the relative entropy between the IMF spectra, the relative entropy between the IMF spectra is calculated to characterize the degree of spectrum differentiation, and the IMF is divided into n combinations. The relative entropy D calculation formula is as follows:

[0193]

[0194] Where x is the amplitude and N is the number of time points.

[0195] (2) Setting power distribution points

[0196] In the calculation results, the relative entropy D(IMF i , IMF i+1 ) must have a process of changing from large to small, so its maximum value point is taken as the high-frequency component P H and low frequency component P L The dividing point K. After obtaining the high-frequency zone and the low-frequency zone from the dividing point K, the p power types and q energy types in the HES selection results are further allocated. For the high-frequency zone, the relative entropies are arranged in order from large to small, and the first p-1 points with larger relative entropy values ​​are taken as power distribution points, and the high-frequency zone is divided into p segments, and the low-frequency zone is the same. Taking into account the life loss of different types of energy storage, the energy storage with a larger rated cycle number is selected as the bearer of the power in the higher frequency band, so as to achieve the effect of overcoming the shortcomings of the small capacity of the power type and reducing the number of cycles of the energy type and extending its service life. The power distribution point settings of M IMFs are as follows: Figure 3 As shown, x, y, x', y' are relative entropy partition points.

[0197] Considering that the residual res mainly contains lower frequency power after decomposition, it is placed in the low frequency area and borne by the energy type with a larger number of cycles.

[0198] (3) Reconstruct the power of each mode of energy storage

[0199] The calculation formula of each mode energy storage power is as follows:

[0200]

[0201] Example

[0202] In this embodiment, the HES mode number n is 2 or 3, and p, q ≥ 1. There are 6 commercial energy storage options available for the HES solution, including 4 energy storage types: lead-carbon battery, lithium-ion battery, sodium-sulfur battery, and flow battery, and 2 power storage types: supercapacitor and flywheel. Finally, the HES solution consisting of "lead-carbon battery + flow battery, supercapacitor" is selected for power allocation.

[0203] 1. Basic data settings

[0204] The actual sampling data of a certain wind and solar power station in Northwest my country on a certain day is taken as the research object, with a wind power installed capacity of 100MW and a photovoltaic installed capacity of 200MW. The sampling interval is 1min, P fuhe like Figure 4 shown.

[0205] Set the polynomial order in the SG filtering algorithm to 6, the w iteration size to 20, the ICPO population size to 4, the number of iterations to 100, and the optimization ranges of Nstd and Ne to [50,600] and [0.15,0.6] respectively. grid The grid connection standards are maximum fluctuations of 15MW and 50MW for 1min and 10min respectively.

[0206] 2. Analysis of power allocation results

[0207] Figure 4 P fuhe After the fluctuation is smoothed by SG filtering, the power P that meets the grid connection standard is obtained. grid , and further obtain the reference power P HES ,like Figure 5 As shown. Among them, P grid The maximum fluctuations of 1min and 10min are 11.69MW and 22.06MW respectively.

[0208] The best fitness value calculated by ICPO is 6.3138, and the best position is (0.2915, 114.563). This optimal combination is used as the ICEEMDAN parameter for P HES Decomposition, the result is Figure 6 shown.

[0209] The relative entropy between the IMFs of the HES scheme - "lead-carbon battery + flow battery, supercapacitor" is calculated, as shown in Table 1.

[0210] Table 1 Relative entropy calculation results

[0211]

[0212] Judging from the results obtained, the relative entropy between IMF1 and IMF2 is the largest, so IMF1 is divided into the high-frequency area, and IMF2-10 is divided into the low-frequency area. The selected HES mode is 3, the energy types are lead-carbon battery and liquid flow battery, and the power type is supercapacitor. The IMF1 component in the high-frequency area is borne by the supercapacitor. In the low-frequency area, the relative entropy between IMF5 and IMF6 is the largest, reaching 126.4. Therefore, considering the life and number of cycles of lead-carbon batteries and liquid flow batteries, liquid flow batteries are selected to bear IMF2-5, and lead-carbon batteries are selected to bear IMF6-10 and the residual component res. In summary, the power distribution results of each mode of energy storage of this HES scheme are shown in Table 2, and the power output is shown in Figure 7 shown.

[0213] Table 2 Power allocation results of “lead-carbon battery + flow battery, supercapacitor”

[0214]

[0215] from Figure 7 It can be seen that supercapacitors have a short working cycle and are charged and discharged very frequently; lead-carbon batteries and flow batteries have long working cycles and are charged and discharged relatively slowly. The former's charging and discharging power is basically between -4 and 4MW, handling high-frequency power in the low-frequency zone, while the latter is basically between -2 and 2MW, handling lower-frequency power in the low-frequency zone. This reflects the characteristics of supercapacitors as power-type energy storage devices that can release a large amount of electricity in a short time, and lead-carbon batteries and flow batteries as energy-type energy storage devices that can be charged and discharged for a long time. Since supercapacitors have the characteristics of reversible charging and discharging, the number of cycles can reach hundreds of thousands, so they can withstand frequent power fluctuations. Lead-carbon batteries and flow batteries have fewer charging and discharging changes, thereby reducing the number of battery cycles and slowing down their attenuation and aging.

[0216] 3. ICPO convergence analysis

[0217] In order to verify the improvement of ICPO convergence after the introduction of Cubic chaotic mapping, the ICPO algorithm is compared with CPO, Particle Swarm Optimizer (PSO), Sparrow Search Algorithm (SSA) and other algorithms. Figure 8 shown. Figure 8 ICPO converges at the 8th iteration, which is faster than CPO and other algorithms, and its fitness function value is also the smallest among the four algorithms, further indicating that ICPO has higher convergence accuracy.

[0218] 4. Analysis of the advantages of ICPO-ICEEMDAN method in optimizing parameters

[0219] In order to verify the advantages of ICPO-ICEEMDAN with optimized parameters, the traditional ICEEMDAN and variational mode decomposition (VMD) are used to analyze the P HES The decomposition result is shown in Figure 9.

[0220] Depend on Figure 9a It can be seen that the number of pseudo modes in the ICEEMDAN decomposition results increases, which increases the number of signal processing, affects the reliability of the analysis results, and is not conducive to the rapid power allocation of HES. Figure 9b After the VMD decomposition result is transformed by FFT (sampling rate 1000), the center frequency of IMF1 is the smallest, but as high as 125Hz. HES When the mutation level is high, VMD has difficulty in converting P HESThe low-frequency components suitable for energy type are decomposed. Another HES solution "lithium battery, supercapacitor" is selected as the VMD power allocation object. The power reconstruction result is as follows Fig.10 As shown in the figure, the number of energy-type charge and discharge times reaches hundreds, which seriously affects the energy storage life. This result also indirectly verifies that the use of three-mode HES with higher scores can better utilize the complementary advantages of energy storage in technology, economy and environmental safety, and reduce energy-type life loss.

[0221] The contents not described in detail in this specification belong to the prior art known to professional and technical personnel in this field.

[0222] Finally, it should be noted that the above description is only a preferred example of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A multi-mode hybrid energy storage power allocation method in a wind-solar storage system, characterized in that: The following steps are included: S1. Perform Savitzky-Golay filtering on the wind-solar composite power to obtain the multi-modal HES reference power P HES , laying the foundation for subsequent power allocation; S2. Set the maximum number of iterations t, population size N, number of noise additions Nstd and the optimization range of amplitude weight Ne, introduce Cubic chaotic mapping as a random source, introduce the number of cycles to reduce the population, establish the ICPO algorithm, and perform Cubic chaos initialization of the ICPO population; S3. Under the ICPO algorithm established in S2, the noise addition process and iteration mechanism are improved, and the ICPO-ICEEMDAN method is obtained by adaptively adjusting the noise level added to the signal. ICPO-ICEEMDAN is used to decompose the P obtained in S1. HES , the fitness function is iterated with the minimum composite entropy of the IMF component as the fitness function. When the number of iterations is less than t, the population position is updated and the iteration continues. When the number of iterations is greater than or equal to t, the optimal parameter combination is obtained, and the optimal IMF component is output based on the parameter combination. S4. Based on the optimal IMF components obtained in S3, the high and low frequency areas are divided according to the relative entropy, and the areas are divided according to the rated number of cycles. P power-type and q energy-type energy storages are responsible for the IMF components corresponding to the high and low frequency areas respectively, realizing the power distribution of the n-mode HES.

2. A multi-mode hybrid energy storage power allocation method in a wind-solar storage system according to claim 1, characterized in that: At time t, wind power, photovoltaic power, composite power, HES output power and direct grid-connected power are represented by P wind (t), P pv (t), P fuhe (t), P HES (t), P grid (t); said S1 comprises the following steps, S1.1 Establish the relationship between composite output and output power P fuhe The calculation formula of (t) is as follows: P fuhe (t)=P wind (t)+P pv (t) (1) S1.2 Calculate the direct grid-connected power P grid (t) SG filtering uses a sliding window method to perform high-order polynomial fitting in each window to smooth the signal and directly connect the grid power P grid The calculation formula of (t) is as follows, and its fluctuation degree meets the preset grid connection standard; In the formula, c j is the weighting coefficient, which is related to the number of polynomial least squares fitting; w is the window size; P ref is the limit value of fluctuation within Δt time; S1.3 Calculation of HES output power P HES The expression of (t) is as follows: P HES (t)=P fuhe (t)-P grid (t) (3)。 3. The multi-mode hybrid energy storage power allocation method in a wind-solar storage system according to claim 1 is characterized in that: S2 The following steps are included: S2.1 Classification of CP defense strategies As the predator gets closer to the center of the crested porcupine, the CP's defense strategy gradually upgrades, from the least aggressive to the most aggressive, which are divided into Zone A - flapping back spines, Zone B - making noise, Zone C - secreting odors, and Zone D - backstab attacks; S2.2 Improve CPO algorithm and establish ICPO algorithm Cubic chaotic mapping is introduced as a random source, and cyclic population reduction technology is introduced to establish ICPO to improve the convergence speed and search quality of the algorithm while avoiding falling into local minima. Specifically: S2.2.1 Chaotic map initialization CPO starts the search process from the initial set of individuals, and the calculation formula is as follows: x i =L+r×(U-L),i=1,2...,N (4) x i+1 =ρx i ×(1-x i 2 ) (5) Where N represents the population size; x i is the i-th candidate solution in the search space; L and U are the lower and upper limits of the search range respectively; r is a random value in [0,1]; ρ is the Cubic mapping control parameter; The cyclic population reduction technology means that not all CPs activate the defense mechanism. Only the threatened CPs will defend. The calculation formula is as follows: Where T is an integer variable that determines the number of cycles, T max is the maximum value; t is the number of iterations; % represents the remainder operation; N min is the minimum number of individuals in the population. As t increases, the size of the new population N' gradually decreases; S2.2.2 Exploration stage Area A: After the CP finds the predator, it flaps its spines to deter it. If the predator still moves towards the CP, the distance between them decreases, and the predator is encouraged to explore area A. If the predator chooses to leave, the distance between them is maximized to explore whether there is a solution in the unknown area. The calculation formula is as follows: In the formula, x t CP is the solution at iteration t, y t i is the vector between the current CP and a randomly selected CP from the population, used to represent the position of the predator at time t, τ1 is a normally distributed random number, τ2 is a random value in [0,1], and r0 is a random value in [1,N]; Area B: CP makes noise to deter predators. When a predator approaches it, the sound becomes louder, otherwise it becomes smaller. The calculation formula is as follows: Where τ3 is a random value in [0,1], r1 and r2 are random numbers in [1,N]; When |U1|=0, the predator stops moving, the distance remains unchanged, and the sound is weakened; when |U1|=1, the predator is still nearby, and the sound is strengthened; Otherwise, emit a general sound; S2.2.3 Development stage Area C: CP secretes odor to prevent predators from approaching. The calculation formula is as follows: Where r3 is a random number in [1,N], δ is the search direction control parameter, and γ t is the defense factor, S t i is the odor diffusion factor; When |U1|=0, the predator stops moving, the distance remains unchanged, and the predator stops emitting odor; when |U1|=1, the predator is still nearby and emits a stronger odor; otherwise, the predator maintains a safe distance from the CP and emits a general odor; Area D: CP uses backstabs to attack predators, analogous to a one-dimensional inelastic collision, and the calculation formula is as follows: In the formula, x t CP is the optimal solution, namely CP; x t i represents the position of the i-th individual at iteration t, i.e., the predator; β is the convergence speed factor, τ4 is a random value in [0,1]; F i t To affect the attack power of the CP of the i-th predator, the relevant calculation formula is as follows: In the formula, is the mass of the ith individual at iteration t, f is the objective function, is the final velocity of the ith individual at time t+1, and is allocated based on the random solution selected from the current population, is the initial velocity of the ith individual at iteration t, τ6 is a random vector in [0,1], ε is a minimum value to avoid being divided by 0, and r4 is a random number in [1,N].

4. The multi-mode hybrid energy storage power allocation method in a wind-solar storage system according to claim 1 is characterized in that: The S3 is specifically: S3.1 Decomposition of P using ICEEMDAN HES S3.1.1 Use empirical mode decomposition to decompose Gaussian white noise and obtain the IMF component of the noise l i(i＝1,2,.,N) , add the original signal x(t) to each component after processing, and get the extended signal x i (t): x i (t)=x(t)+a0l i (15) Where, α0=βstd(x(t)) / std(l1) is the standard deviation of white noise, β is a constant, and std represents the standard deviation operation; S3.1.2 Define the first residual component res1: In the formula, M[x i (t)] represents x i (t) minus the first component IMF1 obtained by EMD decomposition. At this step, the first component IMF1 of the original signal x(t) is as follows: IMF1=x(t)-res1 (17) S3.1.3 Calculate the second residual component res2 and derive IMF2: IMF2=res1-res2 (19) Where, α1 = βstd(r1); S3.1.4 Use the same method to calculate the kth residual and IMF: IMF k =res k-1 -res k (21) In the formula, α k-1 =βstd(r k-1 ),k>1; S3.1.5 Iterate until the maximum number of iterations is met; S3.2 constructs composite entropy as an evaluation index. The composite entropy is composed of the envelope entropy E p and mutual information entropy MI; S3.2.1 Calculation of envelope entropy E p Envelope entropy can reflect whether the frequency characteristics of each IMF are obvious. When the IMF characteristic information is less, the envelope entropy value is larger, otherwise, it is smaller. p It can be expressed as: Where a(j) is the envelope signal obtained after Hilbert modulation of signal x(j); p j is the normalized form of a(j); S3.2.2 Calculation of mutual information entropy MI Mutual information in information theory indicates the degree of correlation between two events and is not easily affected by external factors. The larger the mutual information entropy value MI, the stronger the correlation between the two events. As far as IMF is concerned, the richer the original signal feature information contained is. The calculation formula is as follows: MI(X,Y)=H(Y)-H(Y|X) (24) In the formula, X and Y are different events, H(Y) is the entropy of Y, and X is P HES , Y is IMF, H(Y|X) is the conditional entropy of Y when X is known; S3.2.3 Establishing a composite fitness function The expression of the composite fitness function is as follows: ICEEMDAN improves the stability and robustness of the results by repeatedly adding noise and performing decomposition; S3.2.4 Select Nstd and Ne as optimization objects, use compound entropy as fitness function, and use ICPO to optimize the parameter combination [Nstd, Ne] to improve ICEEMDAN's P HES Decomposition effect.

5. The multi-mode hybrid energy storage power allocation method in a wind-solar storage system according to claim 1 is characterized in that: The S4 is specifically: P HES After ICPO-ICEEMDAN decomposition, the spectrum is obtained by fast Fourier transform; S4.1 Calculation of relative entropy D Considering the relative entropy between the IMF spectra, the relative entropy between the IMF spectra is calculated to characterize the degree of spectrum differentiation. The IMF is divided into n combinations. The relative entropy D is calculated as follows: In the formula, x is the amplitude and N is the number of time points; S4.2 Setting power distribution point The maximum value of relative entropy D is taken as the high frequency component P H and low frequency component P L The dividing point K is obtained from the dividing point K, and after the high-frequency area and the low-frequency area are obtained, the p power types and q energy types in the HES selection results are further allocated; For the high-frequency area, the relative entropies are arranged in descending order, and the first p-1 points with larger relative entropy values ​​are taken as power points to divide the high-frequency area into p segments. The same is true for the low-frequency area. The energy storage with a larger rated cycle number is selected as the bearer of the power in the higher frequency band; the power division point setting of M IMFs is summarized as follows: the maximum relative entropy value is used as the dividing point, the left side contains K IMF powers, which is the high-frequency area and is borne by the power-type energy storage; the right side contains MK IMF powers, which is the low-frequency area and is borne by the energy-type energy storage; The decomposed residual res is placed in the low-frequency area and borne by the energy type with a larger number of cycles; S4.3 Reconstruct the power of each mode of energy storage The calculation formula of each mode energy storage power is as follows: Where x, y, x', y' are the relative entropy division points.

Images