An Optimal Wavelet Packets Basis Method for Cascade Hydro-PV-Pumped Storage Generation Systems to Smooth Photovoltaic Power Fluctuations

: Due to the volatility of natural resources, the power ﬂuctuations of photovoltaic (PV) generation have serious negative impacts on the power quality. This paper reports a novel approach to resolve this problem in a combined cascade hydro-PV-pumped storage (CH-PV-PS) generation system through appropriate power distribution on a small time scale. According to the characteristics of power generation systems and multi-constraints, how to obtain the optimal smoothing e ﬀ ects with the small cost is a key challenge. For this purpose, wavelet packet decomposition is modiﬁed by constructing cost function and ideal power trajectory sequences to obtain a new adaptive power distribution method for the CH-PV-PS generation system considering the boundary conditions of the system in this paper. At the same time, to meet the real-time control requirements of actual systems, an additive function is presented to improve the optimization speed of the adaptive power distribution method. In the meantime, a fuzzy controller is designed to optimize the inevitable delay time of the power scheduling system, and the starting threshold is set to avoid the frequently mode conversion of pumped storage. Finally, the performance of the method is evaluated based on the PV station data in Xiaojin County, Sichuan Province, China. Simulation results based on Matlab-simulink indicate that the proposed method can e ﬀ ectively suppress PV power ﬂuctuations and ensure normal operation of the CH-PV-PS generation system within the multiple constraints. which leads to poor following characteristics of the actual power and reference power value of CHSs and VSPSS. The average


Introduction
With global warming and the depletion of fossil energy, renewable energies, including solar energy, wind energy, and hydro energy, have attracted great attention globally [1]. The global annual growth of energy production by far surpasses the growth rate of the total primary energy supply [2]. However, due to the volatility of natural resources, the power fluctuation of PV generation is a great challenge for the stability and security of the entire power network with the continuous increased penetration of PV power stations in the utility grids. Therefore, the power system needs to be equipped with large-capacity controllable energy absorption and compensation devices, so that the power fluctuation of PV generation can be restrained within a certain range to promote the large-scale development of renewable energy [3][4][5][6]. In order to ensure the stability of power grids and reduce the negative impacts of power fluctuation, the Puerto Rico Electric Power Authority and State Grid Corporation of China (1) This paper reports a new approach constructing a combined cascade hydro-PV-pumped storage (CH-PV-PS) generation system for the first time to smooth PV power fluctuations through appropriate power distribution on a small timescale. Combined with the existing pumped storage power station and cascade hydro power station, the CH-PV-PS generation system is helpful in reducing the cost of smoothing PV power fluctuations and improves the utilization of renewable energy. (2) The optimal wavelet packet method is modified, and a fuzzy controller is designed to achieve better smoothing effects. Compared with the existing technologies, this proposed method can adaptively search for the optimal decomposition layer and the optimal biases according to the PV fluctuation characteristics. At the same time, considering the operating boundary conditions of energy absorption and absorption devices, the characteristics of each device can be utilized and each device can be operated in a better operation state. (3) The efficiency of the proposed smoothing method is compared with three different algorithms.
As a result, the proposed smoothing method is more effective in smoothing PV power fluctuations and ensuring normal operation of the CH-PV-PS generation system within the multi-constraints. The rest of the paper is organized as follows. Section 2 introduces the CH-PV-PS generation system. The smoothing control strategy based on the improved optimal wavelet packets basis method proposed in Section 3. In Section 4, the reference power and ideal power trajectory sequence are discussed. Section 5 designs the fuzzy controller to optimize the delay time. Simulations are carried out in Section 6 and conclusions are given in Section 7.

System Description of the CH-PV-PS Generation System
The structures of the CH-PV-PS generation system studied in this paper are shown in Figure 1, including a centralized PV power station 100 MW, CHSs 150 MW, and a VSPSS 5 MW. The public link point of AC bus voltage is 220 KV, and the booster transformer boosts the voltage to 220 KV to send the electricity to the utility grids. The VSPSS is composed of an FSC, a synchronous motor, and a pump turbine. Energy is transmitted between the stator of the motor and the utility grids through FSC. In the power generation mode, the pump turbine operates in the variable-speed mode according to the operation conditions and power requirements. Besides, the FSC converts the electric energy with different voltage frequencies and phases into the same values as the utility grids through an AC/DC/AC converter. Conversely, in the electric mode, the power direction is the opposite, and the output voltage frequency is adjusted by the FSC, to accurately control the pump turbine to absorb power from the grid. Since the FSC completely isolates the motor from the utility grids, the power station can realize pole-less speed regulation, which improves the range and the speed of power regulation to better suppress the fluctuation of active power in the power system [27]. The energy management system, EMS, is the control center of the CH-PV-PS generation system, which is used to sample PV output power, P PV , run algorithms and send reference power commands to CHSs and VSPSS. P HSObj and P HSAct denote the reference power and actual power of CHSs, respectively. P FSC and P PSAct denote the reference power and actual power of VSPSS, respectively. P G denotes the grid output power, and P G = P PV + P HSAct + P PSAct . The rest of the paper is organized as follows. Section 2 introduces the CH-PV-PS generation system. The smoothing control strategy based on the improved optimal wavelet packets basis method proposed in Section 3. In Section 4, the reference power and ideal power trajectory sequence are discussed. Section 5 designs the fuzzy controller to optimize the delay time. Simulations are carried out in Section 6 and conclusions are given in Section 7.

