1. Load fluctuation

The load volatility indicator can reflect stability and fluctuation within 1 day after load optimization. The mathematical formulas are as follows: min F 1 = ∑ t = 1 23 P n e t t + 1 − P n e t t ∑ t = 1 24 P n e t t ⋅ ε n e t , (1) P n e t t = P l o a d t + ∑ m = 1 N M η c P S , t c t − ∑ m = 1 N M η d P S , t d t , (2) where F 1 and ε n e t are the penalty cost and penalty coefficient for load fluctuation, respectively; P n e t t + 1 is the load value at time t + 1 ; P n e t t is the load value at time t ; P S , t c is the charging power of the energy storage plant at time t ; P S , t d is the discharging power of the energy storage plant at time t ; η c and η d are the charging and discharging efficiency of the energy storage plant, respectively; and N M is the total number of energy storage plants.

2. Energy storage plant costs

During the operation of the energy storage plant, the cost of the energy storage plant consists of two parts, namely, operation and maintenance costs and grid-connected environmental benefits, which are expressed as follows: min F 2 = f s c − f s s = ∑ t = 1 T β s o c ( η d P S t d + η c P S t c ) − ∑ t = 1 T γ s o c η d P d i s , t , (3) where f s c is the operation and maintenance cost of the energy storage plant, f s s is the environmental benefit of the energy storage plant, β s o c is the operating cost coefficient of the energy storage plant, γ s o c is the environmental benefit coefficient generated by the grid-connected energy storage plant, and the meanings of the remaining variables are consistent with the previous paper.

S min ≤ S s t ≤ S max , u t c P S min c ≤ P S t c ≤ u t c P S max c , u t d P S min d ≤ P S t d ≤ u t d P S max d , (4) where S s t is the charge state of the energy storage plant at time t ; S max and S min are the upper and lower limits of the charge state of the energy storage plant, respectively; u t c and u t d are the charge and discharge states of the energy storage plant, respectively; P S max c and P S min c are the maximum and minimum values of the charge power, respectively; and P S max d and P S min d are the maximum and minimum values of the discharge power, respectively.

In the lower model, the economic indicators include the operation and maintenance costs of wind farms, PV plants, and hydropower units; the start-up and shutdown costs of hydropower units; the environmental benefits of grid connection and consumption of wind farms, PV plants, and energy storage systems; and the cost of power purchase and sale of the system. min F 3 = f 1 + f 2 − f 3 + f 4 , f 1 = ∑ t = 1 T μ w P w , t a c c + μ p v P p v , t a c c + μ h y P h y , t a c c , f 2 = ∑ t = 1 T ∑ j = 1 N G s j , t u j , t 1 − u j , t − 1 , f 3 = ∑ t = 1 T β w P p w , t a c c + β p v P p v , t a c c , f 4 = ∑ t = 1 T ρ P C b t P l i n e t Δ t − ρ P C s t P l i n e t Δ t , (5) where F 3 is the system economic index; f 1 is the system operation and maintenance cost; μ w , μ p v , and μ h y are the operating cost coefficients of wind power, PV, and hydropower units, respectively; and P w , t a c c , P p v , t a c c , and P h y , t a c c are the power consumption of wind power, PV, and hydropower units at time t , respectively. f 2 is the start–stop cost of hydropower unit j; s j , t is the start–stop cost of hydropower unit j at time t ; u j , t is the start–stop state of hydropower unit j at time t ; and f 3 is the environmental benefit of wind power, PV, and energy storage consumption. β w and β p v are the environmental benefit coefficients for grid-connected consumption of wind power and PV, respectively; f 4 is the cost of power purchase and sale; ρ P C b t and ρ P C s t are the unit price of power purchase and unit price of power sale to the superior grid at time t , respectively; P l i n e t is the value of exchange power of the contact line at time t ; and Δ t is the interval between the two time sampling.

