Dynamic Feasible Region Genetic Algorithm for Optimal Operation of a Multi-Reservoir System

Seeking the optimal strategy of a multi-reservoir system is an important approach to develop hydropower energy, in which the Genetic Algorithm (GA) is commonly used as an effective tool. However, when the traditional GA is applied in solving the problem, the constraints of water balance equation, hydraulic continuity relationship and power system load demand might be violated by the crossover and mutation operator, which decreases the efficiency of the algorithm in searching for a feasible region or even leads to a convergence on an infeasible chromosome within the expected generations. A modified GA taking stochastic operators within the feasible region of variables is proposed. When determining the feasible region of constraints, the progressive optimal approach is applied to transform constraints imposed on reservoirs into a singular-reservoir constraint, and a joint solution with consideration of adjacent periods at crossover or mutation points is used to turn the singular-reservoir constraints into singular variable constraints. Some statistic indexes are suggested to evaluate the performances of the algorithms. The experimental results show that compared to GA adopting a penalty function or pair-wise comparison in constraint handling, the proposed modified GA improves the refinement of the quality of a solution in a more efficient and robust way. OPEN ACCESS Energies 2012, 5 2895


Introduction
The optimal operation strategy of a multi-reservoir power system is important to develop and maximize utilization of hydropower energy [1][2][3][4].The topic has attracted more interest in recent years for its ability to explore the full potential of power generation by conserving energy or increasing power generation, which is achieved only through properly controlling reservoirs' water release and hydraulic head by optimal operation [5,6].
Seeking optimal operation strategies for reservoirs through the solution of an optimal operation model is the most important part of optimal operation of a multi-reservoir power system.Due to the complexity and nonlinearity of the reservoir system, acquiring the solution of an optimized model has been the focus of research [7].As a result, two typical types of algorithms have been developed, which are the mathematic programming approach and the artificial intelligence optimization approach.
Dynamic Programming (DP), which is based on the principle of optimality, is the most widely used traditional mathematic programming approach in determining the optimal strategy of a reservoir [8][9][10].However, it has computational burden when applied to a high-dimensional reservoir system, which is known as the curse of dimensionality.
With the development of intelligent algorithms and stochastic optimization theory, the Genetic Algorithm (GA) [11] has been widely used in reservoir optimization [12][13][14].However, there are two disadvantages of GA in reservoir operation.Firstly, the crossover and mutation operators might change feasible parents into infeasible offsprings by violating the constraints such as water balance equation, hydraulic continuity relationship, etc., thus weakening the algorithm's searching ability.Additionally, GA may terminate with an infeasible solution as the "best" chromosome in given generations if confronted with high-dimensionality cases or narrow feasible region situations.Secondly, it's hard to choose the proper constraint handling technique.
Among the various kinds of constraint handling techniques, penalty function is the most commonly used way [15][16][17].It distinguishes chromosomes by evaluating their degree of violation of constraints within penalty terms.However, the construction of penalty terms is difficult to perform because they can neither be too large to avoiding over-penalization, nor too small so as to cause under-penalization.Therefore, Deb [18] proposed the Pair-wise Comparison technique to avoid the difficulty of penalty term construction where comparing two chromosomes under some predefined rules to determine which one is superior.The rules are: a feasible chromosome is deemed superior to an infeasible one.In any pair of feasible chromosomes, the one with the larger objective function is superior.While in any pairs of infeasible chromosomes, the one with the smaller number of violations is superior.In order to keep the proportion of feasible and infeasible chromosomes balanced in the whole population, Runarsson and Yao [19] proposed a stochastic ranking technique to randomly rank chromosomes pair-wise by objective function or constraint violation.The repair operator approach offers a new perspective in constraint handling as to covert infeasible chromosomes into feasible ones.Chootinan and Chen [20] established a gradient method to repair chromosomes by reducing their constraint violations, but it requires the constraint functions to be continuous and differential, which limits its application in reservoir optimal operation.Michalewicz and Janikow [21] adopted a specially designed operator for ensuring feasible offsprings generated under the condition of a convex-feasible region, however, it's not suitable for non-convex problems.Teegavarapu and Simonovic [22] introduced repair strategies and heuristic rules to readjust solutions in helping refining their quality.
When applied to seeking the optimal strategy of multi-reservoir system, the traditional GA might be insufficient in generating feasible solutions due to the high dimensionality and complexity of constraints, which can result in getting a sub-optimal strategy.To overcome the traditional GA's shortcomings in searching for the feasible region, this paper presents a modified technique to determine a dynamic feasible region of constraints with the heuristics of reservoir operation and then execute evolving operators within the region for generating more feasible chromosomes.The proposed GA was applied in solving the joint optimal operation strategy of a multi-reservoir system consisting of five reservoirs.Compared with the traditional GA, the results showed that the modified GA is more efficient and satisfactory in seeking an optimal strategy.

