It is challenging to forecast the power generation of wind and solar power due to the inherent intermittency, which is affected by environmental factors such as wind speed and light intensity. The literatures (Lingfors and Widén, 2016; Gholami et al., 2017) shows that the wind and solar output will tend to smooth out with the increase in the cluster scale of wind and solar stations. In this paper, since power stations in the same area usually send power via the same transmission channels, all wind stations in the same area are combined into a virtual wind power plant (VWP) and all solar power stations are combined into a virtual solar power plant (VSP). P t A W = ∑ i = 1 N W P i , t W (1) P t A S = ∑ j = 1 N S P j , t S (2) where N W is the number of wind power stations. N S is the number of solar power stations. t is the period number, t = 1,2 , . . . , T . T is the number of time periods. i is the wind power station number, i = 1,2 , . . . , N W . j is the solar power station number, j = 1,2 , . . . , N S . Moreover, P t A W is the output of the virtual wind plant in period t , MW; P t A S is the output of the virtual solar power plant in period t , MW; P i , t W is the output of wind power stations i in period t , MW; and P j , t S is the output of solar power station i in period t , MW.

The forecast output error of the wind power plant and solar power plant is converted into the output error coefficient. ε t A W = P t A W , a − P t A W , f P t A W , f (3) ε t A S = P t A S , a − P t A S , f P t A S , f (4) where ε t A W and ε t A S are the forecast output error coefficients of the wind power and solar power at time t , respectively; P t A W , a and P t A S , a are the actual output of wind power and solar power at time t , respectively; and P t A W , f and P t A S , f are the forecast output of wind power and solar power at time t , respectively.

Describing the uncertainty of wind and solar output reasonably is important for building a hybrid hydro-wind-solar complementary model. The Gaussian mixture model (GMM) is a kind of non-parametric estimation method. Compared with the parameter estimation method, it can theoretically describe any distribution. Compared with the kernel density estimation method, the GMM can avoid the bandwidth setting effect on the accuracy of the results. Therefore, in this paper, the GMM is adopted to describe the output uncertainty of wind and solar power. Given a random variable x , the expression for the Gaussian distribution is: f x 丨 μ , σ = 1 2 σ 2 π e − x − μ 2 2 σ 2 (5)

The essence of the GMM is a simple linear combination of multiple Gaussian distributions, and its expression can be described as: p x = ∑ n = 1 N ω n φ x 丨 μ n , σ n (6) where p x is the probability distribution of x . ω n is the weight coefficient, ∑ n = 1 N ω n = 1 . Using the maximum likelihood estimation method, the parameters ω n , μ n , and σ n can be determined.

With large-scale integration, the intermittent wind and solar power and the anti-peaking characteristics of wind power will challenge the peak shaving of the system. To meet the peak shaving requirements and balance the influence of wind and solar uncertainty, the flexible adjustment of hydropower output will cause frequent changes in downstream outflows and water level, which will directly affect navigation. Therefore, in this paper, peak shaving is considered the target of the hybrid hydro-wind-solar hybrid system: m i n F x = f 1 x , f 2 x T (7)

f 1 x represents the peak shaving target. Under certain conditions of incoming water, the goal of the hydro-wind-solar coordinate peak-shaving operation is to minimize the fluctuation of the residual load so that the residual load can be as smooth as possible. Therefore, the mathematical expression of the peak shaving target can be set to minimize the residual load variance: f 1 x = min ∑ t = 1 T R t − R ¯ 2 T (8) R t = C t − P t (9) R ¯ = ∑ t = 1 T R t T (10) P t = ∑ m = 1 M P m , t H + P t A W , f + P t A S , f (11) P m , t H = g H m , t , q m , t (12) where R t is the residual load in the t period, MW. C t is the original load of the grid. R ¯ is the average residual load. P t is the sum of the output of the hybrid hydro-wind-solar system at Period t . P t A W , f is the predicted output of the virtual wind plant at period t ., MW. P t A S , f is the predicted output of the virtual solar power plant at period t . P m , t H is the output of the m hydropower station at period t . m is the serial number of hydropower station, m = 1,2 , … , M . H m , t is the water head of hydropower station m at period t . q m , t is the power generation flow of hydropower station m at period t . g · is the output calculation function of the hydropower station.

f 2 x is the navigation target. The goal of navigation focuses on the elevation and variation in the downstream water level. Therefore, the minimum of the downstream water level change variance is taken as the objective function: f 2 x = min ∑ t = 1 T Z d o w n , t − Z ¯ d o w n 2 T (13) where Z d o w n , t is the downstream water level of the downstream hydropower station m at period t , m. Z ¯ d o w n is the mean value of the downstream water level of the downstream hydropower station.