Tidal risk indicators include the system active network loss rate and voltage vulnerability, and the active network loss rate of the grid is given as the percentage of the grid power loss to feed-in tariff. The grid vulnerability indicator is a state vulnerability indicator that starts with the operational security of the grid. It measures the risk resistance of the grid by analyzing the degree of voltage deviation from standard values at each bus to improve power quality. The higher the vulnerability, the more unstable the voltage at that bus is, and the lower the quality of power supplied, the less the risk-resistance is. The expressions are as follows: min F 4 = f P _ l o s s + f v o l , f P _ l o s s = ∑ t = 1 T ∑ i , j ∈ M L G i j U i , t 2 + U j , t 2 − 2 U i , t U j , t ⁡ cos ⁡ δ i j , t ∑ t = 1 T P w , t a c c + P p v , t a c c + P h y , t a c c + P S , t d , f v o l = 1 24 ∑ t = 1 24 1 2 B V t + 1 2 J t , B V t = 1 N ∑ i = 1 N V t , i , V t , i = v t , i − v t , i min v t , i max − v t , i min , v t , i = U t , i U i , o − 1 0.07 , J t = 1 − ∑ i = 1 N p t , i log 2 1 p t , i log 2 ⁡ N 2 π , p t , i = V t , i ∑ i = 1 N V t , i , (6) where F 4 is the tidal risk indicator; f P _ l o s s is the active network loss rate of the system; f v o l is the grid vulnerability; M L is all branches of the network; G i j is the set of conductances of buses i and j; U i , t and U j , t are the voltages of buses i and j at time t , respectively; δ i j , t is the phase angle difference between buses i and j at the beginning and end of the branch at time t ; B V t is the average vulnerability of the grid at time t ; and N is the total number of system buses. v t , i is the vulnerability of any bus i in the network at time t , and 0.07 is defined as the maximum voltage offset; v t , i max and v t , i min , respectively, are the maximum and minimum vulnerability values of all buses at time t before normalization. V t , i is the normalized value of the vulnerability of each bus in time section t . J t is the vulnerability balance of the network at time t , ranging from 0 to 1, with 0 representing absolute balance and 1 representing absolute imbalance, and p t , i is the ratio of the vulnerability of bus i to the total vulnerability of the network at the current time.

1. Grid tide constraints

To improve the training stability of the original GAN model and the quality of the generated sample data, the paper uses a deep convolutional generative adversarial network (DCGAN) to construct the model. DCGAN uses the strong feature extraction ability of a convolutional neural network (CNN) to improve the quality of the generated sample data of the original GAN model. Compared with the original GAN model, the structure of DCGAN removes the fully connected layer and all pooling layers in the network and uses stepwise convolution instead of the pooling layers of the original GAN model.

1. Generator network structure

It contains four neural network layers, and each network layer contains a deconvolution layer, a batch normalization layer, and an activation function. The deconvolution layer changes the input dimension by using multiple deconvolution kernels. Each layer uses ReLU as the activation function, and the last layer uses tanh as the activation function.

2. Discriminator network structure

The sample changes its dimension after passing through the convolutional layers, and finally, the discriminator’s score for that sample is output through a fully connected layer with output dimension 1. All layers use LeakyReLU as the activation function. The structure of the discriminator designed in this paper is shown in Figure 2.

1. Wasserstein GAN-gradient penalty (WGAN-GP)

The Jensen–Shannon divergence is used in original GAN to measure the similarity of two probability distributions, but during the training process, there may be cases where the scatter value is constant or meaningless, which leads to the disappearance of the gradient. To solve the problem of unstable GAN training, WGAN uses the Wasserstein distance to construct the model instead of the discriminator loss function in the original GAN model. The advantage is that when there is no overlap between two probability distributions, the distance between them can still be effectively described, and the Wasserstein distance can be defined as follows: W P x , P x ′ = 1 K sup f L ≤ K E x ∼ P x f x − E x ′ ∼ P x ′ f x ′ , (14) where sup is the minimum upper bound, x is the real sample data, x ′ is the generated sample data, f is the discriminator using parameters for neural network fitting, K is the Lipschitz constant for the discriminator f , and f L ≤ K is the K-Lipschitz limit satisfied by f .