System Description of the CH-PV-PS Generation System
The structures of the CH-PV-PS generation system studied in this paper are shown in Figure 1, including a centralized PV power station 100 MW, CHSs 150 MW, and a VSPSS 5 MW. The public link point of AC bus voltage is 220 KV, and the booster transformer boosts the voltage to 220 KV to send the electricity to the utility grids. The VSPSS is composed of an FSC, a synchronous motor, and a pump turbine. Energy is transmitted between the stator of the motor and the utility grids through FSC. In the power generation mode, the pump turbine operates in the variable-speed mode according to the operation conditions and power requirements. Besides, the FSC converts the electric energy with different voltage frequencies and phases into the same values as the utility grids through an AC/DC/AC converter. Conversely, in the electric mode, the power direction is the opposite, and the output voltage frequency is adjusted by the FSC, to accurately control the pump turbine to absorb power from the grid. Since the FSC completely isolates the motor from the utility grids, the power station can realize pole-less speed regulation, which improves the range and the speed of power regulation to better suppress the fluctuation of active power in the power system [27]. The energy management system, EMS, is the control center of the CH-PV-PS generation system, which is used to sample PV output power, PPV, run algorithms and send reference power commands to CHSs and VSPSS. PHSObj and PHSAct denote the reference power and actual power of CHSs, respectively. PFSC and PPSAct denote the reference power and actual power of VSPSS, respectively. PG denotes the grid output power, and PG = PPV + PHSAct + PPSAct.

Definition of the PV Power Fluctuations
In this study, PV fluctuation is discussed at the minute scale. Before smoothing, PV power fluctuation ∆P bef is defined as the difference between the maximum and minimum of P PV . After smoothing, PV power fluctuation, ∆P aft , is defined as the difference between the maximum and minimum of P G . At time t 0 , ∆P bef and ∆P aft are calculated as follows: (1) where ∆t is the sampling interval and ∆t = 1 s. The predicted PV fluctuation is defined as: where ∆t is the prediction interval, and P pred denotes the predicted PV power. Referring to the "Puerto Rico Electric Power Authority" [28] and the "State Grid Corporation of China", which stipulate that the maximum PV power fluctuation should be less than 10%/min of its installed capacity, the smoothing constraint and the percentage of failure to reach the requirements are defined as: where sum(t) denotes the number of t, which does not satisfy the smoothing constraint; C PV denotes the installed capacity of PV; and ϕ denotes the percentage of failure to meet the requirements. It is well known that the VSPSS power regulation speed is fast, but its capacity is small, which can be used to absorb PV power with a small amplitude and fast fluctuation. The CHSs power regulation speed is slow, but its capacity is large, which can be used to absorb PV power with a large amplitude and slow fluctuation. So, it is necessary to appropriately select the reference smoothing power, which not only meets the utility grid's requirements of the fluctuation rate but also meets the characteristics of VSPSS and CHSs. Therefore, the most difficult factors of smoothing and coordination control are how to distribute system power and consider system boundary conditions. In view of the powerful data processing ability of the wavelet packet decomposition, which can obtain detailed characteristics of the signal, an improved algorithm is researched to solve the above problems.

Smoothing Control Process
The flow chart of the smoothing method is shown in Figure 2. Firstly, PV power fluctuation is calculated. Whether the smoothing algorithm is started or not is decided by PV power fluctuation and the predicted PV power fluctuation. Then, the ideal PV power trajectory sequence is constructed. Referring to the ideal smoothing power trajectory sequence, the short-time PV power, f, is decomposed by the improved optimal wavelet packet basis to obtain the reference smoothing power, P PVObj (t). Then, the ideal power trajectory sequences of CHSs are constructed. Referring to the ideal power trajectory sequences, the difference between P PVObj and f is decomposed by the improved optimal wavelet packet basis to obtain the reference CHSs power, P HSObj (t). Then, the reference power of VSPSS P FSC (t) is calculated. Finally, the delay time is optimized through the fuzzy controller, and the optimized P FSC (t) and P HSObj (t) are sent to the control center as the scheduling instructions. Referring to the ideal smoothing power trajectory sequence, the short-time PV power, f, is decomposed by the improved optimal wavelet packet basis to obtain the reference smoothing power, PPVObj(t). Then, the ideal power trajectory sequences of CHSs are constructed. Referring to the ideal power trajectory sequences, the difference between PPVObj and f is decomposed by the improved optimal wavelet packet basis to obtain the reference CHSs power, PHSObj(t). Then, the reference power of VSPSS PFSC(t) is calculated. Finally, the delay time is optimized through the fuzzy controller, and the optimized PFSC(t) and PHSObj(t) are sent to the control center as the scheduling instructions.

Improved Wavelet Packet Decomposition Algorithm
The decomposition signals with different characteristics can be obtained by selecting a different wavelet packet basis. Therefore, the desired decomposition signals can be obtained by searching for the optimal basis and the optimal layers. Its principle is shown as follows.
Wavelet packet decomposition has a good localization property and fine sampling step size to focus on the details of signals with different frequencies. For any signals, f(t)∈L 2 (R) (square integral function space), the wavelet series can be expressed as: where Cj,k is the dual wavelet integral transformation of f(t) with respect to the orthogonal wavelet function, namely:

Improved Wavelet Packet Decomposition Algorithm
The decomposition signals with different characteristics can be obtained by selecting a different wavelet packet basis. Therefore, the desired decomposition signals can be obtained by searching for the optimal basis and the optimal layers. Its principle is shown as follows.
Wavelet packet decomposition has a good localization property and fine sampling step size to focus on the details of signals with different frequencies. For any signals, f(t)∈L 2 (R) (square integral function space), the wavelet series can be expressed as: where C j,k is the dual wavelet integral transformation of f(t) with respect to the orthogonal wavelet function, namely: where {Ψ j,k (t)} is the Riesz basis of the orthogonal wavelet, Ψ(t), which can be expressed as: Wavelet decomposition can decompose signals at multiple levels to obtain a more precise decomposition method of signal features. In this paper, the orthogonal wavelet function is generated according to the Daubichies (DB) wavelet, which is expressed as: where u n,k (t) is the wavelet packet of the orthogonal scale function, ϕ(t); n∈{2l, 2l + 1}; l = 0, 1, . . . ; h k is the filter coefficient, and its Fourier transform satisfies the following conditions: where N is the number of sampling points of signal f(t). Q(z) satisfies the following conditions: The wavelet filter coefficient, g k = (−1) k−1 h −k+1 , and the decomposed coefficients of the signal, f(t), is d j,n = {d p j,n , p < Z}, where d p j,n is determined by f and u n , k . The discrete data of f(t) is {F k }; k = 0, 1, 2, . . . , N, and its decomposition algorithm is expressed as: Taking two-layer wavelet decomposition and N = 8 as an example, the sampling sequence is decomposed by the wavelet packet according to Equation (11), and its hierarchical structure is shown in Table 2. Table 2. Hierarchical structure of wavelet packet decomposition.
Under an orthogonal wavelet packet basis, the short-time PV output power, f(t), can be decomposed to n-layer wavelet packets, as shown in Figure 3. f(t) corresponds to a decomposition series, x = {x 1 , x 2 , ...,x k }. Here, x k represents the kth decomposition signal, and k = 2N − 1.   Figure 4, which indicate that x4 sequences decomposed by wavelet packet basis DB2 are more suitable as the reference smoothing power than the sequences decomposed by wavelet packet basis DB4. Therefore, for a given signal, it is hoped that a better wavelet packet basis will be selected and the signal decomposed to get the best reference power value for smoothing PV power fluctuation. With the best reference power value, the CH-PV-PS generation system could obtain ideal operation states and smoothing effects. In order to select a better wavelet package basis and the optimal layers, the cost function of a decomposed sequence is firstly defined. Then, to minimize the cost function and obtain the optimal decomposition layer, the bases are searched from the wavelet packet library. For a given signal, the wavelet package basis with the minimum cost function value is the most effective decomposition method, which is called the optimal wavelet package basis. An information cost function, M(x), is defined based on series x = {x1, x2, … xk}. In order to obtain the optimal wavelet packet basis and meet the requirements of real-time engineering control, the information cost function, M(x), should satisfy the following two conditions: (1) Additive conditions as: The wavelet library {u n,k (t)} can form a set of orthogonal basis of L 2 (R) called a wavelet packet basis of L 2 (R), as {u n,k (t)}c. It is assumed that there are h wavelet packet bases in the wavelet packet library, and different wavelet packet basis have different time-frequency localization capabilities, reflecting different signal characteristics. The decomposition effects of different wavelet packet basis are shown in Figure 4, which indicate that x 4 sequences decomposed by wavelet packet basis DB2 are more suitable as the reference smoothing power than the sequences decomposed by wavelet packet basis DB4. Therefore, for a given signal, it is hoped that a better wavelet packet basis will be selected and the signal decomposed to get the best reference power value for smoothing PV power fluctuation. With the best reference power value, the CH-PV-PS generation system could obtain ideal operation states and smoothing effects. In order to select a better wavelet package basis and the optimal layers, the cost function of a decomposed sequence is firstly defined. Then, to minimize the cost function and obtain the optimal decomposition layer, the bases are searched from the wavelet packet library. For a given signal, the wavelet package basis with the minimum cost function value is the most effective decomposition method, which is called the optimal wavelet package basis.   Figure 4, which indicate that x4 sequences decomposed by wavelet packet basis DB2 are more suitable as the reference smoothing power than the sequences decomposed by wavelet packet basis DB4. Therefore, for a given signal, it is hoped that a better wavelet packet basis will be selected and the signal decomposed to get the best reference power value for smoothing PV power fluctuation. With the best reference power value, the CH-PV-PS generation system could obtain ideal operation states and smoothing effects. In order to select a better wavelet package basis and the optimal layers, the cost function of a decomposed sequence is firstly defined. Then, to minimize the cost function and obtain the optimal decomposition layer, the bases are searched from the wavelet packet library. For a given signal, the wavelet package basis with the minimum cost function value is the most effective decomposition method, which is called the optimal wavelet package basis. An information cost function, M(x), is defined based on series x = {x1, x2, … xk}. In order to obtain the optimal wavelet packet basis and meet the requirements of real-time engineering control, the information cost function, M(x), should satisfy the following two conditions: (1) Additive conditions as: (1) Additive conditions as: where M c ({x 1 ,x 2 ...x k }) is the information cost function value of signal decomposed by the wavelet packet basis {u m,k (t)} c . The additive condition indicates that the information cost function value of the signal decomposed by the wavelet package basis {u m,k (t)} c is equal to the superposition of the subsequence cost function value. The cost function with the additive conditions will reduce the computation amount and improve the computation speed of the algorithm.
(2) In order to search for the optimal basis and find the optimal decomposition layer number, the value of the information cost function, M(x), should reflect the degree of convergence between the wavelet packet decomposition signals and the desired decomposition signals. Therefore, it is necessary to define an ideal reference smoothing power trajectory sequence, x e , and an ideal CHSs power trajectory sequence as references, x e = {x e1 , x e2 , ..., x em } (whose calculation procedure is shown in Section 4.1.1). We take m discrete points from signal x n decomposed by the wavelet packet basis at equal distance to get a new decomposition sequence: and the information cost function is defined as: where Q n reflects the degree of convergence between the wavelet packet decomposition signals and the ideal power signals and Q p denotes the energy dissipation function, which reflects the energy loss degree of the decomposition signals. For a period of time, the smaller Q p represents the smaller total active output of VSPSS and CHSs. Let: Obviously, (ln 2 ) is positively correlated with Q c , and it will be replaced by Qc in Equation (15). Therefore, M(x n ) can be expressed as: Energies 2019, 12, 4642 10 of 22 The additive function is defined as: Since the ideal power trajectory stays the same in a certain period of time, m m i = 1 x 2 e is a positive constant. When λ(x n ) is the smallest, M(x n ) is the smallest, and the cost function is converted from Equation (15), which is used to calculate the whole sequence, to Equation (18), which is used to calculate point by point. When a new signal is added, the new signal and the original signal form a new sequence. The cost function value of the new sequence is equal to the cost function value of the original signal plus the cost function value of the new signal, which avoids repeating the solving process and reduces the amount of computation. Therefore, the three categories of signals mentioned in the introduction can be obtained through the above cost function and decomposition algorithm. Furthermore, we can achieve satisfactory smoothing effects.