In this paper, the GMM is used to analyze the distribution of output error coefficients of wind power and solar power plants. The Gaussian mixture model (GMM) is a kind of non-parametric estimation method, it can theoretically describe any distribution. Literature (Pöthkow and Hege, 2013) compared the extract uncertainty contour line effects of the non-parametric models with Gaussian, and observed that Non-parametric models have good feasibility for various data sets. Therefore, it has widely used in uncertainty description for diverse applications (Potter et al., 2009; Mihai and Westermann, 2014), and compactly model relatively complex distributions (Liu et al., 2012).With the historical forecast and actual output coefficient data of the solar and wind power plant, the deviation value between the actual and the forecasted output coefficient can be calculated. The GMM is used to fit the deviation value data of the output coefficient, and the probability density function (PDF) p x of the output error coefficient x is obtained as shown in Figure 1.

From the PDF graph of the output errors of the solar power plant and the wind power plant, the output coefficient deviation value can be obtained given the confidence interval [0.05, 0.95]. ρ t A W , min , ρ t A W , max , ρ t A S , min , and ρ t A S , max correspond to the endpoints of the confidence interval [0.05, 0.95] at each moment, and the upper/lower limit of the wind and solar output at each moment can be calculated: P t A W , min = 1 + ρ t A W , min * P t A W , f (31) P t A W , max = 1 + ρ t A W , max * P t A W , f (32) P t A S , min = 1 + ρ t A S , min * P t A S , f (33) P t A S , max = 1 + ρ t A S , max * P t A S , f (34) where P t A W , min , P t A W , max , P t A S , min , and P t A S , max are the lower limit of the wind power output range, the upper limit of the wind power output range, the lower limit of the solar output range and the upper limit of the solar output range during period t , respectively.

Confidence in chance constraints (14) and (15) represents the probability that the power gird can consume all the wind and solar power, and whether the constraints are established is affected by the wind and solar forecast error. The chance constraint (14) means that when the hydropower is at its maximum value, the probability that the actual output of the hydro-wind-solar system is greater than the forecast output should be greater than the confidence. If the actual output of the wind and solar power is too small, the constraint may be destroyed. Therefore, the minimum output P t A W , min of wind power and the minimum output P t A S , min of solar power within the confidence level ϑ 1 can be found through the CDF of wind and solar power forecast error. The opportunity constraint (14) can be transformed into that when the sum of wind power and solar power output is P t A W , min + P t A S , min , the output of the hydro-wind-solar system meets the requirements: ∑ m = 1 M P m H , max − σ × C t + P t A W , min + P t A S , min ≥ P t (35)

The chance constraint (15) means that when the hydropower is at the minimum value, the probability that the actual output of the hydro-wind-solar system is less than the forecast output should be greater than the confidence. If the actual output of scenery is too large, the constraint may be destroyed. Similarly, the opportunity constraint (15) can be transformed into that when the sum of wind power and photovoltaic output is P t A W , max + P t A S , max , the output of the hydro-wind-solar system meets the requirements: ∑ m = 1 M P m H , min + σ × C t + P t A W , max + P t A S , max ≤ P t (36)

The ε − c o n s t r a i n t method simplifies the problem by retaining only one primary objective function and transforming the remaining objective functions into constraints (Mavrotas, 2009; De Santis et al., 2022); moreover, this method has been widely used in the solution of multi-objective models (Biswas et al., 2018; Fathipour and Saidi-Mehrabad, 2018). Compared with the weighted method, the ε − c o n s t r a i n t method has a smaller search range and better performance in efficiently finding the Pareto front solution. Regarding the dual-objective optimization problem of finding the minimum value, its mathematical description is as follows: min ⁡ ⁡ f 1 x , f 2 x s . t .   A x ≤ b (37) where A is a matrix of appropriate dimensions; and b is a column vector of appropriate dimensions.

For the dual objective, if Objective f 1 x can be regarded as the primary objective function, then Objective f 2 x is transformed into the constraints of f 1 x . The specific operation steps are as follows:

Take f 1 x as the objective function, find its minimum value f 1 min under the original constraints, and find the value of f 2 x at this time, which is regarded as the feasible maximum value f 2 min of f 2 x . The solution model required for this step is: min f 1 x s . t .   A x ≤ b (38)