However, WGAN enforces the Lipschitz continuity condition for weight cropping, which will reduce the training speed and learning efficiency. In order to avoid the WGAN gradient binarization, gradient disappearance, and explosion problems and further enhance the gradient controllability, this paper adds the gradient penalty term to the WGAN loss function to construct WGAN-GP to improve the model training efficiency. The objective function of WGAN-GP is as follows: min G max D V G , D = E x ∼ P x D x − E z ∼ P z D G z + λ E x ^ ∼ P x ^ ∇ x ^ D x ^ 2 − 1 2 , (15) where λ is the gradient penalty term coefficient; ∇ is the gradient operator; x ^ is the random interpolation between the real and generated samples, x ^ = ξ x + 1 − ξ G z ; ξ is the uniform distribution obeying [0,1]; P x ^ is the linear uniform sampling between the real and generated sample sampling points; and ⋅ p is the p -parameter.

2. Conditional GAN (CGAN)

GAN is trained by adding additional labels to the inputs of generators and discriminators. The problem of the unclear and uncontrollable direction of the original GAN generation is solved. The CGAN model is shown in Figure 3. Compared with the original GAN model, CGAN combines noise and additional conditions as inputs on the input side. The additional conditions can be any type of auxiliary information, such as category labels, or other modal data. The objective function of CGAN is as follows: min G max D V D , G = E x ∼ P x D x | c − E z ∼ P z D G z | c , (16) where c is the label of the sample data.

3) Model evaluation metrics

To verify the validity of the model, the average generation effect of the wind and PV scenes was evaluated using the root mean square error (RMSE) and the mean absolute error (MAE). The expressions for RMSE and MAE, respectively, are as follows: R M S E = 1 N ∑ n = 1 N x n ′ − x n 2 , (17) M A E = 1 N x n ′ − x n , (18) where R M S E and M A E are used to measure the error size between the generated samples and the real sample data, respectively (the smaller the value, the higher the accuracy of the generated sample data); N is the total number of scenes generated; and x n and x n ′ are the real samples and the data points corresponding to the generated samples, respectively.

During the population initialization process, the positions of coatis are randomly generated, and the update of individual positions depends on the update of population positions. To expand the effective search space of the algorithm and reduce the search blind spots, a refractive backward learning strategy is used to perturb the population positions. The idea of the refraction inverse learning strategy is derived from convex lens imaging. It extends the search range and enhances the search capability by generating a reverse position from the current coordinates, and its principle is shown in Figure 4.

From the geometric relationship of the line segments in Figure 4, it can be inferred that sin ⁡ α = a + b 2 − x l , sin ⁡ β = x * − a + b / 2 / l * , n = sin ⁡ α / sin ⁡ β , k = l / l * , k n = a + b / 2 − x x * − a + b / 2 , (20) where a , b is the search interval of the solution on the x -axis; l and l * are the lengths of the incident and refracted rays, respectively; and α and β are the angles of incidence and refraction, respectively.

The population location initialization formula for the fusion refraction backward learning strategy is as follows: x i j = l b j + r ⋅ u b j − l b j , x i * j = a j + b j 2 + a j + b j 2 k − x i j k , (21) where x i j is the position of the ith coati in the jth dimension; l b j and u b j are the lower and upper bounds of the decision variables in dimension i , respectively; x i * j is the refractive inverse solution of x i j ; a j and b j are the minimum and maximum values in dimension j on the search space, respectively; and r is a random number obeying uniform distribution in the interval 0,1 .

The COA population location update strategy was based on the coati attacking iguana behavior in Phase I. During this phase, the coati population was assumed to be divided into two parts: one-half of the coatis climbed trees to attack iguanas, and the other waited under the trees for iguanas to fall and preyed on them.

