Planning Hydropower Production of Small Reservoirs Under Resources and System Knowledge Uncertainty

3.2.1. Stage one: understanding the system using a deterministic approach

In Figure 2, a flux diagram of a cascade composed of several reservoirs is shown. The deterministic modeling technique enables us to describe all water fluxes as shown in Figure 3 during every simulation stage.

Note that for each month the storage is calculated for each reservoir taking the difference of total inflows and total outflows. Total outflow is equal to the summation of discharged and spilled water, which is the release of the dam and will flow through to a downstream dam. Total inflow is equal to the summation of the released flow from an upstream reservoir and intermediate flows.

The storage equation is defined in a loop, where the storage at the end of time step “k” is dependent on the storage at the end of time step “k-1_L” the inflow during time step “k_L” and the turbine and spilled flow during “k_L” This is written in the following format:

for k= 1, …, Nand i= 1, …, M, S_i⁰ is given, Q_i⁰ = 0 and SP_i⁰ = 0

where S_i^kis reservoir storage volume of reservoir i at the end of period k, I_i^kis intermediate flow (inflow to the reservoir i during period k, apart from the release from an upstream reservoir). Accordingly, it can be seen in Figure 3 that every stage is composed of three decision nodes. Therefore, to completely determine a given stage, a set of Mequations is sufficient:

Parameters determined in these equations must be within real system design ranges. This means that each variable should respect the following constraints:

where S¯i and S¯iare lower and upper bounds of stored water of reservoir i, Q¯i and Q¯iare discharged water lower and upper bounds of reservoir i. Even though the volume is not forced to follow a certain reference, the last constraint ensures that the final volume of the reservoir is not below the initial value [21].

Equation systems (1) to (2) are sufficient to linearly describe the operational strategy of any cascade system with Mprojects. However, this problem can be simplified further by considering that every branch is represented by one state variable. This suggests a transformation of variable that still represents water fluxes dynamics. The next set of Eqs. (4)–(6) shows how to attain this for the first stage for a three-reservoir system, for example:

Eqs. (7), (8), and (9) are the equivalent for all remaining decision points:

The initial volume (S_i⁰) in each reservoir iand its inflows (I_i^k) are known, and together they build the initial data set for this problem. In order to apply these data, in an iterative optimization algorithm, they are represented in vector form as follows:

Consequently, equations are represented in matrices form; therefore, A represents equality constraints for state variables, and finally, Z is the vector of variables:

The problem formulation is completed by designing the optimization criteria (objective function). The objective function for hydropower energy maximization can be expressed as a product of the head for hydropower generation and the release. Therefore, nonlinear programming (NLP) may be considered for solving this problem. However, this research takes a different approach by applying linear programming to the linear operation model.

Theoretically, hydropower capacity of a storage plant installed at a reservoir can be expressed as

where ηρg, γcan be a constant value, ηis the overall, ρis the water density, and g is the gravity acceleration; Q(t) is the release for the power generation of at time tLand the net head H_n(t) can then be written as

where H(t) denotes reservoir water level, H_tail(t) is the water level of the downstream, and H_loss(t) denotes the head loss at a given instant t. If H_tail(t) and H_loss(t) are negligible in comparison with H(t), the approximation H_n(t) ≅ H(t) is applied, and the objective function for maximizing hydropower energy can be formulated as

where kis the discrete time and Q(k) and H(k) are the releases for hydropower generation and reservoir water level at time k. Nis the total time. Eq. (14) is the nonlinear product of the vectors Qand H. However, in this paper, a linear function is used based on the fact that Eq. (15) is maximized by maximizing each member of Q[Q(1), Q(2), …, Q(J)] and H[H(1), H(2), …, H(N)].

In turn, the nonlinear Eq. (14) can be replaced by the linear Eq. (16) as an objective function Zfor maximizing hydropower energy. Eq. (16) linearly sums up each element of Qand H.

where w_Hand w_Qare penalty factors on water level and the releases, respectively. Eq. (17) can replace Eq. (16) due to the fact that H(j) is directly proportional to S(j), that is, w_Hcan be substituted by w_Sas:

Eq. (17) is the linear combination of reservoir storage S(j) and release Q(j), and it can be represented as an alternative objective function for maximizing hydropower energy.

Subject to

The constraints are linear, the state equation is linear, and the objective function is chosen in linear form. The optimal solution can be obtained using various software tools readily available. Time lag or routing of flows between reservoirs is neglected in this formulation, which is reasonable for monthly time steps. Rainfall and net losses due to seepage, evaporation, and other reservoir losses are subtracted directly from the river inflows.

After formulating the problem deterministically, stage 2 follows which considers the uncertainty in the available resources. As shown in Figure 2, stage 2 is performed in two options according to the quantifiability of the available data.

3.2.2. Consideration of resources uncertainty by a probabilistic approach

To deal with LP under uncertainty, the chance constrained programming was introduced by Charnes and Cooper [22]. It extends the LP to enable the violation of the constraints to a certain extend. The reliability α∈ [0, 1] of not violating a constraint is specified by the decision maker, and thus, it allows for decision maker to directly control the level of risk he/she finds acceptable.

The deterministic reservoir model developed previously will be transformed here in the probabilistic form to deal with some uncertain inputs. The transformation to stochastic optimization is done through the introduction of an additional probabilistic constraint, which is shown below

where S˜ikis the random equivalent of S_i^k, the storage at the end of period k, S_{i, target}is the known decision maker specified target storage level of the reservoir, and α is the decision maker specified reliability of not violating constraint Eq. (19). It takes values between 0 and 1. This formulation is the realistic representation of the cascade control problem since in practice, bounds on storage are often violated by expansion of the conservation pool into the flood control pool. Adopting chance constraints simplifies the stochastic problem to deterministic. The random variables Ĩ_i^kand S˜i0are additive in consecutive periods, which makes the application of the convolution method feasible, whether they are linearly correlated or not [23]. The following procedure is used to do the transformation of the stochastic constraint.

Step 1.Eq. (1) is substituted into Eq. (19) to form:

where Ĩ_i^kis the random equivalent of I_i^k, the inflow during k= 1, …, Nand i= 1, …, M.

Step 2.A deterministic equivalent of the Eq. (20) is found by inversion and rearrangement leading to:

where for k= 1, …, Nand i= 1, …, M, Fxik−1(1−∝)is the inverse value of the cumulative distribution function of the convoluted Ĩ_i^k, evaluated at (1-α). Henceforth, it will be replaced by xik,1−α.

Step 3. The expression for two deterministic chance constraint time steps is given below,

for k= 1 and i= 1, …, M,

for k= 2 and i= 1, …, M,

Substituting for S_i¹

Eq. (21) can thus be expressed in final simplified chance constraint deterministic form as:

for k= 1, …, Nand i= 1, …, M.

It is important to note that the random variables inflow and initial storage are summed here. For the time interval k = 1, the sum is I˜i1+S˜i0, for k = 2 it is S˜i0+I˜i1+I˜i2,…, and for k = N it is S˜i0+I˜i1+I˜i2+⋯+I˜N. In the two-step algorithm developed by Simonovic to determine the best reliability levels for chance constraints, reliability levels are set in the first step, while the optimal open-loop reservoir operating strategy is determined by linear programming in step 2 [14]. This operating strategy sets the discharge from reservoir iin period tequal to x_tⁱ. The discharges x_tⁱ(for all iand t) are then used to transform the probabilistic constraints into deterministic ones. Strycharczyk and Stedinger [15] in their paper stressed out the problem with Simonovic and Marino [14] method that the discharge from a given reservoir iin a given period tis a constant, although the problem is stochastic. In stochastic reservoir management, the content of the reservoir is a random variable since it is fed by a streamflow, which is random. The content of the reservoir at start of period tcan take any value between the dead water level and the maximum reservoir volume. The optimal reservoir release for a certain period is a function of the reservoir content. Therefore, the open-loop operating strategy described in Ref. [14] is not quite acceptable for listing as a stochastic reservoir management problem. Therefore, in this paper we consider the initial reservoir storage level as a random variable also.

