Finding Multiple Optimal Solutions to Optimal Load Distribution Problem in Hydropower Plant

Optimal load distribution (OLD) among generator units of a hydropower plant is a vital task for hydropower generation scheduling and management. Traditional optimization methods for solving this problem focus on finding a single optimal solution. However, many practical constraints on hydropower plant operation are very difficult, if not impossible, to be modeled, and the optimal solution found by those models might be of limited practical uses. This motivates us to find multiple optimal solutions to the OLD problem, which can provide more flexible choices for decision-making. Based on a special dynamic programming model, we use a modified shortest path algorithm to produce multiple solutions to the problem. It is shown that multiple optimal solutions exist for the case study of China’s Geheyan hydropower plant, and they are valuable for assessing the stability of generator units, showing the potential of reducing occurrence times of units across vibration areas.


Introduction
Worldwide, hydropower plants produce about one-quarter of the world's electricity and supply more than 1 billion people with electric power [1].For example, there are more than 11,000 hydropower plants operating in China, providing about 20% of the country total electric power and making hydropower the country's largest renewable energy source [2].It has been shown that improving the operation of hydropower plants can increase the total profit by 1%-3% [3], which is of great significance for making energy production more efficient and reducing greenhouse gas emissions.
Hydropower plant operation involves issues such as optimal load distribution (OLD) and unit commitment [4,5].In this paper, the OLD problem is defined as a deterministic optimization problem with a fixed hydraulic head and a single time step [6], namely economic dispatch.The objective is to minimize the total water discharge by optimizing the load distribution among multiple generator units.Since the OLD is a subproblem involved in unit commitment problems, it is a basic and vital module for a hydropower plant's economic operation and generation scheduling [4].
It is noted that the OLD problem does not take into account the following factors: the amount of water in the basin is limited and therefore the objective usually aims at finding the optimal production levels at each period of an optimization horizon (e.g., a week) in order to achieve the maximum benefits from the limited quantity of water available in the basin; the level of the water in the basin (i.e., the hydraulic head) influences also the amount of electric energy that could be produced by a given amount of water.
There are many situations in which multiple optimal solutions (MOS) (i.e., solutions with the same objective value but with different values for the decision variables) may exist.For example, a general non-convex function, such as sin(x), has multiple minima and maxima.For the OLD problem, the input/output (I/O) function of a generator is often non-convex, piece-wise linear with many pieces having the same slopes (Figure 1) [8].In this case, suppose x * 1 , x * 2 are the optimal loads for generators 1 and 2, and suppose x * 1 and x * 2 on the pieces with the same slopes, then we can freely adjust x * 1 and x * 2 by the same amounts but in the opposite direction while still achieving the same objective value and satisfying the fixed total load constraint.In fact, it was reported that there were MOS for the refill operation of the Three Gorges Reservoir in China [29].MOS have been studied in different contexts.For example, in integer linear programming, all optimal solutions are found by iteratively solving the original problem with new integer constraints added to exclude the optimal solutions found from the previous iterations [30].For general nonlinear programming problems, a dynamical trajectory-based methodology was developed for computing multiple local optimal solutions [31].In the shortest path problem, all the optimal and near-optimal solutions are found by using a near-shortest path algorithm, which is a modified version of backtracking in DP [32].For general optimization problems, the Niche PSO technique was developed to locate and refine MOS [33].The multiple near-optimal solutions, whose objective values are equal or sufficiently close to the optimum, were explored in the deterministic reservoir operation problems [34].However, as far as we know, the problem of finding multiple solutions has not been addressed in the OLD problem.In fact, we will show in the case study that MOS do exist in the OLD problem.
In the rest of this paper, we present the description for the OLD problem in Section 2, followed by our solution approach in Section 3. Using discretization, which is usually used to solve DP models, we reformulate the discretized OLD problem as a shortest path problem.We then develop a special shortest path algorithm to find MOS.The case study of the Geheyan hydropower plant is then shown in Section 4. Section 5 draws the conclusions of this research.
Since the discretization resolution affects the optimality of the solutions found, our method requires a reference or benchmark on the optimum as a guidance for the discretization.Therefore, though the focus of this study is on finding multiple solutions, we also provide models for finding an optimal solution in the Appendix as well.These include a MILP model to find the global optimal solution for the case of piece-wise linear I/O function and a Lagrangian relaxation model for general I/O functions.Additionally, we provides some cases with a nonlinear I/O function to illustrate that MOS might also exist in the OLD problem.