In this paper, we consider the integration of the Levy flight strategy for its location update to improve the local development of tree-climbing coatis. Levy flight is a special type of random wandering process. This form of wandering presents a combination of short-distance exploration and long-distance walking. The location update formula of the improved tree-climbing coatis is as follows: x i t + 1 j = x i t j ⋅ u / v 1 / β + r ⋅ x b s e t t j − I ⋅ x i t j , (22) σ u = Γ 1 + β sin π β 2 Γ 1 + β 2 ⋅ β ⋅ 2 β − 1 / 2 1 / β , (23) where i = 1,2 , … , N 2 , N is the population size, x i t j is the jth dimensional position of the ith coati at the current iteration number, x b e s t t is the best individual position in the current population, I is a random integer in the set 1,2 , t is the current iteration number, and β is a constant, generally set to 1.5, u ∼ N 0 , σ u 2 , v ∼ N 0,1 .

In this paper, a sine–cosine optimization strategy is used to improve the position update of ground coatis and balance the global search and local exploration abilities of COA. It is based on the principle of gradually finding the optimal solution by using the fluctuation properties of the sine and cosine functions. The improved ground coati’s position update equation is as follows: I g j t = l b j + r ⋅ u b j − l b j , (24) x i t + 1 j = x i t j + r 1 ⋅ sin r 2 ⋅ r 3 I g t j − I ⋅ x i t j , r 4 < 0.5 , x i t j + r 1 ⋅ cos r 2 ⋅ r 3 I g t j − x i t j , r 4 ≥ 0.5 , (25) r 1 = a − t × a T , (26) where i = N 2 + 1 , N 2 + 2 , … , N ; I g t j is the current optimal individual; T is the maximum number of iterations; a is a constant, a = 2 ; r 2 and r 3 are random numbers obeying uniform distribution in the interval 0,2 π and − 2,2 , respectively; and r 4 is the selection control factor of the sine and cosine functions, which is the random number in the interval 0,1 .

The optimal solution is updated using Eq. 27 after the first phase of COA. x i t + 1 j = x i t + 1 j , F x i t + 1 j < F x i t j , x i t j , e l s e , (27) where F is the fitness function.

The COA population location update strategy is based on the coati fleeing predator behavior in Phase II. By simulating the predator avoidance strategy, each coati is made to generate a new location randomly near its current location. A spiral search strategy is introduced to avoid COA from falling into a local optimal solution at this phase, and the mathematical formulas are as follows: l b j l o c a l = l b j t , u b j l o c a l = u b j t , (28) x i t + 1 j = x i t j + D ⋅ e b l ⁡ cos 2 π l + 1 − 2 r ⋅ l b j l o c a l + r ⋅ u b j l o c a l − l b j l o c a l , (29) D = 2 ⋅ r ⋅ x b e s t t j − x i t j , (30) where u b j l o c a l and l b j l o c a l are the upper and lower bounds of the jth dimensional decision variables updated with the number of iterations, respectively; D is the distance between the optimal individual and the current individual; and b defines the spiral shape, l ∈ − 1,1 .

Furthermore, the t-distribution variation mechanism is used to perturb the individual positions of coatis to improve the COA finding ability. The t-distribution variational operator combines the Gaussian and Cauchy operators to replace the degree of freedom parameter of the t-distribution with the number of iterations of the algorithm. The t-distribution converges to the Cauchy distribution at the beginning of the algorithm iterations to improve the COA global search capability. As the number of iterations increases, the t-distribution converges to a Gaussian distribution to improve the COA local search level.

The iterative process of t-distribution variation is as follows: x i , n e w t j = x b e s t t j + x b e s t t j ⋅ t i t e r , (31) where x i , n e w t j is the new position of the coati after the perturbation of the t-distribution variation, x b e s t t j is the current optimal individual position of the long-nosed raccoon, and t i t e r is the t-distribution with the number of algorithm iterations i t e r as the degrees of freedom.

The perturbation is updated using Eq. 32 after the perturbation: x i t + 1 j = x i , n e w t j , F x i t + 1 j < F x i t j , x i t + 1 j , e l s e . (32)

The solution flow of the improved COA-based two-tier coordinated optimal scheduling model is shown in Figure 5.

Step 1