The random variables inflow and initial storage have known marginal probability distribution functions (PDF), f(I_k) and f(S_{0, i}), respectively, as a result of fitting a distribution to available historical data. However, the distributions of the sums have to be found. This is accomplished through a step-by-step iterative convolution method:

For k= 1, iterative convolution is used to reduce the two random variables Ĩ_i^kand S˜i0appearing in the equation of state to a single random variable. The new random variable is obtained as the sum of the random inflow and the random initial storage. Since their probability density functions are known, that of the new variable is obtained by convolution, which is expressed in discrete form by

where r_min ≤ r≤ r_max, p_i^Iis the probability distribution of inflow, p_J^Sis the probability distribution of the initial storage. The summation is carried out over the variable I, and hence, expression (26) yields:

The magnitudes of r_minand r_maxare found from r_min= i_min- j_max= a- dand r_max= i_max- j_min= b- cunder the constraints a≤ i≤ band c≤ j≤ d.

From k= 2 to k= N, the recursive equation for convolution can be expressed in genral form as in [24] as follows:

where r_min= i_min+ j_min, r_max= i_max+ j_max, under the constraints: a≤ i≤ band c≤ j≤ d.

The problem formulation becomes similar to linear formulation in the deterministic approach; Eq. (6) with the addition of the deterministic chance constraint and Eq. (25) can be solved with the same linear programming approach as the deterministic model formulation.

3.2.3. Stage two and three: consideration of resources and technical uncertainty using a fuzzy approach

If in stage 2 historical data are scarce, we can apply a fuzzy approach for both resources and technical uncertainty. The reservoir operation optimization model formulation will be expanded to utilize the fuzzy linear optimization approach, and in doing so, it will depart from the classical assumptions that all coefficients of the constraints need to be crisp numbers [25]. In the present section, a FLP formulation based on the work of [25] and further by Tanaka et al. [18] that considers both technological coefficients and resources characterized by uncertainty is presented in brief as follows:

The fuzzy version of the traditional linear programming optimization problem presented in Eq. (29) is:

The manager’s targets and system boundaries are the inequalities in the fuzzy system. The equation expressed that the manager’s targets can be lower than his/her desired level z₀. The same applies to the boundaries that they should be in the tolerance level b. The importance of targets and the constraints is set at the same level. For fully symmetric objective and constraints, Zimmermann [26] formulated the problem in simplified form as follows:

where B=[cA], d=[z0b].

The following expression for the (monotonically decreasing) linear membership function was proposed by Zimmerman for the j^thfuzzy inequality (Bx)_j≤ d_j.

where d_jand p_jare the desired level and the tolerance for violation of the j^thinequality, respectively. A j^thmembership function value of 1 denotes that the targets and constraints are fully satisfied. Violating the constraint outside the tolerance band, p_jgives a membership value of 0 and in-between values are linear. The membership function of the fuzzy set “decision” of model in Eq. (6) including the linear membership functions is shown below. The problem of finding the maximum decision is to choose x^* such that

In other words, the problem is to find the x^* ≥ 0 which maximizes the minimum membership function value. This value satisfies the fuzzy inequalities, (Bx)_j≤ d_jwith the degree of x* [18].

Substituting the expression (31) for linear membership function into Eq. (32) yields

The fuzzy set for decision can be transformed to an equivalent conventional linear programming problem by introducing the auxiliary variable λ:

subject to λ≤1+djpj−(Bx)jpj

It should be emphasized that the above formulation is for a minimization of the objective function and less than constraints and thus should be modified appropriately for other conditions.

Using the fuzzy optimization approach just described, and using the deterministic model given by Eq. (6) with modification for considering linear membership function for “greater than” constraints, the fuzzy formulation becomes:

Subject to

Expanding by substituting for (Bx)_j

all other constraints for k= 1, …, Nand i= 1, …, M