Take f 2 x as the objective function, and find the minimum value f 2 min of f 2 x under the original constraints. The range of f 2 x is f 2 min , f 2 max . The solution model required for this step is: min f 2 x s . t .   A x ≤ b (39)

Take K points ε k k ∈ 1,2 , … , K in f 2 min , f 2 max , and generate K new constraints: f 2 x ≤ ε k .

Take f 1 x as the objective function and add the new constraints to the constraints. After the calculation is performed, K result arrays f 11 , f 21 、 f 12 , f 22 … f 1 k , f 2 k can be obtained, that is, the approximate Pareto maximum excellent Frontier. The solution model required for this step is: min f 1 x s . t .   A x ≤ b f 2 x ≤ ε k (40a)

Yunnan Province is located in southwest China. In this region, with sufficient rainfall occurs, and several large rivers travels through, can be found. Yunnan province is endowed with abundant hydropower resources. In late 2022, the hydropower installed capacity in this region surged to 78.02 GW. This value is close to the hydropower installed capacity of Canada, which is the world’s third largest hydropower capacity. In addition, Yunnan Province is rich in renewable energy. In late 2022, wind power and solar power has surged to 8.83 and 4.17 GW respectively. However, there is very large potential for expanding wind and solar power resources, since the developed capacity accounts for only 5% of the developable capacity. According to the “New Energy Construction Plan of Yunnan Province in 2022” issued by Yunnan Province, Yunnan will accelerate the development of wind and solar power in the next few years. Since wind and solar power generation is depend on the weather, Yunnan is facing challenges of broad-scale renewable energy consumption and power grid peak shaving brought by the broad-scale integration of VREs. Therefore, using the flexibility of hydropower to achieve the multi-energy complementation of hydro-wind-solar power is one of the best ways for Yunnan to achieve large-scale consumption of VREs.

As shown in Figure 2, the Jinghong hydropower station and the Ganlanba hydropower station are both located in Jinghong City, Xishuangbanna Prefecture, Yunnan Province. They are the last two hydropower stations in the Lancang/Mekong River before flow out of China. Due to Jinghong’s rated flow and the restriction of downstream channel conditions, Jinghong’s adjustment of power generation flow will cause large fluctuations in downstream water level and flow, which will adversely affect navigation. To solve this problem, the main development task of the yet-to-be-constructed Ganlanba hydropower station is to provide reverse regulation for the Jinghong hydropower station to meet the downstream navigation requirements. Therefore, the Jinghong-Ganlanba cascade hydropwer station must undertake complex comprehensive utilization tasks. Considering Yunnan’s ever-increasing renewable energy construction plan, this cascade system will also undertake the severe task of VRES complementarity. Therefore, it is required to find a reasonable operation mode considering the conflicting objectives of generation and navigation. In this paper, the Jinghong-Ganlanba cascade hydropower stations are taken as the background. The load curve, interval runoff, initial and final water levels, and other conditions involved in the calculation process refer to the actual operating data of the power grid and power station, as shown in Supplementary Table S1.

In this paper, a non-linear model is built, where a quadratic polynomial fitting function is used to express the water level-storage function and tail water level-discharge function. The model was solved using Lingo18.0.

According to Section 2.1, all wind power plants in the same region are considered as one wind power plant with the installed capacity of 208.5 MW and a solar power plant with the installed capacity of 300 MW. Figure 3 illustrates the typical wind and solar power daily output during the dry season (January, February, and March) of the lower reaches of the Lancang River, which were chosen to evaluate the scheduling situation of cascade hydropower stations under difficult conditions.

Since Ganlanba has not been constructed, three cases are considered to analyze the impact of the construction of the Ganlanba on the overall results and the impact of the navigation requirements on the separate operation of the Jinghong hydropower station: 1) Case 1: only the Jinghong hydropower plant is considered in the hybrid hydro-wind-solar complement operation model, but does not take navigation constraints into account; 2) Case 2: only the Jinghong hydropower plant is considered and takes navigation constraints into account and 3) Case 3: both Jinghong and Ganlanba are considered, and navigational constraints are also taken into account.

Figure 4 and Supplementary Table S2 show the distribution and range of the Pareto front for the three cases which reflect the conflict between the two objective functions. The rectangular area in Figure 4 marks solutions at the turning part of the Pareto front, which are hydropower output processes that can better meet the requirements of peak shaving and navigation at the same time. Under the experimental conditions, the Pareto front turning of the cascade solution is more obvious, indicating that it can moderate the conflict more effectively. Because the navigation constraints limit the output of Jinghong, the range of feasible solutions is reduced, and even some solutions with good peak-shaving effects turn into infeasible solutions.

