Study on the Optimal Operation of a Hydropower Plant Group Based on the Stochastic Dynamic Programming with Consideration for Runoff Uncertainty

Abstract

Hydropower plant operation reorganizes the temporal and spatial distribution of water resources to promote the comprehensive utilization of water resources in the basin. However, a lot of uncertainties were brought to light concerning cascade hydropower plant operation with the introduction of the stochastic process of incoming runoff. Therefore, it is of guiding significance for the practice operation to investigate the stochastic operation of cascade hydropower plants while considering runoff uncertainty. The runoff simulation model was constructed by taking the cascade hydropower plants in the lower reaches of the Lancang River as the research object, and combining their data with the copula joint function and Gibbs method, and a Markov chain was adopted to construct the transfer matrix of runoff between adjacent months. With consideration for the uncertainty of inflow runoff, the stochastic optimal operation model of cascade hydropower plants was constructed and solved by the SDP algorithm. The results showed that 71.12% of the simulated monthly inflow of 5000 groups in the Nuozhadu hydropower plant drop into the reasonable range. Due to the insufficiency of measured runoff, there were too many 0 values in the derived transfer probability, but after the simulated runoff series were introduced, the results significantly improved. Taking the transfer probability matrix of simulated runoff as the input of the stochastic optimal operation model of the cascade hydropower plants, the operation diagram containing the future-period incoming water information was obtained, which could directly provide a reference for the optimal operation of the Nuozhadu hydropower plant. In addition, taking the incoming runoff process in a normal year as the standard, the annual mean power generation based on stochastic dynamic programming was similar to that based on dynamic programming (respectively 305.97 × 10⁸ $\mathrm{k}\mathrm{W}\cdot \mathrm{h}$ and 306.91 × 10⁸ $\mathrm{k}\mathrm{W}\cdot \mathrm{h}$), which proved that the operation diagram constructed in this study was reasonable.

1. Introduction

The optimal operation of hydropower plant groups is an issue for water conservancy and hydropower systems based on the theory of the optimal single hydropower plant operation method and comprehensively considering the water conservancy relationship between hydropower plants [1,2]. The optimal operation is to reorganize the temporal and spatial distribution of water resources according to the comprehensive utilization tasks assigned to the hydropower plant and based on the allocation capacity of the reservoir, achieving the target of promoting benefits and eliminating hazards [3,4]. The research on the optimal operation of hydropower plant groups has a long history. Ponnambalam et al. [5] showed that an implementation of the Karmarkar’s interior-point LP algorithm with a newly developed stopping criterion could solve the optimization problems of large-scale multi-reservoir operations more efficiently than the simplex method. To address the deficiency of the genetic algorithm in the application of reservoir group operation, Oliveira and Loucks [6] used an improved genetic optimization algorithm to effectively formulate the operation rules of reservoir groups. Chandramouli and Raman [7] established and solved the optimal operation model of reservoir groups by using a dynamic programming algorithm and neural network. In order to enhance the performance of evolutionary algorithms, Arunkumar and Jothiprakash [8] introduced chaos technology to generate initial population and applied it to the operation optimization of their reservoir group system.

The optimal operation of a hydropower plant means meeting the requirements of flood control and benefits on the premise of ensuring the safety of reservoir engineering facilities and the ecological environment. [9,10] Generally, under the optimization objective, the hydropower plant operation adopts the historical runoff as the input and obtains the optimal operation decisions for the hydropower plant for each period through intelligent algorithm technology [11,12]. However, there is a certain gap between the determined hydropower plant operation rules and the actual dispatching results [13]. One reason for this is that the stochastic inflow is not considered. The inflow is uncertain under the influence of both anthropic and natural factors [14,15], which, in time, has a certain correlation between different periods as well. Dynamic programming (DP) algorithm [16,17] is a process of searching for the optimal decision step by step, in order, at multiple stages, but it is difficult to solve the stochastic multi-stage problems. Therefore, the stochastic dynamic programming (SDP) algorithm [18,19] can be used to effectively work out operation rules for the runoff uncertainty by combining the probability theory and dynamic programming principle.