Objective Function
Take the maximum of the total hydro-power generation in planning horizon as objective function: where N j,t denotes the power output of reservoir j during period t, j = 1,2,…,n, T is the whole planning horizon, and ∆t is time interval.

Constraints
(1) Water balance equation: , (2) Hydraulic continuity relationship: , , , ( ) (3) Power output function: (4) Storage volume limits: (5) Outflow limits: where V j,t and V j,t+1 are the initial and terminal storage volume of reservoir j in period t, respectively; Q j,t means the inflow of reservoir j in period t, q j,t is the outflow of reservoir j in period t; EL j,t is the sum of evaporation and leakage losses.Ω j refer to a set which contains the reservoirs located upstream of reservoir j; Qu k,t is the lateral inflow generated in the confined catchment between reservoir k and reservoir j; H j,t is the average hydraulic head of reservoir in period t; ( ) j f  is the power output function of a reservoir; , 1 j t V  and , 1 j t V  are the upper and lower limits of storage volume at the end of period t; , j t q is the minimal required discharge meeting the obligations of satisfying ecological and navigational demands in the downstream river channel and , j t q is total discharge capacity through spillways and turbines; , j t N and , j t N are the limits of power output; ND t is the minimal load demand imposed on the reservoirs system.ZE j is the expected terminal water level of reservoir j at the end of the planning horizon.

Structure of Traditional Genetic Algorithm in Multi Reservoir Operation
The Genetic Algorithm is a kind of population-based stochastic algorithm inspired by the process of biological evolution.It follows the principle of "survival of the fittest" to guide the chromosomes toward evolving and picks up the elitists.In an optimization application, these procedures lead the chromosomes to approach toward the global optimum gradually.The structure of the traditional GA in reservoir operation is discussed in detail in this section.
(1) Gene coding and chromosome initialization Water level is always preferred as a real-coded gene, which is denoted in a three-dimensional matrix Chr i,j,t , representing reservoir j's initial water level of period t. i = 1,2…,Pop, where Pop is the chromosomes' size.Chromosomes are often initialized randomly: , , , , , ( ) where , j t Z and , j t Z are water level boundaries, Rnd is a random variable which is subject to uniform distribution in (0,1).
(2) Chromosome evaluation criteria: Two types of criteria for evaluating the chromosome's superiority are discussed: (i) Fitness function Fitness function is most widely used in distinguishing the solutions' differences between chromosomes; the one with larger fitness is deemed as more fit.Usually, the penalty function is incorporated in the fitness function to handle the constraints: where Fit i is the Fitness function and Penal is the penalty term.Among the various forms of structuring a penalty term [23][24][25], penalizing according to the severity of constraints' violation is the most common way, which can be given as: where CountC is the count of violated constraints, INF is the penalty coefficient while Vl j (x) is the violated degree of constraint j.The process of calibrating a suitable INF usually involves trial-and-error and is rather time consuming.Under the situation of seeking optimal operation strategy of water level, the constraints which might be violated include power output limits, released flow limits and system load demand.Hence they are constructed into penalty terms of PenN j,t , PenQ j,t and PenSN t , respectively, with the purpose of reducing the fitness of chromosome that violates related constraints: , where INF 1 ~INF 3 are calibrated independently so as to consider the varied units of different degrees of violation.Finally, the fitness function can be supplemented as:

(ii) Pairwise comparison
To avoid calibrating suitable penalty coefficients, pairwise comparison is carried out to compare chromosomes' objective function E i and constraint violation Vio i , separately.According to Deb [18], Vio i is described as: (3) Crossover operator Crossover operator is a major stochastic searching approach through exchanging gene information between chromosomes as to breed offsprings.A single-point crossover as shown in Figure 1 is adopted: Assuming chromosome i 1 and i 2 are chosen to be crossed over at moment tp + 1. Exchange the two segments of gene after moment tp + 1 entirely, which can be expressed as: This crossover scheme minimizes the chance of violating at its best.It can be seen from Figure 1, at the time horizontal axis, that if both i 1 and i 2 are feasible, the offspring's water level at moments 1~tp − 1 and tp + 2~T are inherited from feasible parents which satisfies the constraints imposed on the singular reservoir, i.e., constraints ( 5) and (6).Similarly, at the vertical axis direction, an overall interchange of gene segments can ensure the offsprings meet the requirement of constraints imposed on a multi-reservoir system as constraint ( 7) at the same moments, only at the moment tp and tp + 1 may a violation occur because of random crossover.Setting the crossover probability to be 1, the child chromosomes ' Chr with population size of Pop are generated after crossover.
(4) Mutation operator Uniform mutation is applied as Figure 2 shown, under a mutation probability of pm to pick a mutation point, replace the former gene with a newly generated one: where Rnd and ' Rnd are independent random variables which subject to uniform distribution in (0,1).The child chromosomes '' Chr with population size of Pop are born after mutation.chromosome's score is set as follows: Pick up num chromosomes randomly from the whole population as rivals, the score equals the rounds that a specific chromosome i wins the competition.
Especially, if constraints are handled with a penalty function, chromosome i's performance is evaluated by a fitness function, so its' score Scr i can be counted by winning in fitness: If the pair-wise comparison is adopted in constraint handling, Scr i will be counted under the predefined rules: (6) Termination criteria When the best chromosome found has survived continually for Snum generations or evolution has taken place for Gen generations, the algorithm terminates and the best chromosome ChrBest is output.

Dynamic Feasible Regional Genetic Algorithm (DFRGA)
Random crossover and mutation in a traditional GA's formation cannot guarantee the offspring to meet the constraints and may change feasible parents into infeasible offsprings.To overcome this shortcoming, we determined a dynamic feasible region of major constraints and executed the operators in this region.
Suppose the crossover or mutation is to be carried at moment tp + 1 and water level (or storage volume) at moment tp are known.Under the requirement of satisfying constraint i during period tp, we can deduce moment (tp + 1)'s feasible region ).Likewise, supposing the state variable is known at moment tp+2, for satisfying constraint at period tp + 1, the feasible region FsbB i can be deduced reversely.Obviously, only those variables located at the intersection of FsbF i and FsbB i can meet both constraints at period tp and tp + 1; the region is denoted as Fsb i : For satisfying all the constraints, the dynamic feasible region SetFeasible is the intersection of Fsb i : The hydraulic continuity relationship in a multi-reservoir system results in the correlating acting of reservoir's operation strategy, which transforms related constraints into a multi-variables coupled situation, increasing the difficulty in determining the feasible region.Progressive optimality [26] is introduced to decouple the complex multi-variables: when the reservoir s is investigated, we hold the other reservoirs' storage strategy, and turn the multi-reservoir coupled constraint into a singular-reservoir constraint which only leaves the strategy of reservoir s unknown.When changing the storage of reservoir s at moment tp + 1, we keep other moments' storage unchanged and turn the multi-variable constraint into the singular-variable , 1 s tp V  's constraint.Finally, only the feasible region of a singular variable , 1 s tp V  is to be acquired.

Feasible Region of Singular-Reservoir Constraint
When V s,tp+1 's feasible region is investigated in singular-reservoir constraint, keeping the other reservoirs' storage unchanged.Then the feasible region can be determined as illustrated in Figure 3 (in a singular-reservoir constraint scenario): Under the condition of a singular reservoir, the violated constraints generally include constraint (5) (outflow limits) and constraint (6) (power output limits).
(1) Outflow limits constraint As the inflow Q s,tp and Q s,tp+1 are known as estimated by flow forecasting, storage V s,tp+1 's feasible region can be deduced from the water balance equation: s tp (22) (2) Power output limits constraint Hydropower output is a function of outflow; here we apply the water consumption rate function to illustrate the power output: )}/ ( ) where f j (q j,t ,H j,t ) is the power output function, ( ) QM  is the maximized discharge through hydraulic turbines, 3  m /s ; Fsb V  for satisfying the power output limits constraint is:

Feasible Region of Multi-Reservoir Constraint
Constraint (7) (system load demand) may be violated by the crossover or mutation operator in a multi-reservoir system.Figure 4 shows the schematic diagram of its feasible region: According to the water consumption rate function, Equation ( 26) can be expressed as: Generally, when the time interval is relatively short the hydraulic head H j,t has little variation, then the maximized discharge through turbines and average water consumption rate can be set as constants through interpolation.Equation (27) becomes: When reservoir s is investigated, reservoirs in the whole system can be classified into two groups according to whether their power output can be influenced by reservoir s's outflow.Group I, whose power outputs are independent of reservoir s's outflow, include reservoirs located upstream of s or in parallel connection with s.In Group II, the reservoirs' power output are influenced by s's outflow (including s itself), which will be denoted by set Ω s .For reservoirs in Group I, we record the sum of power outputs during period tp as Nsum tp since it only contains all known constants, then reservoirs in Group II (Ω s ) have to provide the least power output ND − Nsum tp to meet the load demand: It can be inferred from Equation (29) that the load demand constraint confined at power output initially can be transformed into an outflow related constraint.For reservoir s, seeking bounds of released flow , s tp q can be converted to seeking the variation range of , 1 s tp V  .
The upper limit of q s,tp is the maximized discharge through the hydraulic turbines Qm s .Therefore, the lower limit of V s,tp+1 can be figured out by water balance equation as . The lower limit of q s,tp is the minimal discharge for satisfying the minimal load demand ND tp − Nsum tp .Since the lateral flows and the hydraulic heads of reservoirs are all known, the total power output is only affected by q s,tp , and we adopt the trial and error approach to calibrate it, which generally includes the following steps: Step 1: Assume q s,tp = 0 (the minimal value of q s,tp 's lower bound is 0).
Step 2: Acquiring the sum of power output reservoirs in Ω s as Ncum tp , which can be expressed by: , min{ , } / s tp j tp j j j

Ncum
q Qm p    (30) If Ncum tp ≥ NDtp -Ncum tp , the upper limit of V s,tp+1 can be drawn with the water balance equation: Then terminate the calibration.Otherwise, go to step 3.
Step 3: Carry out statistics on the difference between maximized discharge through turbines and the outflow of the reservoir: Mark the reservoirs which have spillage through spillway in period tp as Ω spi , whose difference qDec j < 0. Reversely, the reservoirs haven't spillage in period tp are marked as Ω Unspi .
Step 4: Increase the outflow of reservoir s: Move to step 2. Similarly, the upper and lower bound of V s,tp+1 as ( ) ( ) FsbB V  can be determined through period tp + 1's constraint in a reverse way.Finally, the feasible region of constraint imposed on multi reservoir is the intersection of the two parts:

Modified Operators
To breed more feasible offsprings, we carry out the crossover and mutation in the feasible region SetFeasible as shown in Equation ( 21).
(1) Modified crossover operator: 2  (2) Modified mutation operator: ' is a random value in the feasible region.
Other procedures remain unchanged, and the fitness functions are applied to distinguish the differences between chromosomes.

Evaluation Indexes
To test the performance of the proposed model, we ran the stochastic algorithms Tr times independently in numeric experiments.Some statistic indexes have been used to evaluate the performances of algorithms, which include: (1) Electricity generation Average electricity generation: Maximized deviation: (2) Standard deviation of electricity: (3) Convergence ratio : where Slast i is consecutive number of generations when the optimum keeps unchanged in ith experiment until termination.
(4) Ratio of termination with feasible solution: where Ω Fsb is a set of feasible solutions.
(5) Average time consumption: where TCst i is the computation time of the ith experiment.
(6) Average feasible chromosome ratio: The average proportion of feasible chromosomes per generation in the ith experiment is: Then the statistical average proportion in total Tr times is: where η i,j is the proportion of feasible chromosomes in the jth generation of ith independent experiment, Gen i is the number of generations of the ith experiment.

Case Study
The multi-reservoir system of the Qingjiang cascades (consisting of the Shuibuya, Geheyan, and Gaobazhou reservoirs) and Three Gorges cascades (consisting of the Three Gorge and Gezhouba reservoirs), located at the middle part of Yangtze river in Hubei province of China, as shown in Figure 5, is the largest reservoir project system ever built in China till now.It is of great importance in relieving flood disaster and supplying clean energy, where it has the most typical structure as a complex and diverse system (i.e., it has various of regulation abilities, operational targets and districts coverage of energy supplying).Hence it is a representative case to explore its potential in power generation through optimal operation.To seek an optimal schedule in a 10 days' planning horizon with an interval of 1 day, under given conditions of forecasted inflow, lateral inflow and load demand from the power system, the constraints of the reservoirs are listed in Table 1.The reservoirs keep water along the deep-narrow valley and they are typical river-channel-type reservoirs.As the planning horizon is such a short period, leakage and evaporation losses in river-channel-type reservoirs can be negligible in this case.GA with penalty function (PFGA), GA with the pair wise compare approach (PCGA) and dynamic feasible regional GA (DFRGA) are used to seek the joint optimal operation strategy of the Qingjiang and Three Gorges multi-reservoir system separately.To analyze the sensitivity of algorithms to chromosome size, we investigated three groups of population size Pop = 50,100 and 150 respectively.Crossover probability is set to be 1, while mutation probability is 0.1.The rivals' size is set num = pop/2.Maximized generation Gen is set to be 100 and Snum = 5.The number of independent experiments Tr is set at 50.
All experiments are programmed by Visual Basic 6.0 under the WINXP operating system on a system equipped with an Intel core 2 duo E7500 CPU and 2G of RAM.After simulations, the evaluation indexes in Table 2 are drawn.Table 2 shows DFRGA is better than the other two schemes in all indexes related to electricity generation (objective function), which demonstrates that the modified algorithm performs significant improvements in algorithm's searching ability and robustness.We name the solution which has the median objective function among total Tr experiments as the "median solution".Figure 6 shows the feasible chromosome ratio variation of the "median solution" during evolution.The results of the experiments (Table 2 and Figure 6) indicate that: (1) The proposed DFRGA improved the optimizing ability of GAs greatly; it has increased electricity generated by 1.43% and 1.72% relative to PFGA and PCGA, respectively, under various size scenarios which may be attributed to DFRGA's better capability for searching the feasible region.Additionally, by decreasing the deviation in electricity by 83.94% and 85.23% on average relative to PFGA and PCGA, respectively, under various size scenarios DFRGA shows a more robust and steady performance.
(2) DFRGA has increased the proportion of feasible chromosomes during evolution significantly.Compared with other schemes, DFRGA has increased the feasible proportion by 60.17% on average.Searching in the feasible region assures a higher proportion of feasible chromosomes generated, which results in enhancing the convergence ratio by 71.33% and raising the ratio of termination with feasible solution by 26.33%.
(3) DFRGA has superior performance in breeding feasible chromosomes which validates the effectiveness of the core improvement.As Figure 6 shows, after several generations' evolution (about 20), DFRGA can reach a higher ratio of feasible chromosomes compared with the other two schemes, and keep a high level of enduring proportion.In contrast, PFGA and PCGA are so inferior in breeding feasible chromosomes and they may terminate in infeasible solutions when chromosome size is small.Besides, these two schemes must evolve at least 30 generations to provide a feasible proportion of no more than 30%.
(4) DFRGA has gained high efficiency in optimizing compared with the other two schemes.Although this algorithm introduces extra region estimation approaches, it consumes less time for an early convergence, as shown by Figure 6.
Figure 7 shows the power output series of reservoirs in the "median solution".Additionally, the average power output Navg denotes the average total power output of the whole reservoir system during the planning horizon: Figure 7 demonstrates that the DFRGA has achieved the largest average power output in the planning horizon among the schemes, and reached the largest objective function.Meanwhile, the DFRGA has attained a most uniform power output process which reflects the effectiveness of joint operation in a multi-reservoir system.

Conclusions
Seeking an optimal operation strategy of reservoirs is effective in developing hydropower and making full utilization of water resources.When a Genetic Algorithm is adopted in handling the problem, the constraints of water balance equation, hydraulic continuity relationship and power system load demand might be violated by a stochastic evolving operator, which decreases the efficiency of algorithm in searching for the feasible region or even lead to a convergence on an infeasible solution within the expected generations.These fatal defects to block the algorithm from reaching the global optimal strategy.To overcome these shortages, a dynamic feasible regional genetic algorithm (DFRGA) is proposed in this paper.It deduces the feasible region of constraints with knowledge of reservoir operation, and conduct evolving operators in the dynamic region so as to produce more feasible offsprings.
The Three Gorges cascade and Qingjiang cascade reservoir system is chosen as a case to validate the proposed algorithm.Compared with traditional GA, results indicate DFRGA has improved the solution quality by an increasing power generation by 1.43%, enhanced the robustness by decreasing the deviation by 83.94%, and shortened the converged generations and time consumption.All of these demonstrate DFRGA is an efficient and robust modified algorithm in reservoir operation.

Figure 3 .
Figure 3. Schematic diagram of feasible region in the singular-reservoir constraint scenario.

Figure 4 .
Figure 4. Schematic diagram of the feasible region in the multi-reservoir constraint scenario.
randomly generated values in the feasible region.tp + 1 is the crossover point, ( ) VZ  is the conversion function to convert storage into water level.

Figure 5 .
Figure 5.The location map of Qingjiang and Three Gorges multi-reservoir system.

Figure 7 .
Figure 7.The power output series of "median solution" in different GA algorithms. ' ( )

Table 1 .
Constraints of reservoirs in planning horizon.

Table 2 .
Statistical result of different GA's schemes in numerical experiments.