The Ideal Smoothing Power Trajectory Sequence
According to the grid-connection requirements of PV power stations, the PV power fluctuation within one minute shall not exceed 10% of its rated capacity [29]. In this paper, PV power fluctuation within one minute shall not exceed 10% as the smoothing constraint condition. Taking a PV power station with a rated capacity, P N , as an example, the ideal smoothing power trajectory sequence is studied.
The short-time PV power, P PV , is sampled with T S as the period to get the sampling sequence, x s = {x s1 , x s2 , ..., x sn }. Here, the PV power fluctuation per minute must not exceed 10% of its rated capacity as the constraint condition, and the maximum change rate is expressed as r max = 0.1P N /60. The constraint of the ramping rate is defined as: r PV < r max . During the (i − 1)T s -iT s , the ideal smoothing power trajectory sequence function is calculated as follows: where i = 1, 2, ... n, and L denotes the smoothing coefficient, which can be adjusted to change the smoothing effects, 0 < L < 1. Each time a P i (t) is calculated, x si+1 = P i (iT s ) is re-assigned, and the x s sequence is refreshed at the same time. During the 0-nT s , the ideal smoothing power trajectory sequence function is expressed as: where M discrete points of P(t) are taken at equal distance to get the ideal smoothing power trajectory sequence as:

The Reference PV Smoothing Power
Referring to the ideal smoothing power trajectory sequence, the reference smoothing power is decomposed by the improved optimal wavelet basis as follows: Firstly, the c wavelet packet base {u m,k (t)} c (c = 1, 2, ..., h) is selected to decompose the PV power signal, f(t), to the first layer, and the cost function value of the low-frequency wavelet packet decomposition in the first layer is calculated. Here, the upper nodes are called the parent node and the lower nodes are called the child node. Furthermore, the upper low-frequency signal is decomposed to the next layer, and the cost functions of the parent node and child node are calculated and compared. If the child node's cost function value is less than the parent node's cost function value, the decomposition continues. Otherwise, the cost function values of all low-frequency signals are summed up as λ(x n ) c . The group of wavelet packet bases that makes λ(x n ) c the smallest is called the optimal basis. Under the optimal wavelet basis decomposition, the parent node of the last layer is taken as the reference smoothing power, denoted as P PVObj (t). Then, the ramping rate constraint is introduced to revise the reference smoothing power, which is P Pvobj (t) = {P Pvobj (t), r PV < r max }. During the whole decomposition process, only the low-frequency part of the signal is decomposed, so the operation speed of the algorithm is greatly improved.

Ideal Power Trajectory Sequence of CHSs
CHSs have an optimal operation efficiency according to the working conditions, and the power of CHSs under the optimal operating efficiency is defined as P Hmax . When CHSs participate in PV power smoothing, the ideal power trajectory is defined by the quadratic function under the constraint of the ramping rate, {R CH < b, R CH1 < c}. The maximum climbing rate of its regulating power shall not exceed b. Here, b = 0.1P Hmax /60. The magnitude of the second derivative should not exceed c. Here, c = b/60. Taking the output power of the CHSs with P Hmax as an example under a certain working condition with the optimal operating efficiency, the ideal power trajectory sequence of CHSs is defined.
Firstly, we calculate the power trend of CHSs P HSObj = P PVObj − P PV + P Hmax . Then, P HSObj is sampled to obtain the sampling sequence, x sh = {x sh1 , x sh2 ..., x shn }, whose sampling period is T S . Obviously, x sh represents the power variation trend of CHSs. Then, during (i−1)T s -iT s the ideal power trajectory sequence function is defined as: where i = 1, 2, ..., n. Each time P hi (t) is calculated, x shi+1 = P hi (iT s ) is re-assigned, and the x sh sequence is refreshed. The ideal power trajectory sequence function of CHSs during 0-(n−1)T S is calculated as: where M discrete points of P h (t) are taken at equal distance to get the sequence as: Then, the ideal power trajectory sequence of CHSs is defined as:

The Reference Power of CHSs and VSPSS
In order to ensure that the VSPSS power does not exceed its rated capacity, it is necessary to re-define the information cost function and the additive function. So, we redefine the information cost function as: The additive function is redefined as: where N PSC is the time proportion when the actual operation power of VSPSS exceeds its rated capacity, and k is the penalty coefficient. Then, let: Referring to the ideal power trajectory sequence and the new cost function, f 1 is decomposed by the improved optimal wavelet basis. The reference CHSs power under the optimal wavelet packet basis decomposition is obtained, which is the parent node of the last layer defined as P HS_Obj (t). In the actual operation, in order to ensure the qualified operation efficiency of CHSs, its output reference power is defined as: P HSObj (t) = P Hmax (t) + P HS_Obj (t), where P Hmax (t) is the current optimal operation power of CHSs. Then, the ramping rate constraint is introduced to revise the CHSs reference power, P HSObj (t) = {P HSObj (t), r PV (t) < r max (t)}. The reference output power of VSPSS is obtained as: In order to avoid an excessively fast conversion of the VSPSS mode between pump and power generation, conversion constraints are introduced as follows: P Fsc (t) = (P Fsc (t), 0) (P Fsc (t) ≥ P CLmt (t), P Fsc (t) < P CLmt (t)) .
(31) P FSC (t) under P CLmt is set to 0 MW; the conversion threshold, P CLmt , is defined as: where σ is the inertia factor; P Cbase is the conversion base power; δ is the penalty coefficient; and N c (t 0 -t) denotes the number of conversions during t 0 -t. Considering the capacity constraint of VSPSS, its final reference power instruction is: We set P FSC (t) above C PS to C PS MW, where C PS is the rated capacity of VSPSS. The reference power distribution principle of CHSs and VSPSS is shown as Figure 5. Due to part of the PV power fluctuations with high frequency being absorbed by the VSPSS and part of PV power fluctuation with low frequency being absorbed by CHSs, the output power of the CH-PV-PS generation is smoothed. Then, various constraints ensure the normal operation of the system and improve the service life of the system.

Fuzzy Control to Optimize the Delay Time
In practical engineering applications, because the dates of PV output power are obtained point by point through real-time sampling, the smooth algorithm can only deal with the current period dates and previous dates. For example, the wavelet packet is used to decompose and reconstruct the data within the current 5 minutes. By setting a certain delay time, d, the smoothing target power in the future d period can be obtained, which is the wavelet packet decomposition signal of the last d period. In the other words, it will be used as the target power value during the future d time. In this way, the smoothing target power values of each period are obtained, as shown in Figure 6. It is assumed that the energy absorption and compensation devices can fully follow the target instruction, and the smoothing effects under different delay times are shown in Figure 7.
It can be seen from Figure 6 and Figure 7 that the smoothing algorithm has obvious delay effects. Moreover, the value of the delay time has a great influence on the smoothing effects. Within a certain range, the smoothing effects increase with the decrease of the delay time. However, if the delay time is too short, it will increase the difficulty of calculation and make the operation of VSPSS and CHSs more difficult. According to the current fluctuations, a certain starting threshold can be set to reduce the delay effects on the smoothing algorithm. In addition, a fuzzy controller is designed to optimize the delay time in this paper. Fuzzy control is an effective tool to solve the effects of time delay through language control rules under conditions without an accurate mathematical model of the controlled object [28]. Moreover, according to the original fluctuation and fluctuation after smoothing, a better delay time is set to improve the smoothing effects. Due to part of the PV power fluctuations with high frequency being absorbed by the VSPSS and part of PV power fluctuation with low frequency being absorbed by CHSs, the output power of the CH-PV-PS generation is smoothed. Then, various constraints ensure the normal operation of the system and improve the service life of the system.

Fuzzy Control to Optimize the Delay Time
In practical engineering applications, because the dates of PV output power are obtained point by point through real-time sampling, the smooth algorithm can only deal with the current period dates and previous dates. For example, the wavelet packet is used to decompose and reconstruct the data within the current 5 min. By setting a certain delay time, d, the smoothing target power in the future d period can be obtained, which is the wavelet packet decomposition signal of the last d period.
In the other words, it will be used as the target power value during the future d time. In this way, the smoothing target power values of each period are obtained, as shown in Figure 6. It is assumed that the energy absorption and compensation devices can fully follow the target instruction, and the smoothing effects under different delay times are shown in Figure 7.
It can be seen from Figures 6 and 7 that the smoothing algorithm has obvious delay effects. Moreover, the value of the delay time has a great influence on the smoothing effects. Within a certain range, the smoothing effects increase with the decrease of the delay time. However, if the delay time is too short, it will increase the difficulty of calculation and make the operation of VSPSS and CHSs more difficult. According to the current fluctuations, a certain starting threshold can be set to reduce the delay effects on the smoothing algorithm. In addition, a fuzzy controller is designed to optimize the delay time in this paper. Fuzzy control is an effective tool to solve the effects of time delay through language control rules under conditions without an accurate mathematical model of the controlled object [28]. Moreover, according to the original fluctuation and fluctuation after smoothing, a better delay time is set to improve the smoothing effects.  Table 3. The designs of the membership function and fuzzy rules are as follows: When the original fluctuation and the fluctuation after smoothing are both large, we take the smaller delay time. When the original fluctuation and the fluctuation after smoothing are both small, we take the larger delay time. The priority level of the fluctuation after smoothing is larger than the original fluctuation, and the weighted average method is designed to remove fuzziness.  Table 3. The designs of the membership function and fuzzy rules are as follows: When the original fluctuation and the fluctuation after smoothing are both large, we take the smaller delay time. When the original fluctuation and the fluctuation after smoothing are both small, we take the larger delay time. The priority level of the fluctuation after smoothing is larger than the original fluctuation, and the weighted average method is designed to remove fuzziness.  Figure 8, and the fuzzy control rules are shown in Table 3. The designs of the membership function and fuzzy rules are as follows: When the original fluctuation and the fluctuation after smoothing are both large, we take the smaller delay time. When the original fluctuation and the fluctuation after smoothing are both small, we take the larger delay time. The priority level of the fluctuation after smoothing is larger than the original fluctuation, and the weighted average method is designed to remove fuzziness.  Table 3. The designs of the membership function and fuzzy rules are as follows: When the original fluctuation and the fluctuation after smoothing are both large, we take the smaller delay time. When the original fluctuation and the fluctuation after smoothing are both small, we take the larger delay time. The priority level of the fluctuation after smoothing is larger than the original fluctuation, and the weighted average method is designed to remove fuzziness.