Problem Description
The problem of optimal load distribution among hydropower plant units is to minimize the total water discharge while ensuring the power system reliability [2,4,6].This problem is modeled as follows: where n is the number of generator units, x i and Q i (x i ) are the electric load and the water discharge for generator unit i respectively.L is the required total electric load for the hydropower plant.R i is the feasible load for generator unit i and is often equivalent to The minimum l i , maximum u i and infeasible range (z l i , z u i ) for unit i impose the operational constraints upon the real-time unit commitment and load distribution.The infeasible range (z l i , z u i ) is defined as the vibration area where the generator shaft tends to vibrate violently [35].It should be noted that this problem is independent of time.
The water discharge function Q i (x i ) is often non-convex and hence the OLD problem shown in (1) often has a non-convex objective function [36].In addition, the feasibility space is discontinuous due to the vibration range of the loads.Therefore, the OLD problem is often very difficult to solve even for finding just a single optimal solution.In fact, traditional methods such as Lagrangian relaxation, MILP or DP often can only find a near-optimal solution.On the contrary, our method will find multiple optimal or near-optimal solutions using discretization and a shortest path algorithm.

Solution Approach for Finding Multiple Solutions
We notice that the OLD problem has only one join constraint on the decision variables, i.e., n i=1 x i = L.This special structure allows us to consider applying dynamic programming algorithms.We first review the general DP framework in Subsection 3.1 to show the general idea.However, this general framework and traditional DP algorithms might not work for the OLD problem due to the non-convexity of the objective function and a discontinuous feasible region.This leads to the idea of using discretization that allows us to find not only a single optimal solution but a multiple of them.

Review of General Dynamic Programming (DP) Framework
DP, based on Bellman's principle, is a classic optimization method that improves computational efficiency greatly by saving some intermediate results into memory, which can be reused [37].A DP model can be formulated by the following procedures [14,15]: (1) Setting the decision variable x i as the load allocated to the i th generator.
(2) Setting the stage variable as the accumulative load of the first i generators: (3) Establishing the state equation as: (4) Define f i (s i ) as the minimum total water discharge when the accumulative electric load from the first generator to the i th one is s i , i.e., We also defined We can find f i (s i ) recursively using the following Bellman equation: where f 0 (s 0 ) = 0 is defined as the boundary condition.
In the cases when the I/O function is quadratic and there is no constraint on the loads, we can find the closed form for f i (s i ) and the optimal load distribution x i .However, this is not the case for the OLD problem since the discharge function is non-convex and the feasible range is discontinuous.We can also try approximate dynamic programming [38] to approximate f i (s i ).However, this method might be time consuming and inaccurate because of the discontinuous feasible region involves in this problem.
We will discretize the feasible space of the loads and show that the discretized OLD problem is equivalent to a shortest path problem that can be solved efficiently.The discretization resolution is adjusted such that the optimum total discharge found by the discretized OLD is relatively close to the optimum total discharge found using Lagrangian relaxation or MILP methods (see Appendix A).We modify the shortest path algorithm slightly to keep track of all the optimal paths.Although discretization generally provides only approximate solutions, near-optimal solutions are often accepted in practice.In addition, for this particular case study, we will show that the solutions obtained were actually the global optima solutions.

Discretized OLD as a Shortest Path Problem
The decision variable x i in the OLD problem is continuous and can take any value in its feasible domain R i .Let us discretize the interval [0, L] into discrete values {p 0 , p 1 , . . ., p K } with p j = j∆, ∀j = 0, 1, . . ., K and ∆ = L/K.Suppose we can impose a further constraint that x i ∈ {p 0 , p 1 , . . ., p K }, and consider a new optimization problem: where S i is defined to be the feasible set for x i in the discretized OLD problem, i.e., It is noted that an optimal solution of the new problem shown in Equation ( 5) might be different from that of the original (continuous) problem shown in Equation (1).In addition, the optimum of the new problem should not be smaller than that of the original problem.However, if the discretization resolution is fine enough and if multiple solutions do exist, the gap could be small.In fact, we will show in the case study that the discretized problem can produce optimal solutions.
We will show how the discretized OLD problem can be reformulated as a shortest path problem.In the two dimensional coordinate system, we draw coordinates (i, p k ) where i = {0, 1, . . ., n} and k = {0, 1, . . ., K}.Let (x 1 , x 2 , . . ., x n ) be a feasible solution to the discretized OLD problem, i.e., n i=1 x i = L and x i ∈ S i .In Figure 2, we draw the points (0, 0), (1, s 1 ), (2, s 2 ), . . ., (n, s n ) on the two dimensional coordination system, where s i = i j=1 x i .If we connect the points sequentially, we obtain a path from (0, 0) to (n, L).In other words, each of the paths from (0, 0) to (n, L) would be a solution to the total load constraint n i=1 x i ≡ L. We also set the cost of the arc from coordinate (i − 1, s i−1 ) to coordinate (i, s i ) is set as follows: Then, the total cost of the path would be the total discharge n i=1 Q i (x i ) when (x 1 , x 2 , . . ., x n ) is a feasible solution.The discretized OLD problem is equivalent to finding a shortest path from the origin (0, 0) to the destination (n, L).We then can apply different shortest path algorithms to solve the discretized OLD problem efficiently [39].Figure 3 shows that the discretized OLD problem becomes the shortest path problem from the origin (0, 0) to (n, L).The nodes of the network represent the accumulative loads, the edges represent the load distribution and the costs of the edges are the water discharges.
Figure 3.The discretized OLD is equivalent to a shortest path problem for finding an optimal path from the origin (0, 0) to the point (n, L) that passes through the points on the grid.