Initialization of basic parameters. Initialize the basic data and basic parameters of the algorithm, such as the number of iterations of Improved Coati optimization algorithm and the initial population size.

Step 2

Population initialization. Initialize the coati population, introduce the reflexive learning mechanism, and use the fitness function to evaluate the merits of the individuals generated in Step 1 to find the optimal individual positions.

Step 3

Optimization search in Phase I. Introduce the Levy flight and sine–cosine optimization algorithm in Phase I of ICOA to improve the original ICOA strategy and find the optimal individual position.

Step 4

Optimization search in Phase II. In Phase II of ICOA, fuse the spiral search mechanism and the t-distribution variation mechanism to find the optimal individual position.

Step 5

Upper and lower model interactions. Transfer the charging and discharging strategies and load profiles of the upper-tier energy storage plants to the lower-tier plants. Use the upper-tier transfer results to initialize the lower-tier parameters. Optimize the scheduling strategies of wind power and PV, and feedback the output of wind power, PV, and hydropower units to the upper tier.

Step 6

Iterative optimization. Judge whether the output requirements are met, and if the termination iteration requirements are met, output the optimized scheduling plan before the output day; if not, return to Step 3.

A modified IEEE 30-bus system is used for arithmetic analysis. The system includes four hydropower units, two wind farms, two PV plants, and one energy storage plant. The system structure is shown in Figure 9. The energy storage power station is connected to bus 2, the PV power station is connected to bus 7, and the wind farm is connected to bus 8. The summer solstice day with the maximum system load operation mode is taken as the typical day of load, and the wind power and PV generate the output profile, as shown in Figure 10. The parameters of the upper model are set as follows: the charging and discharging efficiency of both energy storage power plants is 95%; the maximum charging and discharging power is 50 MW; the upper and lower limits of the charge ratio are 0.9 and 0.2, respectively; the rated capacity of the wind farm is 100 MW; and the rated capacity of the PV power plant is 80 MW.

The optimal output strategy of the energy storage system and the system load fluctuation are obtained from the upper-tier optimization. The load profiles before and after optimization are shown in Figure 11. It can be observed that the peak-to-valley difference of the optimized system load profile becomes smaller. The load distribution before and after optimization is provided in Table 5. The difference between peaks and valleys of load changes decreased from the initial 170.04 MW to 140.05 MW, representing a decrease of 17.63%. In particular, the peak load is reduced from 283.4 MW to 261.5 MW, and the valley load is increased from 113.36 MW to 121.45 MW. The peak–valley difference before and after optimization changes from 60% to 53.56%, and the load fluctuation rate decreases from 9.45% to 7.97%. The optimization strategy based on the load profile can effectively reduce the overall load fluctuation. Energy storage plants can smooth the overall load profile of the system, reduce the peak-to-valley difference of the load, and realize “peak shaving and valley filling.”

The charge/discharge and charge state changes in the optimized energy storage system are shown in Figure 12. The scheme controls the charge and discharge states of the energy storage system in order to compensate for the anti-peaking characteristics of wind power. It can reduce the output of the hydropower unit and bring some environmental benefits. As can be observed from Figure 12, the storage system starts charging at low load from 1:00 a.m. to 6:00 a.m. and from 3:00 p.m. to 6:00 p.m. to increase the “valley” of the load profile. In the peak hours of 10:00 a.m. to 2:00 p.m. and 7:00 p.m. to 10:00 p.m., the storage system starts to discharge to reduce the “peak” of the load profile. At the same time, the energy storage system is also charged during the “non-valley” hours from 7:00 a.m. to 11:00 a.m. during daytime and during the “non-peak” hours from 10:00 p.m. to 12:00 p.m. in order to increase the consumption capacity of PV. The charging ratio always meets the energy storage operating constraints during the scheduling process. The upper-tier optimal scheduling can effectively relieve the peaking pressure of hydropower units and leave sufficient margin for improving the consumption rate of wind power and PV.