Figure 5 illustrate the differences of solutions on the Pareto front for Case1 and show the changes in the output and water abandonment under Case1, respectively, by selecting different downstream water level variance limits ( f 2 = 1.4, 1.0, 0.6, 0.2, 0.01, 0.0001). When the target requirements for f 2 tighten, Jinghong’s discharge adjustment ability is restricted, which result in Jinghong abandoning water for peak shaving. For Case1, the peak shaving and navigation task requirements cannot be achieved at the same time.

In the Pareto front of Case 1, three scenarios of f 2 = 1 , f 2 = 0.5 , and f 2 = 0.01 are selected as the solution that prefers a better peak shaving effect, the neutral case, and the solution that prefers a better navigation effect, respectively. These three points can effectively compare the gaps in peak regulation and navigation between these cases, as shown in Supplementary Table S3.

Supplementary Table S3 summarizes the eigenvalues of each case. It is obvious that under the same f 2 , Case 3’s peak shaving impact is superior based on the peak-to-valley difference of the residual load. Since the target f 1 tends to reduce the hydropower output at the load valley, the hydroelectric plant can only supply reserve capacity at the load valley leaving a constant minimum load at 6178.285 MW. Figure 6 show the output and water abandonment of Case 3 under three f 2 conditions. From Figure 5, 6, it can be concluded that: 1) in Figure 6, the output of Ganlanba slightly changes under different f 2 conditions, which indicates that Ganlanba has a weak adjustment ability and is mainly responsible for navigation tasks. Therefore, Figure 5, 6 show that with the addition of Ganlanba, Jinghong’s adjustment capacity can be unleashed, and a better result can be obtained. 2) Figure 5, 6 also show that Ganlanba’s participation can significantly reduce the abandoned water for peak shaving, which may successfully address conflicting objectives of peak shaving and navigation.

Figure 7 shows the downstream water levels before and after the addition of Ganlanba and displays the affected navigation due to Ganlanba’s involvement. When the downstream water level variance is set to 0.01, the corresponding solution of Case 3 is located at the turning point of the Pareto front. In Figure 7, the Pareto front of Case 2 also has a corresponding solution, which can be compared. Without Ganlanba, the navigable water level is maintained at 537 m for 1/3 of the time, which indicates that the downstream channel barely meets the navigation requirements. For Case 3, the Ganlanba’s addition lowers the expense of maintaining navigation, and it maintains the navigation target’s ideal value between 529.45 m and 529.78 m. The safety of downstream navigation has been further assured compared to the minimum water level requirement of 525 m downstream of the Ganlanba.

In conclusion, the Ganlanba hydropower plant can successfully reduce the tension between peak shaving and navigation, address the issue of water abandonment for peak shaving, and raise the level of safety for downstream navigation.

Literature (Shen et al., 2021) proposes that the load fluctuation coefficient can be used to evaluate the peak shaving effect, and the revised index mainly reflects the overall smoothness of the load curve, which is better than the visualization effect of the residual load mean square error. Therefore, in this paper, this index is used to reflect the peak shaving effect in different scenarios: α = 1 T ∑ t = 1 T R r , t − R r ¯ 2 R r ¯ (40b) where, R r , t is the residual load of the scenario at the period t ; and R r ¯ is the average residual load of the scenario r , with r = 1,2 , … , 1000 .

In Case 3, the actual operation of the cascade is simulated under two conditions: one condition with complementary constraints for the solar and wind output (Eqs 33, 34) and the other condition without such constraints. When there is no hydro-wind-solar complement constraint, the residual load in each scenario is: R r , t = C t − ∑ m = 1 M P m , t H − P r , t (41) P r , t = P r , t A W + P r , t A S (42) where, R r , t is the residual load of the r th scenario at period t without hydro-wind-solar complement constraints; P m , t H is the output of hydropower station m at period t ; P r , t is the wind and solar power output of the scenario r at period t ; P r , t A W , P r , t A S represent the wind and solar power output of scenario r at period t respectively.

When considering the hydro-wind-solar complementary constraint, here are the residual loads for each scenario: R r , t = R t , P r , t < P r , t max   a n d   P r , t > P r , t min R t + P r , r − P r , t max , P r , t > P r , t max R t + P r , t − P r , t min , P r , t < P r , t min (43) P t max = P t A W , max + P t A S , max (44) P t min = P t A W , min + P t A S , min (45) where P t max is the upper limit of the range for the sum output of wind and solar; and P t min is the lower limit of the range for the sum output of wind and solar.