Finding MOS
The solution approach for finding MOS in the OLD problem is same as the standard shortest path algorithm.Let us define q(i, p k ) as follows: which is a discretized version of f i (p k ) as shown in the general DP framework in Equation (4).We define q(i, p k ) = ∞ if the problem is infeasible.
The idea of a shortest path algorithm is to solve the subproblems with small i and k and then use the results to solve subproblems with larger i, k until i = n and p k = L.In addition, this algorithm keeps track of all the optimal subpaths to produce multiple solutions for the shortest paths.
We have: This is basically the Bellman equation or the principle of optimality [37].This equation is equivalent to saying that the optimal path from the origin (0, 0) to (i, p k ) can be found by identifying the points among the set of (i − 1, p k − x i ), i.e., the best points with the minimum total cost as shown in Figure 3.The following shortest path algorithm shows the formal procedure for finding an OLD solution: • For i = 1: We set q(1, p k ) for k ∈ {1, 2, . . ., K} as follows: • For i = 2 to n: We set q(i, p k ) for k = {1, 2, . . ., K} as follows: The number of operations for each i that is in the range [2, n − 1] is, at most, the number of links between node (i, k) and (i + 1, k ), which is equal to K(K+1)

2
. The number of operations for i = 1 and i = n are at most K each.Thus, the total number of evaluations is at most (n − 2)K 2 + 2K.
To find of all the optimal solutions, we can modify the traditional shortest path algorithm slightly to keep track of the optimal subpaths.The idea is to use the principle of optimality in dynamic programming as follows: if an optimal path from (0, 0) to (i, p k ) contains the point (i − 1, p k ), then the set of all the optimal subpaths to (i − 1, p k ) is also contained in the set of all the optimal subpaths to (i, p k ).
Let D(i, p k ) be the set of all the optimal solutions of problem (7) at each i ≥ 1 and p k (each optimal solution is a vector of length i).Let {x * i1 , . . ., x * im } be the set of optimal solutions for the recursive problem (8).Then, with i = 1, we have: ∅ otherwise and with i ≥ 2, we have: The formal modified shortest path algorithm for finding MOS of the OLD is as follows: • For i = 1: We set q(1, k) for k ∈ {1, 2, . . ., K} as follows: ∅ otherwise • For i = 2 to n: We set q(i, p k ) for k = {1, 2, . . ., K} as follows: q(i, p k ) = min where {x * i1 , . . ., x * im } are the optimal solutions of problem (8).Otherwise, we set D(i, p k ) = ∅.
Notice that the algorithm works with all types of I/O functions and the discontinuous feasible region involved in this problem.The algorithm's computational complexity includes two parts: (1) the traditional computational complexity for the shortest path, which is (n − 2)K 2 + 2K; and (2) the computational complexity for procedure of trace back, M , which is the number of multiple optimal solutions.Since the M is not often a large number, the proposed approach could offer a practical solution in terms of computation time.

Multiple Solutions Space
The multiple solutions found by the special shortest path algorithm are scattered points in space.However, the space of all the optimal solutions could form continuous regions.In the case of piece-wise linear discharge functions, we can find this space by using marginal utility analysis.Among all the optimal solutions for which the constraints x i ∈ R i are not binding, the marginal utilities must be the same, i.e., ∂Q 1 ∂x 1 = ∂Q 2 ∂x 2 = • • • = ∂Qn ∂xn .Let x * i1 , x * i2 , . . ., x * im be a particular set of optimal unit loads.Then we can perturb these solutions at the same time such that their sum does not change.For example, we can adjust the optimal solution to x * i1 + , x * i2 − , where is small enough, while keeping all other unit loads unchanged.The range of is chosen so that the new loads x * i1 + and x * i2 − still lie within their original linear pieces of the piece-wise discharge functions.