Simulation and Discussion
Experimental studies were carried out to verify the comprehensive effects of the proposed method on the CH-PV-PS generation system, which was composed of a 100 MW PV power station, the 150 MW CHSs, and a 5 MW VSPSS. The output powers of the CH-PV-PS generation system measured in extreme typical weather were used to build a generation system simulation platform.
The optimal wavelet packet decomposition algorithm (OWPDA) proposed in this paper was compared with the filtering algorithm (FA) [30] and the wavelet packet decomposition algorithm (WPDA) [26] to verify its superiority and feasibility. The influence of CHSs on the smoothing effects was assessed by four comparison experiments, which were respectively implemented on the system with CHSs and without CHSs. The fuzzy controller was designed to optimize the delay time of all the algorithms. The smoothing targets refers to the standard "Puerto Rico Electric Power Authority" [29] and "technical regulations on PV access to utility grids", which stipulates that the active power change of large-scale PV power stations should not exceed 10% of the installed capacity per minute and should not exceed 1/3 of the installed capacity within 10 min. Maximum fluctuation, average fluctuation per minute, and the percentage of failure to reach the requirements were calculated to verify the smoothing effects. Average cumulative errors were calculated to verify the rationality of the system power distribution. Figure 9 shows the smoothing effects of FA without CHSs. The maximum fluctuation before smoothing was 19.08 MW, the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 15.41 MW, and the average fluctuation after smoothing was 5.34 MW, and the percentage of failure to reach the requirements after smoothing was 12.12%. The average cumulative error per minute between the reference power and the actual power of VSPSS was 35.75 MW.

Simulation and Discussion
Experimental studies were carried out to verify the comprehensive effects of the proposed method on the CH-PV-PS generation system, which was composed of a 100 MW PV power station, the 150 MW CHSs, and a 5 MW VSPSS. The output powers of the CH-PV-PS generation system measured in extreme typical weather were used to build a generation system simulation platform.
The optimal wavelet packet decomposition algorithm (OWPDA) proposed in this paper was compared with the filtering algorithm (FA) [30] and the wavelet packet decomposition algorithm (WPDA) [26] to verify its superiority and feasibility. The influence of CHSs on the smoothing effects was assessed by four comparison experiments, which were respectively implemented on the system with CHSs and without CHSs. The fuzzy controller was designed to optimize the delay time of all the algorithms. The smoothing targets refers to the standard "Puerto Rico Electric Power Authority" [29] and "technical regulations on PV access to utility grids", which stipulates that the active power change of large-scale PV power stations should not exceed 10% of the installed capacity per minute and should not exceed 1/3 of the installed capacity within 10 min. Maximum fluctuation, average fluctuation per minute, and the percentage of failure to reach the requirements were calculated to verify the smoothing effects. Average cumulative errors were calculated to verify the rationality of the system power distribution. Figure 9 shows the smoothing effects of FA without CHSs. The maximum fluctuation before smoothing was 19.08 MW, the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 15.41 MW, and the average fluctuation after smoothing was 5.34 MW, and the percentage of failure to reach the requirements after smoothing was 12.12%. The average cumulative error per minute between the reference power and the actual power of VSPSS was 35.75 MW.  By analyzing the above dates in combination with Figure 9, some reasons why the smoothing effects of FA is poor in the system as shown in Figure 9b can be concluded as follows:

Filtering Algorithm
Due to the influence of climate, PV power fluctuation has different characteristics, and the uncertainty in the filtering time constant means that the fluctuation of the reference smoothing power is unqualified as shown in Figure 9a.
In the process of system power distribution, the dynamic response of the speed and maximum power (5 MW) of VSPSS are not considered, which increases the average cumulative error per minute between the reference power and the actual power of VSPSS as shown in Figure 9c, and has negative effects on the service life of the generation system.
When the PV power fluctuates greatly, the capacity of VSPSS is insufficient as shown in Figure  9c. In order to improve the smoothing effects, the installed capacity of VSPSS should be increased or additional energy absorption and compensation devices should be introduced. Figure 10 shows the smoothing effects of FA after the introduction of CHSs. According to the operation conditions, the optimal power value of CHSs is assumed to be 60 MW. The maximum fluctuation before smoothing is 19.08 MW, the average fluctuation before smoothing is 10.34 MW, and the percentage of failure to reach the requirements before smoothing is 24.84%. The maximum fluctuation after smoothing is 15.87 MW, and the average fluctuation after smoothing is 4.02 MW, and the percentage of failure to reach the requirements after smoothing is 9.39%. After the introduction of CHSs, due to the fluctuations with large amplitude and slow frequency being absorbed by CHSs, as shown in Figure 10c, the smoothing effects are improved, as shown in Figure  10b. However, the system power distribution based on FA does not consider the characteristics and operation boundary conditions of the system as shown in Figure 10c By analyzing the above dates in combination with Figure 9, some reasons why the smoothing effects of FA is poor in the system as shown in Figure 9b can be concluded as follows: Due to the influence of climate, PV power fluctuation has different characteristics, and the uncertainty in the filtering time constant means that the fluctuation of the reference smoothing power is unqualified as shown in Figure 9a.
In the process of system power distribution, the dynamic response of the speed and maximum power (5 MW) of VSPSS are not considered, which increases the average cumulative error per minute between the reference power and the actual power of VSPSS as shown in Figure 9c, and has negative effects on the service life of the generation system.
When the PV power fluctuates greatly, the capacity of VSPSS is insufficient as shown in Figure 9c. In order to improve the smoothing effects, the installed capacity of VSPSS should be increased or additional energy absorption and compensation devices should be introduced. Figure 10 shows the smoothing effects of FA after the introduction of CHSs. According to the operation conditions, the optimal power value of CHSs is assumed to be 60 MW. The maximum fluctuation before smoothing is 19.08 MW, the average fluctuation before smoothing is 10.34 MW, and the percentage of failure to reach the requirements before smoothing is 24.84%. The maximum fluctuation after smoothing is 15.87 MW, and the average fluctuation after smoothing is 4.02 MW, and the percentage of failure to reach the requirements after smoothing is 9.39%. After the introduction of CHSs, due to the fluctuations with large amplitude and slow frequency being absorbed by CHSs, as shown in Figure 10c, the smoothing effects are improved, as shown in Figure 10b. However, the system power distribution based on FA does not consider the characteristics and operation boundary conditions of the system as shown in Figure 10c,d, which leads to poor following characteristics of the actual power and reference power value of CHSs and VSPSS. The average cumulative error of CHSs and VSPSS per minute are both large, which are 24.09 and 8.21 MW, respectively. Therefore, the smoothing effects are poor. cumulative error of CHSs and VSPSS per minute are both large, which are 24.09 and 8.21 MW, respectively. Therefore, the smoothing effects are poor.  Figure 11 shows the smoothing effects of WPDA without CHSs. The maximum fluctuation before smoothing was 19.08 MW, and the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 14.73 MW, and the average fluctuation after smoothing was 4.10 MW, and the percentage of failure to reach the requirements after smoothing was 11.8%. The average cumulative error per minute between the reference power and the actual power of VSPSS was 47.93 MW. Compared with the FA, WPDA can obtain more signal features, and the reference smoothing power is more reasonable. However, due to the limited capacity (5 MW) of VSPSS, the PV power fluctuations with large amplitudes were not absorbed as shown in Figure 11c, which leads to poor smoothing effects as shown in Figure 11b.  Figure 11 shows the smoothing effects of WPDA without CHSs. The maximum fluctuation before smoothing was 19.08 MW, and the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 14.73 MW, and the average fluctuation after smoothing was 4.10 MW, and the percentage of failure to reach the requirements after smoothing was 11.8%. The average cumulative error per minute between the reference power and the actual power of VSPSS was 47.93 MW. Compared with the FA, WPDA can obtain more signal features, and the reference smoothing power is more reasonable. However, due to the limited capacity (5 MW) of VSPSS, the PV power fluctuations with large amplitudes were not absorbed as shown in Figure 11c, which leads to poor smoothing effects as shown in Figure 11b.  Figure 11 shows the smoothing effects of WPDA without CHSs. The maximum fluctuation before smoothing was 19.08 MW, and the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 14.73 MW, and the average fluctuation after smoothing was 4.10 MW, and the percentage of failure to reach the requirements after smoothing was 11.8%. The average cumulative error per minute between the reference power and the actual power of VSPSS was 47.93 MW. Compared with the FA, WPDA can obtain more signal features, and the reference smoothing power is more reasonable. However, due to the limited capacity (5 MW) of VSPSS, the PV power fluctuations with large amplitudes were not absorbed as shown in Figure 11c, which leads to poor smoothing effects as shown in Figure 11b.  Figure 12 shows the smoothing effects of WPDA after introducing CHSs. The maximum fluctuation before smoothing was 19.08 MW, and the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 12.79 MW, and the average fluctuation after was 3.64 MW, and the percentage of failure to reach the requirements after smoothing was 7.27%. The average cumulative error per minute between the actual power value and the reference power value of VSPSS was 0.98 MW, which has a good following performance, and that of CHSs was 10.64 MW. As shown in Figure 12d, the fluctuation speed of CHSs reference power is too fast, which leads to a poor following performance. Compared with FA, WPDA distributes system power according to the response bandwidth of the CH-PV-PS generation system, so the reference smoothing power is more reasonable, as shown in Figure 12a, and the smoothing effects are improved, as shown in Figure 12b. However, due to the system response bandwidth varying with the external operation conditions, the fixed orthogonal basis of decomposition, fixed decomposition layers, and fixed frequency band are not reasonable, and WPDA fails to consider the system boundary condition, as shown in Figure 12c and Figure 12d, which not only affects the smoothing effects but also reduces the efficiency and service life of the system, so this smoothing method applied to the system is unsatisfactory.  Figure 12 shows the smoothing effects of WPDA after introducing CHSs. The maximum fluctuation before smoothing was 19.08 MW, and the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 12.79 MW, and the average fluctuation after was 3.64 MW, and the percentage of failure to reach the requirements after smoothing was 7.27%. The average cumulative error per minute between the actual power value and the reference power value of VSPSS was 0.98 MW, which has a good following performance, and that of CHSs was 10.64 MW. As shown in Figure 12d, the fluctuation speed of CHSs reference power is too fast, which leads to a poor following performance. Compared with FA, WPDA distributes system power according to the response bandwidth of the CH-PV-PS generation system, so the reference smoothing power is more reasonable, as shown in Figure 12a, and the smoothing effects are improved, as shown in Figure 12b. However, due to the system response bandwidth varying with the external operation conditions, the fixed orthogonal basis of decomposition, fixed decomposition layers, and fixed frequency band are not reasonable, and WPDA fails to consider the system boundary condition, as shown in Figure 12c,d, which not only affects the smoothing effects but also reduces the efficiency and service life of the system, so this smoothing method applied to the system is unsatisfactory.  Figure 13 shows the smoothing effects of the improved OWPDA after introducing CHSs. Similarly, the optimal power of CHSs was assumed to be 60 MW according to the operation conditions. Because the optimal operating efficiency was considered, the CHSs always work near the optimal efficiency point as shown in Figure 13d. The maximum fluctuation before smoothing was 19.08 MW, the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 1.90 MW, and the average fluctuation after smoothing was 0.67 MW, and the percentage of failure to reach the requirements after smoothing was 0.3%. The actual power value and reference power value of VSPSS, and the actual power value and reference power value of CHSs both have excellent following characteristics as shown in Figure 13c and Figure 13d, and the average cumulative errors per minute are 1.17 and 2.38 MW, respectively. 6.3. The Improved Optimal Wavelet Packet Basis Decomposition Algorithm Figure 13 shows the smoothing effects of the improved OWPDA after introducing CHSs. Similarly, the optimal power of CHSs was assumed to be 60 MW according to the operation conditions. Because the optimal operating efficiency was considered, the CHSs always work near the optimal efficiency point as shown in Figure 13d. The maximum fluctuation before smoothing was 19.08 MW, the average fluctuation before smoothing was 10.34 MW, and the percentage of failure to reach the requirements before smoothing was 24.84%. The maximum fluctuation after smoothing was 1.90 MW, and the average fluctuation after smoothing was 0.67 MW, and the percentage of failure to reach the requirements after smoothing was 0.3%. The actual power value and reference power value of VSPSS, and the actual power value and reference power value of CHSs both have excellent following characteristics as shown in Figure 13c,d, and the average cumulative errors per minute are 1.17 and 2.38 MW, respectively. The comparisons of the performance indexes are shown in Table 4. Compared with the performance indexes of FA and WPDA, the improved OWPDA gives full play to the characteristics of CHSs and VSPSS. As shown in Figure 13c and Figure 13d, PV power with large amplitude and slow fluctuation is absorbed by CHSs and PV power with small amplitude and fast fluctuation is absorbed by VSPSS, therefore, the system power distribution is more reasonable. Due to the improved method selecting the decomposition layer number and decomposition basis adaptively by referring to the ideal smoothing power trajectories, the reference smoothing power is more reasonable as shown in Figure 13a. Moreover, this method also takes into account the operation boundary conditions of the generation system. Therefore, it not only significantly improves the smoothing effects, as shown in Figure 13b, but also improves the operation efficiency and service life, and ensures normal operation of the CH-PV-PS system.  The comparisons of the performance indexes are shown in Table 4. Compared with the performance indexes of FA and WPDA, the improved OWPDA gives full play to the characteristics of CHSs and VSPSS. As shown in Figure 13c,d, PV power with large amplitude and slow fluctuation is absorbed by CHSs and PV power with small amplitude and fast fluctuation is absorbed by VSPSS, therefore, the system power distribution is more reasonable. Due to the improved method selecting the decomposition layer number and decomposition basis adaptively by referring to the ideal smoothing power trajectories, the reference smoothing power is more reasonable as shown in Figure 13a. Moreover, this method also takes into account the operation boundary conditions of the generation system. Therefore, it not only significantly improves the smoothing effects, as shown in Figure 13b, but also improves the operation efficiency and service life, and ensures normal operation of the CH-PV-PS system.

Conclusions
In this paper, a new power fluctuation smoothing method for the CH-PV-PS generation system was presented to smooth PV power fluctuation, which not only increased the utilization of clean and renewable energy but also had good smoothing effects. PV power fluctuations were completely absorbed and compensated, which resulted in complementary energy advantages. To solve the problem of how to allocate the system power more reasonably according to the characteristics of the generation system and PV power fluctuations, an improved OWPDA was proposed. Considering the PV power fluctuation characteristics, system operation boundary, and response speed, the cost function was introduced in the algorithm and the ideal power trajectory sequence was constructed to achieve a more reasonable power distribution. Moreover, considering the real-time control requirements for engineering applications, the starting threshold of the algorithm was set, and a function with an additive characteristic was introduced to improve the optimization speed of the adaptive power distribution method, and a fuzzy controller was designed to optimize the inevitable delay time of the control system. The analysis results indicate that the proposed control algorithm in this paper is better than FA and WPDA, which not only significantly improves the smoothing effects but also improves the operation efficiency and service life of the system. This power fluctuation smoothing method has important engineering application prospects, but the capacity optimization configuration of the CH-PV-PS generation system needs further study.