The power output of each new energy generation system is shown in Figure 13. It can be observed from Figure 13 that the highest wind power grid consumption rate is concentrated during 3:00 a.m.–4:00 a.m. and at 9:00 p.m. It is mainly due to the charging of energy storage plants in the upper optimization during this time, which provides space for the wind power to be consumed. At the same time, the wind power consumption increases at 11:00 p.m. The wind power output reaches its peak at 11:00 p.m., so the wind power consumption capacity is also increased at that time after optimized scheduling. PV grid-connected consumption is concentrated during the period from 6:00 a.m. to 5:00 p.m. during the daytime. The number of starts and stops of hydropower units was reduced to 0. It can maintain a good economy while supplying smooth power to the power system. After optimization, the consumption rate of wind power and PV reaches 66.82%.

The comparison of active network loss and grid vulnerability is shown in Figure 14. The active network loss of the system is significantly reduced compared to the pre-optimization period. The average network loss is reduced from 7.43 MW to 3.59 MW, indicating a reduction of 51.68%. The system voltage vulnerability is reduced from the initial voltage vulnerability index of 0.391 to 0.389, indicating a reduction of 0.51%.

To verify the economy and effectiveness of the optimization model proposed in this paper, three strategies are selected for comparative analysis.

Case 1

Two-tier optimization model of a power system with wind power/PV/hydropower and energy storage.

Case 2

Two-tier optimization model of a power system with wind power/PV/hydropower and energy storage with all wind power and PV consumed and energy storage considered.

Case 3

Two-tier optimization model of a power system with wind power/PV/hydropower and energy storage without access to energy storage.

All three strategies are based on the new energy generation output and typical daily load profile. The upper and lower optimization objectives of Cases 2 and 3 are the same as those of Case 1. The comparison of the calculation results of the three strategies is shown in Table 6. As can be observed from Table 6, the system load fluctuation rate of Case 1 is the smallest, at only 7.97%. The strategy gives full play to the “low storage and high discharge” function of the grid-connected energy storage power plant. It effectively reduces the peak-to-valley difference of load on the grid. It not only brings certain environmental benefits but also reduces the difficulty in system peak regulation. In Case 2, it results in the full consumption of wind and PV when the overall system load fluctuates due to the wind power anti-peak characteristics. The fluctuation rose from the initial 9.45% to 12.81%, and the load fluctuation rate increased by 35.56%. In Case 3, the fluctuation of the load remains consistent with the fluctuation of the initial load since there is no participation of energy storage plants in the upper objective function.

Table 7 provides a comparison of the economic indicators under different operation strategies, where the system operation cost of Case 1 is 34052.35 yuan, which is 22.47% lower compared to 43921.92 yuan of Case 3. The system benefit of Case 2 exceeds its operation and maintenance costs, making the system operation cost result in a negative value in numerical terms, i.e., a profitable trend. This is due to the setting of simulation parameters that can provide environmental benefits to the system.

The optimal consumption strategies and consumption rates of wind power and PV power generation are shown in Figure 15. It can be observed that it fails to reduce the load fluctuation of the system because there is no energy storage involved in the system operation in Case 3. It does not ensure the stability of the system, and there is not enough space for the consumption of wind power and PV. The consumption rate of wind power and PV can only reach 54.30%, which is 17.90% lower than the 66.14% in Case 1. It also reduces the revenue of energy storage plants in the system. Among the three strategies, it can better take into account the fluctuation of the system while maintaining a good consumption rate in Case 1, which can effectively improve the economy of the system.

Table 8 provides a comparison of the tidal current risk indicators between the three strategies. It can be observed that Case 2 voltage vulnerability index is the lowest, which is reduced by 0.51% and 3.73% compared to 0.389 and 0.402 in Cases 1 and 3, respectively. The active network loss ratio is 3.75% in Case 3, which is reduced by 13.99% and 21.88% compared to the active network loss index of 4.36 and 4.80 in Cases 1 and 2, respectively. It achieves a balance between active network loss and voltage vulnerability (Case 1) and provides a reference for the scheduling of safe and stable operation of the system.