Geheyan Hydropower Plant
As shown in Figure 4, the Geheyan hydropower plant is located about 9 km upstream of the Qing River from Changyang county, Hubei Province in China.The plant has a large hydropower generation capacity of 2.94 × 10 9 kWh per year with a firm output of 2.41 × 10 5 kW.The area of the Geheyan reservoir basin is 14,430 km 2 .Corresponding to normal and minimum water levels of 200 m and 160 m, respectively, the normal and dead water storage are 3.12 × 10 9 m 3 and 1.98 × 10 9 m 3 .The reservoir is used for inter-year flow regulation.Based on these records, a piece-wise linear approximating technique is used to describe the I/O curve.It very often uses this piece-wise linear I/O function in practical operation of hydropower plants [5].

Global Optimal Solution
The MILP model (Appendix B) is used to find the optimal solution (reference or benchmark), in which the I/O curve of generation unit is a piece-wise linear function and a penalty function method (by modification of the I/O curve) is used to avoid vibration.Lingo 9.0 [40] is used to solve the MILP model.The global optimal solutions are listed in Table 1 when the total electric load varies from 500 MW to 1200 MW.For simplicity, an assumption is made that the outputs of these units decline from unit 1 to unit 4.
Table 1.The results of global optimal solution to the OLD problem by using MILP.

Output of units (MW)
Total discharge The Lagrangian relaxation method (Appendix A) is also applied to find lower bounds of the optimal total water discharges.This method is often applicable in cases when the problem size is large (i.e., with large n and with more complex I/O function).The total water discharge, Q * L , are shown in Table 1.For example, a lower bound total water discharge of 517.1 m 3 /s is found when the load is set to 500 MW, which is quite close to the true optimum of 518 m 3 /s.We can also see that in some cases, e.g., with L = {550, 600, 900}, the optimal objective value of the OLD problem is equal to one obtained using the Lagrangian relaxation method.Nevertheless, both the MILP model and the Lagrangian relaxation method produce only a single solution.The focus of our paper is on finding MOS as we will show in the next subsection.

Multiple Solutions
The multiple-solution approach is implemented to find all solutions to the OLD problem of the Geheyan hydropower plant.The computation time is several seconds for a specific load.More than one solution is obtained, all of which have the same total water discharge as that of MILP.Some typical cases are described as follows.
(1) In the case that total electric load is equal to 500 MW (L = 500), the optimal solution consumes water discharge of 518 m 3 /s and distributes them to unit 1 and unit 2 (Table 1).
Table 2 shows single or multiple solutions under different discrete intervals.When the discrete interval of DP is 100 MW, a single solution is found but it is not optimal compared to the optimal water discharge found from MILP (518 m 3 /s).When the discrete interval is 10 MW, a single optimal solution is located.When the discrete interval decreases to 1 MW, multiple solutions are found and they all reach the optimal discharge.With a smaller interval (0.1 MW), even more optimal solutions are identified.
Table 2.The relationship between discrete the interval and solutions in load of 500 MW.The multiple solutions space can be found by the marginal utility principle,

Discrete interval
∂xn , e.g., ∂Q 1 ∂x 1 = 263−259 255−250 = ∂Q 2 ∂x 2 = 259−255 250−245 .Figure 6 shows the existence of MOS on a line segment.(2) When total load is set to 650 MW, four multiple solutions spaces are identified as follows.(3) When the total load is set to 1100 MW, four generation units are used.By using the multiple solution approach, each unit's output is either 295 MW, 285 MW, 275 MW, 265 MW or 255 MW.That is, if a solution sampled from above values satisfies the OLD constraints, it is an optimal solution.It is shown that the multiple solutions exist in the form of a few scattered points.