The runoff process during the process of hydropower plant operation is uncertain, which limits the traditional algorithms and determined dynamic programming, while the SDP algorithm fully takes the uncertain factors into consideration [20,21]. Little [22] applied the Markov theory to construct the SDP model of a hydropower plant. Howard [23] proposed a method coupling the Markov theory with the DP (MDP), which avoids prematurely falling into the local solutions when compared with the traditional model for the multi-year issues. Rossment [24] combined the Lagrange multiplier with DP theory to establish a chance-constrained programming model. Bras et al. [25] applied the SDP method to the Aswan Dam, and the results showed that the SDP model is effective when considering runoff prediction information. Saadat M and Asghari K [26] made an improvement of the SDP algorithm by properly adjusting the interval indices of the reservoir storage capacity, which resulted in the values of the objective function increasing by 30%. Therefore, it is critical to explore the issue of operation discission-making for cascade hydropower plants under the uncertainty of the incoming runoff.

The development of large-number cascade hydropower plants in the whole basin plays a great role in promoting the economic development of the Southwest of China. In previous research, most of them were certain operation, and the input used was the observed runoff [3]. However, in this paper, the distribution of uncertain incoming runoff was constructed using the copula joint function and Gibbs sampling based on the large number of simulations of observed runoff data. Moreover, Lancang River, as a transboundary river, is located in Southwest China, and the water resources are sufficient; while on the other hand, the incoming runoff in this area is uncertain as well. Therefore, it is of great significance to investigate the region’s hydropower plant group while taking the uncertainty of incoming runoff into consideration.

2. Study Area and Data

The source of the Lancang River is located in the Qinghai Province (Figure 1). Flowing out of China from Mengla county of Yunnan Province, it goes through five Southeast Asian countries: Myanmar, Laos, Vietnam, Cambodia, and Thailand. The Lancang River has a total length of about 2139 km in China, an area of 164,800 km² and an annual outbound water volume of 765 × 10⁸ m³ [27]. There are seven cascade hydropower plants on the Lancang River, including the Gongguoqiao (GGQ), Xiaowan (XW), Manwan (MW), Dachaoshan (DCS), Nuozhadu (NZD), Jinghong (JH), and Ganlanba (GLB) power plants. The lower reaches of the Lancang River include three hydropower plants: namely, NZD, JH, and GLB. The NZD hydropower plant has a multi-year regulation capacity, and its comprehensive utilization tasks are mainly for power generation, irrigation, flood control, shipping, ecology, tourism, etc. The JH hydropower plant is located 5 km away from the northern suburb of Jinghong city, and its main task is power generation, taking shipping, flood control, tourism, etc. into account. GLB is the reverse regulation power plant of the JH hydropower plant, and its main task is power generation, taking the needs of shipping and ecological water into account [28].

In this study, the natural monthly runoff of the NZD and JH hydropower plants from 1953 to 2013 and the natural daily runoff (1953–2013) collected from the Yunjinghong hydrological station and the Huaneng Lancang River Hydropower Co., Ltd. (Kunming, China) were used to investigate and derive the hydropower plant operation. And the location of hydrological stations and hydropower plants were shown in Figure 2.

3. Methods

In this section, we run through the details of the proposed methodology, including the copula joint function, Gibbs simulation model, and the constructed stochastic operation model. Figure 3 demonstrates the framework of the main research content in this paper.

3.1. Construction of Joint Distribution Model for Monthly Runoff

3.1.1. Marginal Distribution Model of Incoming Runoff

In this paper, generalized extreme value distribution (Gev), logarithmic normal distribution (Logn), gamma distribution (Gamma), normal distribution (Norm), and weibull distribution (Weibull) were selected as the candidate marginal distributions of measured incoming runoff sequence [29,30,31]. The above probability distributions were chosen because they have the advantages of strong adaptability and a good fit, and they are widely used in the construction of marginal distribution models for frequency analysis of hydrological events. The probability density functions are listed here as follows:

Probability density function for Gev distribution:

$$f(x/\mu ,\sigma ,\xi )={\sigma }^{-1}{[1+\xi (\frac{x-\mu }{\sigma })]}^{-\frac{1}{\xi }-1}\exp \left\{{[1+\xi (\frac{x-\mu }{\sigma })]}^{-\frac{1}{\xi }}\right\}$$

(1)

Probability density function for Logn distribution:

$$f(x/\mu ,\sigma )=\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{{(\ln x-\mu )}^2}{2{\sigma }^2})$$

(2)

Probability density function for Gamma distribution:

$$f(x/\alpha ,\beta )=\frac{1}{{\beta }^{\alpha }\mathrm{\Gamma }(\alpha )}{x}^{\alpha -1}\exp (\frac{-x}{\beta })$$

(3)

Probability density function for Norm distribution:

$$f(x/\mu ,\sigma )=\frac{1}{\sigma \sqrt{2\pi }}\exp (-\frac{{(x-\mu )}^2}{2{\sigma }^2})$$

(4)

Probability density function for Weibull distribution:

$$f(x/\lambda ,k)=\frac{k}{\lambda }{(\frac{x}{\lambda })}^{k-1}{e}^{-{(x/\lambda )}^k},x\ge 0$$

(5)

where μ is average, σ is standard deviation, ξ, α, and k represent shape parameters, β and λ represent scale parameters, and x is the independent variable.

3.1.2. Joint Distribution Model of Adjacent Monthly Incoming Runoff

The introduction of the copula function caused a technological innovation in the field of hydrology, which could solve the multivariate frequency problem [32,33,34].

Supposing X, Y represent the adjacent monthly incoming runoff, the corresponding designed values are x, y, and the marginal distributions are F_X (x), F_Y (y). According to Sklar’s theorem, the joint distribution of X, Y can be described by the 2-D copula function C:

$$F(x,y)=C_{\theta }(F_X(x),F_Y(y))$$

(6)

where, F(x,y) is the joint distribution of X, Y; θ is the parameter of the copula function.

There are a great many types of copula functions, among which the Archimedean copula is widely used in the field of hydrology. In this chapter, 3 commonly used Archimedean copula functions (Gumbel copula, Clayton copula, and Frank copula) were selected to construct the joint distribution model of the monthly runoff in the NZD reservoir [35,36,37], the function is described as follows:

Gumbel copula:

$$C(u,v)=\exp \left\{-{[{(-\ln u)}^{\theta }+{(-\mathrm{lnv})}^{\theta }]}^{1/\theta }\right\}$$

(7)

Clayton copula:

$$C(u,v)={(u^{-\theta }+v^{-\theta }-1)}^{-1/\theta }$$

(8)

Frank copula:

$$C(u,v)=-\frac{1}{\theta }\ln [1+\frac{e^{-\theta u}-1)(e^{-\theta v}-1)}{(e^{-\theta }-1)}]$$

(9)

where u and v represent marginal distribution function, namely $u=F_X(x),v=F_Y(y)$.

3.2. Stochastic Simulation Model of Runoff

Gibbs sampling is a kind of Markov method, and the sampling process is as follows: each variable is sampled from the conditional distribution of other variables in a fixed order, a Markov chain converging to the target probability distribution is constructed, and samples are extracted from the chain [38,39].

Calculating conditional distribution by copula function:

$$F_i(x_{i+1}/x_i)=\partial C_{i,i+1}(x_i,x_{i+1})/\partial x_i,i=1,2\cdots ,12$$

(10)

The Gibbs sampling method is used for random simulation of the monthly runoff, and the steps are as follows:

(1) Generate 2 random numbers, a₀, $a_1\in (0,1)$.

(2) $P\left(X_{1,1}\le x_{1,1}/X_{2,1}=a_0\right)=\partial C\left(x_{1,1},x_{2,1}\right)/\partial x_{2,1}=a_1$, $X_{i,j}$ represents the monthly runoff in the i month of the j year, according to the formula of conditional distribution, derive v, namely $x_{1,1}$.

(3) Generate 12 random numbers $a_2, a_3,\cdots ,a_{13} \in (0,1)$ and solve the following equations:

$$\{\begin{array}{l}P(X_{2,1}\le x_{2,1}/X_{1,1}=x_{1,1})=\partial C(x_{1,1},x_{2,1})/\partial x_{1,1}=a_2,\\ P(X_{3,1}\le x_{3,1}/X_{2,1}=x_{2,1})=\partial C(x_{2,1},x_{3,1})/\partial x_{2,1}=a_3,\\ \hfill ⋮\hfill \\ P(X_{12,1}\le x_{12,1}/X_{11,1}=x_{11,1})=\partial C(x_{11,1},x_{12,1})/\partial x_{11,1}=a_{12},\\ P(X_{1,2}\le x_{1,2}/X_{12,1}=x_{12,1})=\partial C(x{1}_{2,1},x_{1,2})/\partial x_{12,1}=a_{13},\end{array}$$

(11)

The monthly runoff in all 12 months of the first year and January of the second year are obtained ($x_{1,1},x_{2,1},\cdots ,x_{12,1},x_{1,2}$).

(4) Repeat step (3) for n times to obtain n-year simulated monthly runoff.

3.3. Construction of Stochastic Operation Model and Solution

3.3.1. Objective Function

Using 5000 groups of stochastic simulated runoff of the NZD hydropower plant, the stochastic optimal operation model of cascade hydropower plants in the lower reaches of the Lancang River was constructed with the objective of maximizing the power generation.

$$F=max\left\{{\varphi }_t+E_{t+1,n}\right\}$$

(12)

where F is the generation expectation from time t to time n; ${\varphi }_t$ is the power generation at time t; $E_{t+1,n}$ is the generation expectation from t + 1 to n, namely, the generation expectation in the remaining period.

3.3.2. Constraint Condition

The main constraints of these models are described as follows [3,26]:

• Water balance constraint:

  $$V_t=V_{t-1}+(q_t-Q_t)\cdot \Delta t-S_t$$

  (13)
• Discharge constraint:

  $$Q_{t,min}\le Q_t\le Q_{t,max}$$

  (14)
• Water level constraint:

  $$Z_{t,min}\le Z_t\le Z_{t,max}$$

  (15)
• Output constraint:

  $$N_{min}\le N_t\le N_{max}$$

  (16)
• The minimum ecological flow constraints:

  $$Q_t\le bas{e}_{\left(i,t\right)}$$

  (17)
• The minimum shipping flow constraints:

$$Q_t\le Q_t{}^e$$

(18)

where $V_{t-1}$ is the initial storage capacity, and $V_t$ is the final storage capacity; $q_t$ represents incoming runoff, $Q_t$ represents the discharge, and $S_t$ represents waste water; $Q_{t\min }$ and $Q_{t\max }$ are the minimum and the maximum discharge; $Z_t$, $Z_{t,\min }$, and $Z_{t,\max }$ are, respectively, the hydropower plant’s water level at month t, dead water level, and the highest water level of hydropower plant i; $N_t$, $N_{\min }$, and $N_{\max }$ are, respectively, the hydropower plant output at month t, firm output, and the installed capacity; $bas{e}_{\left(i,t\right)}$ is the ecological base flow at month t of hydropower plant i; and $Q_t$ is the minimum flow for shipping.

3.3.3. Solution Algorithm

The SDP algorithm [40,41,42] was used to solve the stochastic operation model. When the SDP is applied to reservoir optimal operation, the runoff process is described by Markov chain, combined with the runoff forecast in the current period, and the incoming runoff in the future period is expressed by transfer probability. Therefore, the power generation in the current period can only be calculated according to the forecast runoff in the current period, and the power generation in the future period is represented by the power generation expectation.

The operation period (monthly) was taken as the stage variable, the water level of the NZD reservoir was divided with the accuracy of 1 m, the water level $Z_{b.t}$ at the beginning of the month was taken as the state variable, and the water level $Z_{e,t}$ at the end of the month was taken as the decision variable. The water balance equation was used as the state transition equation.

Setting the runoff in period t as $q_{t,1},q_{t,2},\cdots q_{t,n}$ (n = 10), and the transfer probability matrix from period t to period t + 1 is as follows:

$$P_{ijt}=\left[\begin{array}{l}{P}_{11}{}_t P_{12}{}_t \cdots P_{1n}{}_t\\ P_{21}{}_t P_{22}{}_t \cdots P_{2n}{}_t\\ ⋮ ⋮ ⋮\\ P_{n1}{}_t P_{n2}{}_t \cdots P_{nn}{}_t\end{array}\right]$$

(19)

Assuming that the incoming runoff $q_{t,i}$ at k period has been given, the formula of optimal power generation expectation is as follows:

$$f_t^{*}(Z_b{}_{,t},q_{t,i})=max\left\{{\mathrm{\Phi }}_t(Z_{b,t},Z_e{}_{,t},q_{t,i})+{\displaystyle \sum _{j=1}^n{p}_{ijt}{f}_{t+1}^{*}(Z_{e,t},q_{t+1,j})}\right\} (t=1,2,\cdots n)$$

(20)

where $f_t^{*}(Z_b{}_{,t},q_{t,i})$ represents the generation expectation at time t; ${\mathrm{\Phi }}_t(Z_{b,t},Z_e{}_{,t},q_{t,i})$ represents the power generation at time t; $Z_b{}_{,t}$ is the water level at the beginning of the month; $Z_e{}_{,t}$ is the water level at the end of the month; $q_{t,i}$ is the amount of runoff at time t; $i\in (1,10)$, $p_{ijt}$ represents the transfer probability matrix from period t to period t + 1;$j\in (1,10)$,$f_{t+1}^{*}(Z_{e,t},q_{t+1,j})$ represents the generation expectation at time t + 1; and $q_{t+1,j}$ is the amount of runoff at time t + 1.

Regardless of the initial state and initial decision, for the state formed by the previous decision, the other decisions must constitute the optimal strategy: that is, the n-dimensional optimization problem must be transferred into n 1-dimensional problems. Compared with the DP algorithm, the difference of the SDP is that the incoming runoff is a stochastic process.

4. Results and Discussion

4.1. Stochastic Simulation of Runoff

4.1.1. Marginal Distribution for Incoming Runoff

It was concluded from the cumulative frequency (in Figure 4) that the 5 selected distribution curves (Weibull, Log-normal, Gamma, Normal and Gev) fit well with the curve of the experience frequency. The K-S test, AIC, and RSME were adopted to quantitatively evaluate the fitting effect of the marginal distribution, and the results are shown in

Table 1. K-S, AIC, and RMSE of incoming runoff in the NZD reservoir.

Index	Distribution	January	February	March	April	May	June
K-S test	Weibull	0.240	0.749	0.314	0.854	0.539	0.420
	Logn	0.449	0.895	0.369	0.774	0.954	0.975
	Gamma	0.531	0.958	0.315	0.888	0.970	0.954
	Norm	0	0	0	0	0	0
	Gev	0.602	0.999	0.235	0.994	0.863	0.598
AIC	Weibull	296.555	257.53	275.016	303.779	397.945	493.377
	Logn	286.587	255.956	262.653	300.312	392.247	483.081
	Gamma	286.931	255.069	263.75	299.536	391.855	484.338
	Norm	288.843	254.317	266.92	299.446	393.985	490.282
	Gev	288.81	256.575	263.048	301.98	394.265	485.293
RMSE	Weibull	0.062	0.035	0.054	0.029	0.038	0.044
	Logn	0.04	0.029	0.038	0.037	0.025	0.022
	Gamma	0.041	0.025	0.039	0.032	0.021	0.025
	Norm	0.046	0.02	0.045	0.024	0.025	0.038
	Gev	0.306	0.316	0.302	0.302	0.275	0.262
Index	Distribution	July	August	September	October	November	December
K-S test	Weibull	0.364	0.229	0.471	0.373	0.15	0.104
	Logn	0.535	0.915	0.697	0.977	0.85	0.67
	Gamma	0.431	0.722	0.624	0.907	0.606	0.498
	Norm	0	0	0	0	0	0
	Gev	0.291	0.293	0.547	0.517	0.216	0.23
AIC	Weibull	531.949	588.492	560.242	516.981	472.988	371.753
	Logn	523.338	578.995	547.004	505.39	451.806	351.651
	Gamma	524.208	580.171	548.8	506.957	455.697	354.275
	Norm	527.805	586.518	555.915	513.035	468.245	361.636
	Gev	525.192	581.272	548.725	507.028	451.248	350.342
RMSE	Weibull	0.052	0.05	0.051	0.056	0.074	0.075
	Logn	0.038	0.021	0.027	0.024	0.031	0.033
	Gamma	0.041	0.028	0.033	0.032	0.041	0.04
	Norm	0.049	0.047	0.047	0.049	0.064	0.056
	Gev	0.281	0.255	0.266	0.269	0.265	0.29

(the optimal marginal distribution is bolded). It is clearly found from Table 1 that for the K-S test, there were 8 months where the corresponding optimal distribution was Logn distribution (the optimal distribution of January, February, and April was GEV distribution, and that of May was Gamma distribution); similarly, for the AIC and EMSE there were 7 months and 9 months, respectively, where the optimal distribution was Logn distribution. Overall, the Logn was the best distribution of the five distributions. Therefore, the Logn distribution was selected as the monthly marginal distribution of the incoming runoff sequence.

4.1.2. Copula Joint Distribution for Incoming Runoff

The copula function is used to derive the joint function between adjacent monthly runoff. The parameters of three candidate copula functions were estimated by the maximum likelihood method, and the copula function was selected using AIC criterion. The results are shown in

Table 2. Parameters and AIC of joint distribution of inflow in the NZD reservoir.

Monthly	Parameters	Copula
		Clayton	Frank	Gumbel
January–February	θ	2.84	12.11	2.99
	AIC	−74.30	−94.26	−85.33
February–March	θ	1.04	4.95	1.68
	AIC	−24.75	−29.79	−25.47
March–April	θ	0.94	4.93	1.78
	AIC	−19.82	−29.76	−32.66
April–May	θ	0.72	3.54	1.45
	AIC	−15.08	−15.07	−11.72
May–June	θ	0.60	2.07	1.29
	AIC	−6.74	−4.90	−7.35
June–July	θ	0.74	2.92	1.40
	AIC	−10.76	−12.24	−12.72
July–August	θ	0.78	2.61	1.37
	AIC	−7.86	−9.88	−11.78
August–September	θ	1.09	3.60	1.50
	AIC	−20.51	−16.76	−18.56
September–October	θ	0.87	2.94	1.34
	AIC	−13.81	−11.10	−9.92
October–November	θ	1.65	3.78	1.26
	AIC	−34.69	−17.03	−3.12
November–December	θ	5.13	12.59	2.49
	AIC	−109.38	−91.74	−62.23
December–January	θ	0.13	0.55	1.04
	AIC	1.43	1.55	1.62

. The bold AIC values in the table indicate the optimal copula function. The Frank copula was selected in January–February, February–March, and December–January; the Gumbel copula was selected in March–April, May–June, June–July, and July–August; and the Clayton copula was selected for other months.

The contour map of the cumulative distribution function of the optimal copula function for each adjacent month is shown in Figure 5. In Figure 5, it is easy to obtain the probability that the incoming water in a certain month is within a certain magnitude and that the adjacent month’s incoming water is within a certain magnitude. Through the joint distribution, the corresponding joint probability of the two adjacent months’ incoming water can be obtained. It should work as well.

4.1.3. Stochastic Simulation of Incoming Runoff

Combined with the lognormal distribution and copula function, the Gibbs method was used to simulate monthly runoff. The simulated and measured runoff of the NZD reservoir are shown in Figure 6. It can be seen from the figure that the simulated monthly runoff distribution was highly consistent with the measured runoff distribution. The maximum value appeared in August and the minimum value appeared in March. The maximum value of the simulated runoff in each month was higher than the measured value, and the months with the largest deviation occurred in the wet season, especially in August. However, from the median value, it could be found that the median value of the simulated runoff and measured runoff were basically consistent. The multi-year average of the simulated runoff (537.2 10⁸ m³/s) was 0.7 10⁸ m³/s less than the measured runoff (537.9 10⁸ m³/s). The standard deviation of the simulated runoff (73.27 10⁸ m³/s) was reduced by 4.71 10⁸ m³/s compared with the standard deviation of the measured runoff (77.98 10⁸ m³/s). According to the measured annual mean runoff of ±1 standard deviation as a reasonable range, the simulated annual runoff was statistically analyzed (in Figure 7). The proportion of the simulated annual runoff exceeding this range was 13.98%, and the proportion lower than this range was 14.9%. Most of the simulated runoff values were within this range, accounting for 71.12%, showing the uncertainty of the incoming water process.

4.2. Process of Incoming Runoff

4.2.1. Transfer Probability Matrix Based on Measured Runoff

The monthly runoff of the NZD hydropower plant was sorted in descending order, 10 representative runoff QR (QR1 representing the average of the runoff in $P\in (0,10\%]$, QR2 representing for the average of runoff in $P\in (10\%,20\%]$, etc.) were selected, and the runoff in each month under the same QR was plotted on the same curve (Figure 8).

Based on the monthly representative flow of the NZD hydropower plant (

Table 3. Representative runoff in each month for measured runoff of the NZD hydropower plant (10⁸ m³).

Month	Representative Runoff QR
	1	2	3	4	5	6	7	8	9	10
June	67.2	60.9	53.1	51.5	48.1	44.9	41.5	37.3	34.1	27.6
July	112.9	102.4	95.3	88.0	82.7	75.9	74.3	70.1	64.9	56.2
August	141.6	122.4	112.9	104.5	96.9	91.4	84.6	75.6	71.2	53.3
September	119.0	106.1	94.6	89.0	83.8	77.0	71.4	68.0	60.4	49.6
October	89.8	74.6	69.3	65.1	62.0	57.5	54.1	49.6	45.7	35.5
November	54.1	43.6	40.4	38.1	37.3	34.7	31.0	29.2	27.6	20.4
December	28.9	26.5	24.2	23.2	22.6	21.5	20.8	19.5	18.4	16.0
January	21.1	19.3	18.1	17.4	16.9	16.8	16.4	15.1	14.2	12.0
February	17.4	16.0	15.4	14.9	14.5	14.2	13.7	13.0	11.9	10.5
March	16.8	16.2	15.0	14.4	14.0	12.9	12.6	12.1	11.7	10.8
April	20.7	19.7	18.8	18.3	17.5	17.0	16.1	15.3	13.5	11.8
May	35.5	30.7	29.7	27.6	26.5	25.3	23.5	21.6	19.0	15.3

), the one-step transfer probability matrix of monthly runoff is derived. It could be seen from the transfer probability matrix (Figure 9) that because the measured runoff was limited (only 61 years from 1953 to 2013), the derived transfer probability matrix did not have obvious statistical regularity. For example, when the incoming runoff in June dropped into QR1, there were only three representative runoffs in July (QR1, QR4, and QR5), while the probability of incoming runoff for other representative runoffs was 0. The transfer probability matrix of runoff was the important input of the stochastic optimal operation model. However, the transfer probability matrix based on measured runoff had the phenomenon of too many 0 values, which made it difficult to obtain reasonable optimal operation results. Therefore, the Gibbs method was used to generate a large number of samples to deal with the impacts of insufficient data on the transfer probability matrix.

4.2.2. Transfer Probability Matrix based on Simulated Runoff

In this paper, 5000 groups of stochastic runoff of the NZD hydropower plant were simulated based on the Gibbs method to derive the transfer probability matrix. The representative runoff in each month is shown in

Table 4. Representative flow in each month for simulated runoff of the NZD hydropower plant (10⁸ m³).

Month	Representative Runoff QR
	1	2	3	4	5	6	7	8	9	10
June	65.9	59.5	54.5	51.0	47.9	44.8	41.7	38.6	34.4	20.0
July	108.7	99.3	92.9	87.8	83.3	79.0	75.1	70.5	64.1	39.8
August	139.0	122.4	112.7	104.9	97.7	91.1	84.7	77.7	68.7	36.2
September	114.6	103.0	94.9	88.4	82.8	77.8	72.9	67.5	60.8	36.6
October	83.7	75.4	70.0	65.5	61.6	57.9	54.4	50.3	45.5	28.3
November	51.4	45.5	42.1	39.1	36.7	34.4	32.0	29.4	26.4	15.9
December	29.0	26.7	25.1	23.9	22.9	21.9	20.8	19.6	18.0	12.6
January	20.6	19.4	18.5	17.9	17.2	16.6	16.0	15.2	14.3	10.0
February	17.3	16.3	15.6	15.0	14.5	14.0	13.5	13.0	12.2	8.4
March	16.8	15.8	15.0	14.5	13.9	13.4	12.9	12.3	11.5	8.1
April	21.1	19.6	18.7	17.9	17.2	16.5	15.8	15.1	14.1	9.6
May	34.9	31.6	29.5	27.8	26.2	24.8	23.4	21.7	19.6	10.2

, and the one-step transfer probability matrix from June to May of the next year was obtained (Figure 10). Moreover, the phenomenon of centralized probability distribution and the transfer probability of 0 did not appear during the process of deriving the transfer probability matrix based on the simulated runoff. It was concluded that the transfer probability matrix based on simulated runoff could solve the problem of deficiency in measured runoff data.

4.3. Result of Hydropower Plants Operation

In this paper, the transfer probability matrix of the simulated runoff of the NZD reservoir was used as the input of the reservoir operation model, and the SDP algorithm was used to solve the model to obtain the operation results of the NZD reservoir under 10 representative runoffs. At the same time, taking the incoming runoff of a normal year as an example, the operation results of the hydropower plant operation model based on the DP and SDP algorithms were compared and analyzed to verify the rationality of the SDP. In order to show the water level process of the NZD hydropower plant under different incoming runoff in the current period after being dispersed, the corresponding water levels at the beginning and end of the month of the same representative flow were drawn on the same curve to obtain the yearly optimal operation water level diagram (Figure 11).

Taking the measured data as the input of the operation model, the optimal operation rules of power generation and the corresponding stage power generation expectation under the discrete water level and the incoming runoff at the beginning of each month were solved (Figure 11). Since the starting and ending period of annual operation of the NZD hydropower plant was from June to May of the next year, and the decision water level of May was fixed (dead water level), the operation rule was a vertical line. Similarly, in June of the next year, because the initial water level was fixed, the decision water level was a horizontal line. For other months, as the water level at the beginning of the month of the NZD hydropower plant increased, the increase trend of the water level at the end of the month was obvious; and as the range changed from QR1 to QR10, the corresponding water levels at the end of the month of the same water level at the beginning of the month showed an increase tendency, which proved the rationality of the water level operation diagram. The incoming runoff of the NZD hydropower plant in flood season (July to September) was large and unstable, so the curves of the July, August, and September operation diagrams were scattered and irregular. In the dry season, the curve of the corresponding water level in the operation diagram was denser, the curve space was closer and closer, and the curves even overlapped together, which was regular. The optimal operation diagram was used to obtain the corresponding end water level of a certain beginning water level by determining the representative range of any incoming runoff, which was convenient for operation in practice.

In order to examine the rationality of the operation rules formulated by the SDP, the deterministic dynamic programming (DP) was selected to obtain the operation results as a compared case. Taking the incoming runoff process of the NZD hydropower plant in a normal water year as an example, the SDP and DP algorithms were used to solve the operation model. The operation processes of the SDP and DP are shown in Figure 12. The annual power generation based on the SDP and DP were 305.97 × 10⁸ $\mathrm{k}\mathrm{W}\cdot \mathrm{h}$ and 306.91 × 10⁸ $\mathrm{k}\mathrm{W}\cdot \mathrm{h}$, respectively. The total power generation and operation process of the SDP was same as those of the DP, which proved to be reliable for the SDP.

5. Conclusions

There are many uncertain factors faced during the actual operation of hydropower plants, among which the uncertainty of incoming runoff contributes the most. Therefore, taking into consideration the randomness of incoming runoff, more risk information about the actual operation of hydropower plants would be provided for decision-makers to make reasonable decisions. In this paper, the copula joint function was selected to construct the joint distribution of adjacent monthly runoff of the NZD hydropower plant, and the runoff simulation model was constructed based on the Gibbs method. While considering the randomness of incoming runoff, the stochastic optimal operation model was constructed and solved by the SDP algorithm.

The monthly incoming runoff of 5000 groups of the NZD hydropower plant were simulated based on the copula function and the Gibbs method. Setting the measured annual mean runoff as ±1 standard deviation as a reasonable range, the simulated runoff dropping into this range accounted for 71.12%, and those exceeding this range accounted for 13.98% and lower than this range was 14.9%, which showed the uncertainty of the incoming runoff.

The monthly runoff of the NZD hydropower plant was sorted, and 10 representative flows (frequency of 10%~100% with range of 10%) were selected. Due to the insufficiency of measured runoff, there were too many 0 values in the derived transfer probability, but after the simulated runoff series were introduced, the results significantly improved.

Taking the transfer probability matrix of simulated runoff as the input of the stochastic optimal operation model of cascade hydropower plants, the operation diagram containing the future-period incoming water information was obtained, which could directly provide a reference for the optimal operation of the Nuozhadu hydropower plant. In addition, taking the incoming runoff process in a normal year as the standard, the annual mean power generation based on stochastic dynamic programming was similar to that based on dynamic programming (respectively 305.97 × 10⁸ $\mathrm{k}\mathrm{W}\cdot \mathrm{h}$ and 306.91 × 10⁸ $\mathrm{k}\mathrm{W}\cdot \mathrm{h}$), which proved that the operation diagram constructed in this study was reasonable.