Through Eqs 40–45, the residual load under 1,000 simulation scenarios was carried out, and the distribution of the load fluctuation coefficients is shown in Figure 8. It is obvious that considering hydro-wind-solar complementarity can effectively control the residual load fluctuation in real operation where actual wind and solar power deviate from the forecast value. In summary, hydro-wind-solar complement operation can effectively reduce the adverse impact of wind and solar uncertainty on the peak shaving results.

Under the condition of fixed hydropower output and without providing wind and solar reserve capacity, the navigation effect will not be affected by VRE when actual wind and solar output deviate from the forecast value. In the case of hydro-wind-solar complement operation, the part of the output deviation between the upper and lower limits of the reserve capacity needs to be absorbed by changing the hydropower output. Due to the small adjustment capacity of Ganlanba, in this situation, Jinghong is required to undertake the adjustment tasks. Since the navigation is determined by the discharge of Ganlanba, the picture A of Figure 9 shows the maximum and minimum water levels of the Ganlanba reservoir and the upper limit and lower limit water levels of the Ganlanba reservoir under 1,000 scenarios. In all scenarios, Ganlanba reservoir’s water level will not exceed its limit, which means that when wind and solar forecast errors cause the upper reservoir outflow to change, Ganlanba has sufficient adjustment capacity to address it. The picture B of Figure 9 shows the downstream water level variance of Ganlanba Reservoir under 1,000 scenarios. The mean square deviations of the downstream water level in most scenarios are smaller than the predicted situation, and a few are slightly larger than the predicted situation. However, the difference is small and negligible. This demonstrates that Ganlanba can maintain navigational conditions despite VREs’ forecast errors.

That is, wind and solar uncertainties have little effect on navigation.

In this section, different wind and solar power installed capacity scales are selected to explore their impact on the simulation results. There shows four wind and solar power scales. Scenario 1 elects solar power scales of 300 MW and wind power scales of 208.5 MW, which are the wind and solar installed capacity in 2020. Referring to the “List of New Energy Projects in the 14th Five-Year Plan of Yunnan Province,” Xishuangbanna will add 1.43 GW of solar power during the “14th Five-Year Plan” period. However, no wind power projects are planned. Scenario 2 sets the estimated installed capacity of solar and wind power in 2025 as 1730 and 208.5 MW. Scenario 3 sets the projected installed solar and wind capacity in 2030 as 3160 and 208.5 MW, and the growth rate is consistent with that in Scenario 2. Scenario 4 doubles the wind-power capacity installed based on Scenario 3, and sets the estimated installed capacity of solar and wind power as 3160 and 417 MW.

Figure 10 shows the residual load curves of each scenario, and Figure 11 shows the hydropower output gap between Scenario 1 and Scenario 2, 3 and 4, in which a negative value indicates that hydropower output will decrease as solar/wind power output increases, while a positive value indicates the opposite. As shown in these two figures, for Scenario 2 and 3, when solar power output increases, the residual load reduce from 8:30 to 20:00, and hydropower will reduce output during the load saddle periods and uses the saved water to increase the output during the load peak periods. This is because solar power is mainly generated during the daytime which is the load saddle, and hydropower adapts to achieve hydro-solar complementarity and improve the peak shaving effect. The simulation results of Scenarios 2 and 3 in Supplementary Table S4 also prove that increasing the output of solar power has a positive effect on the peak shaving of the model.

As shown in Table 4, the peak-to-valley difference of the residual load increases from Scenario 3 to Scenario 4. This is because hydropower has reached the best peak shaving performance in Scenario 3. When the wind power capacity is increased from Scenario 3 to Scenario 4, the hydropower output changes little with the change of reserve capacity, as shown in Figure 11, and helps little in further peak shaving. Therefore, with the increase in wind power capacity, although both the residual peak load and valley load decrease, as the anti-peak characteristic of wind power shows in Figure 3., the peak-valley difference will increase.

Figure 11 shows that the hydropower output changes little, which means that the outflow of hydropower changes little. From Scenario 1 to Scenario 4, the increase in wind power and solar power scale will not have a large negative impact on navigation.

In summary, during the dry season, the increase in the scale of solar power output has a positive impact on the peak shaving of the hydro-wind-solar complementarity system; in contrast, wind power tends to have a negative impact on peak shaving.