Non-Linear I/O Function
A second degree polynomial interpolation technique is used as a nonlinear approximation to the I/O curve.That is, the top three nearest points to (P i,k−1 , Q i,k−1 ), (P i,k , Q i,k ) and (P i,k+1 , Q i,k+1 ) are used to construct a quadratic function.This quadratic function is used to estimate the Q i (x i ).Based on the above approximation for the I/O curve, some cases are shown as follows.
(1) In the case that total electric load is equal to 500 MW, the optimal solution consumes water discharge of 518 m 3 /s and distributes them to unit 1 and unit 2. Indeed, the multiple optimal solutions are the same as that from the piece-wise linear approximation (Figure 6) because the nonlinear and linear approximation are the same within the range of [240 MW, 255 MW].(2) In the case that total electric load is equal to 650 MW, (296.7 MW, 286.7 MW, 66.6 MW) is the unique optimal solution.(3) In the case that total electric load is equal to 1100 MW, a few scattered points, the same as that of piece-wise linear I/O curve, are the multiple optimal solutions.
The modified DP algorithm is also used in an economic dispatch problem, EXAMPLE 3D in [4] (See Appendix C), to illustrate that the multiple solutions can also exist for a theme plant with nonlinear I/O functions.When a discrete step of 0.1 MW is used for the computation (programming with a 32-bit floating-point data), some multiple optimal solutions have been found, such as (725.It should be noted that these solutions are not the extract optimal but near-optimal solutions, because they depend on the data length of the computer.

Application
In the practical operation of a hydropower plant, the total load L often changes from time to time.In these cases, a small increment in load results in a new dispatch, which may end with the most economical values as determined by the model, but is not acceptable from an operational point of view [21].Typically, some MOS require one or more generators to drop or raise their output dramatically.We often want to avoid adjusting, especially through the vibration area since this may greatly affect production safety.Therefore, multiple solutions from OLD are valuable, since they provide more alternative schemes from which decision maker can choose as shown in the following examples: (1) Improving unit stability: Multiple solutions provide choices to reduce the readjustment efforts when certain conditions change.For example, the total load L often changes from time to time as shown in Table 3, which lists historical data from a practical operation of the Geheyan hydropower plant.If the decision maker has only a single solution, the solutions with different values of L

Finding an Optimal Solution Using Lagrangian Relaxation
Let β be the Lagrangian multiplier for the constraint n i=1 x i = L, the Lagrange dual problem become: This is equivalent to: For each fixed β, the objective function can be evaluated by solving n subproblems min x i ∈R i (Q i (x i ) − βx i ).Each of these subproblems is often very easy to solve since it has a single scalar variable x i .There are many possible ways to solve the subproblems depending on the specific form of the load function Q.For example, if Q is differentiable, a gradient descent method can be applied.The Lagrange dual problem can be solved using a subgradient method since at each fixed value β, we have efficient ways to evaluate the dual objective.The subgradient can be set as (L− n i=1 x i ), and indicates the degree of violation on the constraint ( n i=1 x i = L).The subgradient iterative algorithm can be employed and β can be updated as follows: where α is the gradient step.There are many ways and rules for setting α.The simplest way is to fix α with some values say 0.9.
In the case that I/O function is piece-wise linear, which is often true in practice, the subproblem can be solved easily by comparing the objective values at the ends of the linear pieces.Therefore, the Lagrange dual problem can be transformed into a linear programming problem and can be solved efficiently.In addition, if the I/O function is convex piece-wise linear, the optimal load distribution problem can be reformulated into a linear programming problem.

Finding a Single Solution Using Mixed Integer Linear Programming
We assume the I/O function Q i (x i ) is piece-wise linear and suppose the function is formed by K i discrete points (P i,k , Q i,k ), k = 1, . . ., K i .For any x i ∈ [P i,1 , P i,K i ], the discharge value Q i (x i ) is often described using conditional statements, i.e., we need to know the particular piece P i,k , P i,k+1 that x i belongs to, to apply the right linear function for Q i (x i ).In mathematical programming, these conditional statements can be replaced by binary variables.We will show that the pair of values (x i , Q i (x i )) can be represented by the variables (x i , q i , λ, y) that satisfied the following set of constraints:

Figure 1 .
Figure 1.An example showing that multiple solutions might exist in the OLD problem.

Figure 2 .
Figure 2. A solution to the discretized OLD problem can be viewed as path from the origin S(0, 0) to the destination E(n, L).

Figure 4 .
Figure 4. Map showing the location of the Geheyan Hydropower plant in China.

Figure 5 .
Figure 5. Efficiencies of generator in the Geheyan Hydropower plant for 110 m water head.

Figure 6 .
Figure 6.Multiple solutions to total electric load of 500 MW (output of unit 1 ≥ output of unit 2).

Figure 7 .
Figure 7. Multiple solutions space to total electric load of 650 MW (output of unit 1 ≥ output of unit 2 ≥ output of unit 3).