Published online 27 February 2007 in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/hyp.6539
Using a hybrid genetic algorithm – simulated annealing
algorithm for fuzzy programming of reservoir operation
Yu-Chen Chiu,
1Li-Chiu Chang
2and Fi-John Chang
1*
1Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei, Taiwan, ROC 2Department of Water Resources and Environmental Engineering, Tamkang University
Abstract:
We present a novel approach for optimizing reservoir operation through fuzzy programming and a hybrid evolution algorithm, i.e. genetic algorithm (GA) with simulated annealing (SA). In the analysis, objectives and constraints of reservoir operation are transformed by fuzzy programming for searching the optimal degree of satisfaction. In the hybrid search procedure, the GA provides a global search and the SA algorithm provides local search. This approach was investigated to search the optimizing operation scheme of Shihmen Reservoir in Taiwan. Monthly inflow data for three years reflecting different hydrological conditions and a consecutive 10-year period were used. Comparisons were made with the existing M-5 reservoir operation rules. The results demonstrate that: (1) fuzzy programming could effectively formulate the reservoir operation scheme into degree of satisfaction ˛ among the users and constraints; (2) the hybrid GA-SA performed much better than the current M-5 operating rules. Analysis also found the hybrid GA-SA conducts parallel analyses that increase the probability of finding an optimal solution while reducing computation time for reservoir operation. Copyright 2007 John Wiley & Sons, Ltd.
KEY WORDS genetic algorithm; simulated annealing; hybrid GA-SA; reservoir operation; fuzzy programming
Received 5 September 2005; Accepted 28 April 2006
INTRODUCTION
Water shortages in Taiwan are a regular occurrence owing to a monsoonal rainfall of wet and dry periods. Because approximately 75% of the rainfall is concentrated in typhoons occurring from May to October and this rainfall runs off in a matter of days, water supplies depend on reservoir retention. An increasing population and cou-pled economic growth requires advanced management techniques to assure water is available when needed. There are 66 reservoirs in Taiwan, with most classified as small or middle sized, having a storage capacity of less than 5 ð 106 m3. These reservoirs are experiencing a serious reduction in storage capacity caused by sedi-mentation related to landslides during typhoons. Unless reservoir operation can be improved, water shortages are expected to occur more frequently. Reservoir operation has emerged as a critical issue in the sustainable utiliza-tion of water resources in Taiwan.
Reservoir operations in Taiwan must consider rainfall seasonality and typhoon influences, multiple water uses, and legal and other constraints. Because the decision-making process involves uncertainty and inaccuracies, water resource managers mainly rely on operation rules, which were obtained through simulation and optimization techniques. Applications of optimization in reservoir operation are not new. Dynamic programming (DP),
* Correspondence to: Fi-John Chang, Department of Bioenvironmental Systems Engineering and Hydrotech Research Institute, National Taiwan University, Taipei, Taiwan, ROC. E-mail: changfj@ntu.edu.tw
originally developed by Bellman (1957) and applied by Hall and Buras (1961), has been used to solve the problems of reservoir operation. Many methods have been proposed to improve and/or implement DP for reservoir operation (e.g. Datta and Houck, 1984; Trezos and Yeh, 1987; Archibald et al., 1996; Chang et al., 2002; Chaves et al., 2003). The major disadvantage of DP is that computation times grow exponentially with the dimension of the state space (number of decision variables). Dimensionality issues limit the use of DP in exploring an optimal long-term operating schedule for effective water distribution.
In recent years, new methods of optimization have emerged. Search methods based on natural evolution mechanisms, such as genetic algorithms (GAs) and sim-ulated annealing (SA) have been shown to have advan-tages over classical optimization methods and have become widely used for solving a number of hydrolog-ical and water resource problems (Oliveira and Loucks, 1997; Pardo-Ig´uzquiza, 1998; Wardlaw and Sharif, 1999; Chang and Chang, 2001; Chang et al., 2005). GAs are based on the process of genetic change in liv-ing organisms, whereas SA is based on thermody-namic principles. The main advantages of these algo-rithms are their flexibility, ease of implementation, broad applicability and the potential of finding a near-optimal solution for many problems. Both GAs and SA support general-purpose optimization methods. GA approaches provide a global search through a system-atic and parallel manner. SA search methods, which
can theoretically converge to the global optimum solu-tion with unit probability (Kirkpatrick et al., 1983), pro-vide a local solution. Nevertheless, both methods usually required longer computation times (Wang and Zheng, 1998).
In this study, we propose the use of fuzzy programming to formulate the reservoir operation for multiple water uses in different hydrological conditions and a hybrid GA–SA algorithm for optimizing the operation schemes. The idea of constructing a hybrid GA– SA algorithm is mainly inspired by combining these two algorithms in order to facilitate the exhaustive and parallel treatment of optimization problems, to increase the probability of finding an optimal solution, and to reduce the computa-tion time required to reach a solucomputa-tion.
METHODOLOGIES
Reservoir operation in Taiwan must consider hydrological conditions, especially typhoon influences, multiple water uses, and other constraints. To deal with uncertainty and inaccuracies in the decision-making process of reservoir operation, fuzzy programming is proposed to formulate the problem into a degree of satisfaction ˛ among the users and constraints. The hybrid GA–SA algorithm is then used to search the optimal reservoir operation. The GA can perform a global search (population) and the SA algorithm can perform a local search (single point). To demonstrate this idea, the GA, SA and hybrid GA–SA algorithms were used to search 12 monthly water releases of Shihmen Reservoir in Taiwan. Monthly inflow data for three years reflecting different hydrological conditions were used to investigate these algorithms’ performances (efficiency and effectiveness). For further evaluation of the hybrid GA– SA method in searching for long-term reservoir operation performance, we implemented the method over a consecutive 10-year period.
Because the M-5 rule curves are the official way for guiding reservoir operation, it is crucial to learn about its performances in various hydrological conditions as a base for comparison. Consequently, a reference operational scheme was identified from the existing M-5 operational rules. In the following sections, fuzzy set theory, which is used to transform the objective and constraints of reservoir operation into fuzzy programming, is first briefly introduced and the optimizing algorithms (GA, SA, and hybrid GA–SA) are given.
Fuzzy set theory
Because reservoir operation involves targets and values that are often subjective and change with time, evaluating all objectives and constraints is a complex undertaking, if not impossible. To deal with this uncertainty, fuzzy set theory was used to identify targets and values that were then subject to natural system search techniques. Fuzzy set theory, developed by Zadeh (1965), provides a mathe-matical basis for representing linguistic vagueness. In the last decade, the area of fuzzy set theory has seen a surge
in research, particularly in hydrology, where research has focused on rainfall–runoff relationships (Maskey et al., 2004), rainfall forecasting (Yu et al., 2000), and reser-voir operation (Chang et al., 2002; Dubrovin et al., 2002, Mousavi et al., 2004).
Bellman and Zadeh (1970) developed a fuzzy program-ming approach to address optimization problems with fuzzy objectives and constraints to develop a decision-making model. In this approach, the fuzzy objective and constraints are transformed through membership func-tions into crisp objectives and constraints (Zimmermann, 1976; Wang, 1997). The optimum solution in fuzzy pro-gramming will be the intersection of the degrees of mem-bership of the objective Gxand constraints Cxand is defined as follows:
Dx D Gx ^ Cx DminGx, Cx 1 Because many possible solutions are possible, fuzzy set theory, specifically fuzzy programming, was used to determine the optimal solution with the maximum membership value:
Zx Dmax Dx 2 The fuzzy programming can then be expressed as
CTX ½ Z AiX bi Xi ½0
3
where X is the decision vector, CT and Ai are objective and constraints matrices respectively, bi is the maxi-mum desired levels, and Z is a given aspiration level. The degrees of membership of the objective 0x and constraints ix can be defined by the following mem-bership functions: 0x D 1 CTX > Z 1 Z Cp TX 0 Z p0C TX Z 0 CTX < Z p0 4 ix D 1 AiX < bi 1AiXbi pi biAiX biCpi 0 AiX > biCpi i D1, 2, . . . , m 5 where p0 and pi are tolerances of the objective and the ith constraint.
Let ˛ represent the satisfaction of the objective and constraints; the above problem can then be integrated as a maximizing satisfaction and transformed into the following mathematical programming:
max ˛ s.t. 0D1 Z CTX p0 ½˛ i D1 AiX bi pi ½˛, i D1, 2, . . . , m 6
˛, i 2[0, 1] X ½0
Therefore, the task of fuzzy programming is to determine ˛, the degree of satisfaction for the problem.
The optimum solution in fuzzy programming will then be searched by the GA, SA, and hybrid GA–SA optimizing algorithms.
Genetic algorithm methods
The GA is a search and optimization technique based on the mechanics of natural selection and genetics (Hol-land, 1992). According to Goldberg (1989), the major differences between GAs and traditional optimization methods are: (i) GAs use an evolving coding parame-ter set, not the parameparame-ters themselves; (ii) GAs search a population of points, not a single point; (iii) GAs use probabilistic transition rules, not deterministic rules. The GA is an adaptive and robust method for searching the optimum solution to a complex problem, although it may not necessarily lead to the best possible solution.
Reproduction. Reproduction is the procedure by which
chromosomes are chosen in the next generation. A popu-lar approach is fitness proportionate selection (Goldberg 1989), in which the probability Piof an individual i being selected is given by PiD fi n iD1 fi 7
where fi is the fitness of i and n is the population size.
Crossover. As in living organisms, the GA attempts
to preserve useful mutations so that, as new strings are created, the best are preserved. The crossover operator exchanges building blocks between two strings that perform well. The number of strings in which material is exchanged is controlled by the crossover probability, usually in the range of 0.5–1.0. In this study, the blend crossover (BLX-˛) (Eshelman and Schaffer, 1993; Chang and Chen, 1998) is used.
Mutation. In living organisms, mutation changes the
characteristics of genetic material in a chromosome to sustain genetic diversity in the population. Mutation is made occasionally with small probability. A random position of a random string is selected and is replaced by another character from the alphabet; for example, in the binary coding, this simply means changing ones to zeros and vice versa. A mutation operator of the real-coded GA is that specific element of a chromosome randomly that randomly jumps in the search space if it has an exactly equal chance of undergoing the mutative process.
Simulated annealing
SA is analogous to the annealing process used in metallurgy, where a metal object is heated to near its
melting temperature and then cooled slowly. This allows metal atoms to align themselves, crystallize and attain a minimum energy state. The SA algorithm, first proposed by Kirkpatrick et al. (1983), is a multipurpose heuristic searching technique for optimizing functions of many variables, and will evaluate a sequence of local optima in search of the global optimum. The SA algorithm superimposes a neighbourhood structure on the finite but large space of feasible solutions. Given a current solution Xcur, a candidate solution Xcan is drawn randomly from the corresponding neighbourhood. This new solution Xcan will be accepted to substitute the original solution Xcur subject to either improvement of the objective or a random experiment with an acceptance probability given by the Metropolis rule expC/ˇ, where C D CXcan CXcuris the difference of the cost function values and ˇ is the control parameter corresponding to temperature. According to the Metropolis rule, the SA would accept a worse solution in a given fraction of cases, and the fraction of worse solutions gradually approaches zero as the temperature decreases. If the new solution is accepted, then it will become the current point; otherwise, a new starting point from the current one is attempted. The number of searching points at a temperature can be set and seen as the length of the Markov chain.
The solution points in SA satisfy the Boltzmann dis-tribution, so that at any given temperature the probability of obtaining a lower cost ending point is higher when the iteration is sufficiently large. However, in any imple-mentation of the algorithm the Markov chain is of finite length. Therefore, asymptotic convergence can only be approximate.
The main advantage of SA is its general ease of solu-tion. Solving a specific problem by SA, the user needs to determine the annealing schedule TK, control parameter ˇ, length of Markov chain Lmax, and a neighbourhood structure.
Hybrid genetic algorithm–simulated annealing algorithm
The GA– SA process consists of coupling a GA with SA in an iterative manner. In the initial GA phase, genetic operators are used to identify an interim solution, which is then used in SA analysis. The GA generates a set of initial random solutions and then uses selection, crossover and mutation operators for few generations Ninitial.
After Ninitial generations, all of the latest solutions are sent to the SA for further improvement. Because SA is a single-insertion neighbourhood scheme, the SA is executed for each solution from the GA. Once the SA is performed for all solutions in the latest generation of the GA, the best solutions obtained from SA are the solutions of the GA for the next generation. The GA and SA exchange continues until the required number of cycles is reached or the system has frozen. A flow chart of this two-stage searching technique is shown in Figure 1. Several important factors which could significantly influence the searching effectiveness and
Figure 1. The flow chart of the hybrid GA–SA algorithm
efficiency through the searching process are given as follows.
Changed temperature. If the initial temperature T0 is set to a high value, then the search route will result in high degree of randomness and will produce many redundant iterations in the course of SA (Dougherty and Marryott, 1991). Such a strategy is not used in our GA approach. Instead, the selection scheme of the GA is used. The hybrid algorithm starts with a randomly generated set of solutions (population X0). Next, the GA is applied for few generations Ninitial to produce a new set of solutions (population X1) by rejecting the higher cost offspring so that the average cost must be less than the random cost
X0. Then, the SA is performed from X1 with a given length Lmax of iteration for each solution until the next generation X0
1is created. The initial temperature is set as a small value, T0D0.01, and the subsequent temperatures are set as suggested by Lin et al. (1993):
TkD
T1 k D1 T1k1 k ½2
8
where T1 Dhighest cost lowest cost/(populationsize/ 2), k is the temperature stage and D 0.8 is the decre-ment ratio of temperature.
Stopping rule (system frozen). When the average cost
of each generation is about equal to that of the best gen-eration, the system tends to convergence (or is frozen). In this study, we set the lowest temperature equal to 0.001. As Tk is less than this temperature, the system is assumed to converge and the computing process can then be stopped. In fact, the objective function does not fluctuate widely in computing process; the search process was terminated after five times of temperature decrease.
APPLICATION
Reservoir system
The metropolitan area of Taipei, population 6.5 mil-lion, is located on the Tanshui River and its tributaries in northern Taiwan (Figure 2). A major tributary of the Tanshui River is the Tahan River, on which Shihmen Reservoir is located. Shihmen Reservoir is managed for multiple purposes, including municipal and industrial water supply, agriculture, flood control, and power pro-duction. The operations of Shihmen Reservoir are based on the M-5 rule curves, as shown in Figure 3. The upper, lower, and critical limit curves are used for irrigation, power generation, and flood control respectively. Based on the operating rule curves, the operating policy can be briefly described as follows.
1. When the water level is above the upper limit, the water release for hydropower generation should be increased to keep the water level near the upper limit. 2. When the water level is between the upper and lower limits, the release, including public use, irrigation water supply and hydropower generation, is under normal operating conditions.
Figure 2. Location of Tanshui River basin, Taiwan
0 50 100 150 200 250 9 1 0 2 3 4 5 6 7 8 10 11 12 month
effective storage (million m
3)
upper lower critical Figure 3. The M-5 rule curve of Shiman Reservoir
3. When the water level is between the lower and critical limits, public use and irrigation water can be supplied as usual, but the water release for hydropower generation must be cut back by 20%.
4. When the water level is below critical limits, irrigation water must be cut back by 30% and the water release for hydropower generation must be reduced to the magnitude of public use plus irrigation water requirements only.
It appears that the operating rule defines the release within each year as a function of the existing storage level and overall release target amounts.
The monthly inflow of Shihmen Reservoir for three different hydrological years (1988, 1990, and 1993) was used in this analysis. Table I summarizes the monthly inflow and future water demand (i.e. 2010) where the demand was estimated based on the recorded monthly water use in 2001 multiplied by 1.1.
Modelling
This analysis focused on storage and release, where storage is provided for flood control, power generation and recreation and releases are made to meet water supply demands. The dominant demands are for agriculture and
Table I. The monthly inflow and water demand Monthly inflow (106m3) Monthly demand
(106 m3) 1988 1990 1993 Jan 37Ð24 59Ð47 34Ð12 65Ð58 Feb 36Ð43 46Ð67 24Ð61 92Ð00 Mar 52Ð38 78Ð76 70Ð94 119Ð37 Apr 133Ð77 273Ð45 96Ð24 107Ð21 May 82Ð94 100Ð32 84Ð22 109Ð11 Jun 87Ð73 219Ð50 159Ð29 108Ð82 Jul 40Ð12 111Ð68 69Ð76 132Ð68 Aug 62Ð67 637Ð19 70Ð90 131Ð63 Sep 114Ð88 571Ð39 34Ð24 112Ð34 Oct 258Ð45 76Ð47 40Ð97 111Ð39 Nov 92Ð73 48Ð57 27Ð06 78Ð12 Dec 36Ð35 33Ð00 26Ð14 50Ð37 Total 1035Ð69 2256Ð47 738Ð49 1218Ð62
public water supply use. The objective function and constraints were formulated as:
Objective function. min 1 5 12 iD1 SiXi Xi 2 C4 5 12 iD1 RiDi Di 2 ni 9 max ˛ D max 0, min 1 1 1p 0 1 5 12 iD1 SiXi Xi 2 C 4 5 12 iD1 R iDi Di 2 ni 1 Sibi pi i D5, 6, 7, 8, 9 11
where Siis the effective storage of the ith month, Riis the the release of the ith month, Xi is the the low rule in the M-5 curve of the ith month, Diis the downstream demand of the ith month, and ni is the number of cumulative shortage months.
Constraints.
water balance RiDSi1SiCqi flood seasons 0 Si XCi other seasons 0 Si Smax
0.9S1 S121.1S1
for i D 5, 6, 7, 8, 9
where qi is the monthly inflow of the ith month, XCi is the upper rule in the M-5 curve, and Smaxis the effective storage of the reservoir.
As mentioned above, Shihmen Reservoir is a mul-tipurpose reservoir; however, if one looks through the operation policy, it is easy to see that the main considera-tion is water supply for a different purpose. Consequently, the relative weights for the purpose of water storage ver-sus water release are set as 1 : 4, as suggested by Chang and Chen (1991). Because reservoir operation involves targets and values that are often subjective and change with time, evaluating all objectives and constraints is a complex, if not impossible, undertaking. To deal with this uncertainty, fuzzy theory may provide the most appropri-ate methodology for modelling reservoir operation.
The fuzzy programming is rewritten as max ˛ s.t. 1 5 12 iD1 SiXi Xi 2 C4 5 12 iD1 RiDi Di 2 ni Z C 1 ˛p0
water balance RiDSi1SiCqi for i D 1, 2, . . . , 12 flood seasons SibiC1 ˛pi for i D 5, 6, 7, 8, 9 other seasons 0 SiSmax for i D 1, 2, 3, 4, 10, 11, 12 0.9S1S121.1S1
0 ˛ 1 10
Z is a given aspiration level and is set to be zero, p0 is the tolerance of objective and is set as 12, and bi and pi are the maximum elevation at which water should not be higher than this level and the tolerance of the ith resource constraint respectively. Let bi be the upper curve of the M-5 rule curves and let pi be the difference between upper curve and full storage in this study. Therefore, the task of fuzzy programming is to determine ˛ as
For the purpose of comparison, the currently used M-5 operating rule curves for the operation of this reservoir are also evaluated.
Results
The performance of the proposed algorithms and the M-5 rule curve for the operation of Shihmen Reservoir were evaluated by using the degree of satisfaction ˛ and the generalized shortage index GSI (Hsu, 1995), given by GSI D 100 N N iD1 DPDi 100 ð DYi 2 12 where N is number of sample years (N D 1), DYi is days of the ith year (DYi D365 or 366) and DPDiis the deficit percentage index of the ith year, given by
DPD DDDR% ð NDC 13
where DDR is the daily deficit rate and NDC is number of days in a continuous deficit. In our case, because we use monthly data, the DYi is changed to the number of months in a year, i.e. DYiD12, and DDR is changed to represent the monthly deficit rate and NDC is changed to represent the number of months in a continuous deficit. The GSI is an indicator of the social tolerance limits to water shortage that measures the yearly deficit rate. Under this criterion the model that produces the minimum GSI is the one with the better performance.
Table II gives the values of the parameters used in the GA–SA. These parameter values are also used in GA and SA analysis. Because the number of generations Ngin the GA phase could influence the convergence and quality of final solutions, Ng was determined in the following manner.
(1) The GA is performed for the first 30 generations with its maximum, minimum, and averaged values of the objective in each generation recorded (i.e. population size of 300). The results are shown in the Figure 4. It appears that the values of the objective are increasing rapidly in the first 10 generations and
Table II. Parameters used in the GA– SA
Algorithm First phase Second phase
GA SA
Parameters set Population size S 300 Length of Markov chain Lmax 100
Crossover probability Pc 0Ð8 First stage
Mutation probability Pm 0Ð05 Initial temperature T0 0.01°C
Number of generations Ng 10 Other stage
Initial temperature T1 (highest cost lowest cost)/(300/2)
Decrement ratio of temperature 0.8
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 25 30 objective value a Generation number (a) 1988 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 25 30 objective value a Generation number (b) 1990 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 5 10 15 20 25 30 objective value a Generation number (c) 1993
maximum objective value average objective value minimum objective value
maximum objective value average objective value minimum objective value
maximum objective value average objective value minimum objective value
Figure 4. The performances of the first 30 generations by the GA
increasing slightly after that. As a result, the number of generations Ng in the GA phase was set to less than 10.
(2) Let X0
k be the set of solutions obtained by SA in the kth stage and used to perform the GA operation. After Ng generations of the GA, we obtain a new set of solutions denoted as XkC1.
Let ˛X0
k and ˛X 0
kC1 be the average cost of X 0
k and XkC1 respectively.
The degree of improvement I in the average objective is defined as I D˛XkC1˛X 0 k ˛X0 k 14 As mentioned above, the values of the objective usu-ally are rapidly increasing in the first 10 GA generations and slightly increasing after that. In order to make the hybrid GA– SA procedure more efficient and effective, the following strategy is proposed:
if I ½ 1.5 or NgD10, then stop the operator when k D1first stage;
if I ½ 0.05 or NgD10, then stop the operator when k >1.
The results are summarized in Table III and Figures 5 and 6 for three different hydrological conditions. From these results we note that the GA, SA and GA–SA provide consistently superior performance for both degree of satisfaction ˛ (above 60%) and the GSI index when compared with the M-5 rule curve. The hybrid GA–SA method gives the highest satisfaction and lowest GSI values in all cases. On close examination of the results, we found that the traditional M-5 rule curves (the official way for guiding reservoir operation) did not provide suitable operating schemes in the three years investigated, where the GSI values in 1988 (average hydrological condition) and 1993 (drought year) are very high and the degree of satisfaction ˛ D 0 in 1990 (wet year) and 1993. The results suggest that the M-5 rule curves should be adjusted to avoid a continuous water deficit in a drought year (even an average hydrological year) and a high risk of dam safety in a wet year. Table III also shows that there is high correlation between the degree of satisfaction and the GSI value: the higher the degree of satisfaction, the lower the GSI value. This result provides good evidence that the degree of satisfaction proposed in this study is a suitable criterion to evaluate reservoir operation with regard to water deficit and dam safety.
Figure 5 provides a comparison of methods, GA and SA, for the 1988, 1990 and 1993 hydrologic
Table III. The search results of different methods
Year Item M-5 rule GA SA GA– SA
1988 ˛ 0Ð46 0Ð85 0Ð83 0Ð86 GSI 193Ð63 5Ð68 9Ð24 4Ð81 1990 ˛ 0 0Ð94 0Ð96 0Ð95 GSI 2Ð45 0Ð47 0Ð78 0Ð20 1993 ˛ 0 0Ð63 0Ð66 0Ð73 GSI 721Ð26 32Ð12 32Ð79 29Ð12
conditions. The GA has a stepped response in the first 100 generations and then reaches a steady state. The SA shows variability at the outset then climbs to a steady state. Figure 6 presents the searching performance
through the GA– SA interaction stages. The ordinate is the degree of satisfaction. The upper abscissa indicates the number of temperature decreasing stages by SA and the lower abscissa represents the number of generations by the GA. After five stages of SA and less than 60 gener-ations of the GA, the objective values reach a stable (con-verged) condition in all three cases. The hybrid GA–SA gives interesting and promising results, where the objec-tive value is consistently increased in each GA–SA inter-acting stage. It appears that, in the first two or three GA–SA interacting stages, the GA contributes more than SA in terms of increasing the objective values, whereas the SA dominates the improvement of system perfor-mance in the later interacting stages in 1988 and 1990. 0.86 0.84 0.82 0.8 0.78 0.76 0.74 objective value a 0 50 100 150 200 250 300 350 400 0 50 100 150 200 250 300 350 400 number of generation (a1) GA in 1988 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 objective value a 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 number of generation (a2) SA in 1988 0 2 4 6 8 10 12 14 number of generation (b2) SA in 1990 0.94 0.93 0.92 0.91 0.9 0.89 0.88 0.87 0.86 0.85 objective value a number of generation (b1) GA in 1990 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 objective value a 0.7 0.6 0.5 0.4 0.3 0.2 0.1 objective value a 0 0 0.5 1 1.5 2 2.5 3 3.5 number of generation (c2) SA in 1993 0.74 0.72 0.7 0.68 0.66 0.64 0.62 0.6 objective value a 0 50 100 150 200 250 300 350 400 number of generation (c1) GA in 1993 × 105 × 104 × 105
0.88 0.86 0.84 0.82 0.8 0.78 0.76 0.74 0.72 0.7 0 0 4 0 1 2 5 10 20 30 14 24 34 44 54 10 20 30 40 50 60 objectgive value a 0.95 0.94 0.93 0.92 0.91 0.9 0.89 0.88 0.87 0.86 objectgive value a 0.74 0.72 0.7 0.68 0.66 0.64 0.62 0.6 0.58 0.56 objectgive value a Number of generation (a) GA-SA in 1988 Number of generation (b) GA-SA in 1990 Number of generation (c) GA-SA in 1993 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 Number of Temp. decreasing stage
Number of Temp. decreasing stage
Number of Temp. decreasing stage
5th Temp. stage Initial stage 1st Temp. stage 2nd Temp. stage 3rd Temp. stage 4th Temp. stage initial stage 1st Temp. stage 2nd Temp. stage 3rd Temp. stage 4th Temp. stage 5th Temp. stage 2nd Temp. stage 3rd Temp. stage 4th Temp. stage 5th Temp. stage initial stage 1st Temp. stage
Figure 6. The performance of the GA–SA through the searching pro-cesses in (a) 1988, (b) 1990 and (c) 1993
The number of searches (simulation times) is listed in Table IV. The searching times of the three methods are roughly estimated as follows.
1. GA: number of generations ð population size, i.e. 400 ð 300 in all three cases.
2. SA: number of iterations until a stable solution is reached which shows in the abscissa of part B Figure 5(a2), 5(b2), and 5(c2) for 1988, 1990, and 1993, respectively.
Table IV. The searching times of three methods in different years
Year Number of searches
GA SA GA– SA
1988 12 ð 104 44.9 ð 104 16.8 ð 104
1990 12 ð 104 13.6 ð 104 16.5 ð 104
1993 12 ð 104 32.4 ð 104 16.1 ð 104
3. GA–SA: number of generations ð population size + number of temperature decreasing stages ð number of iterations in each stage ð population size (e.g. 60 ð 300 C 5 ð 100 ð 300 D 16.8 ð 104 for 1988). Table IV shows the searching times of the three methods in different years. The searching times of the hybrid GA– SA are slightly greater than the GA (the total computation time of the GA–SA is mainly contributed by SA), whereas SA usually takes more time to reach a stable solution. It is worth mentioning that the searching times of the hybrid GA– SA can be dramatically reduced if the population size is reduced. For instance, if one executes an SA search based on those of better solutions obtained from the GA process, e.g. upper 100 solutions in population size of 300, one might obtain almost the same final objective values, whereas the number of searches could be reduced to 68 000 (60 ð 300 C 5 ð 100 ð 100 D 6.8 ð 104) times for 1988.
To testing and evaluate the hybrid GA–SA method further, we implemented a search on the long-term reservoir operation performance. Consecutive 10-year monthly historical inflow data were used. The search procedure, as above, is performed annually. Each year, 12 optimizing monthly releases and water storage levels were searched, and the last water storage level of the year was then propagated into the next year as its initial water storage level. For the purpose of comparison, the M-5 rule curves were also performed based on same inflow conditions and water requirements (as mentioned above). The results are presented in Table V and Figure 7. It appears that the hybrid GA– SA method outperforms the M-5 rule curves in terms of a much higher degree of satisfaction and lower GSI values in all 10 years. From Figure 7, we can easily see that the M-5 simulation method would produce a very high GSI value if the previous year was a drought year. This is mainly because the previous year, if it is a drought, will have a low water storage level at the end of that year and the low water level will propagate to the current year for the first four or five drought months. This would certainly exaggerate the drought phenomenon and greatly increase the GSI values.
CONCLUSIONS
In this study, Shihmen Reservoir operation was first trans-formed by fuzzy programming to identify the optimal
Table V. The results of M-5 rules and GA–SA in 10 years 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 Inflow (106 m3) 1036 1488 2256 1138 1910 738 1885 885 1456 1656 GSI M-5 rules 194 432 257 81Ð0 143 9Ð84 487 7Ð78 255 186 GA–SA 4Ð81 1Ð80 0Ð60 8Ð19 0Ð05 35Ð8 6Ð00 22Ð7 4Ð55 6Ð87 ˛ M-5 rules 0Ð76 0 0 0Ð8 0 0Ð54 0 0Ð69 0 0 GA–SA 0Ð86 0Ð86 0Ð94 0Ð78 0Ð95 0Ð64 0Ð82 0Ð79 0Ð69 0Ð85 0 100 200 300 400 500 600 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 GSI 0 500 1000 1500 2000 2500 Inflow (million m 3) GSI of M-5 rules GSI of GA-SA Inflow
Figure 7. The inflow series and GSI values by the M-5 rules and the GA–SA method in 10 years
degree of satisfaction among the objectives and con-straints. Three natural-based evolutional search methods (GA, SA and a hybrid GA– SA) were used to search this complex and non-linear system to determine the optimal schedule of water release and storage for three years with different hydrological conditions. Based on our compu-tation results, the following conclusions can be drawn. 1. The degree of satisfaction proposed in fuzzy
program-ming is a suitable criterion to evaluate the reservoir operation with regard to water deficit and dam safety. 2. All three natural-based search methods have a superior performance, with regard to the high degree of satis-faction and lower GSI index, than the current M-5 rule curves approach in all cases.
3. The GA is more efficient than SA, whereas SA might have a better performance than the GA after a long search iteration. The proposed hybrid GA–SA method achieves the highest satisfaction and lowest GSI values in all cases. The hybrid GA– SA can avoid a long and ineffective searching process in SA at the beginning stages and improve the solution quality of the GA. Further evaluation of the hybrid GA-SA method was undertaken by searching the long-term reservoir operation performance; we implemented the method using a con-secutive 10-year period and compared the performance with the M-5 rule curves simulation. The results again demonstrate that the hybrid GA–SA method indeed out-performs the M-5 operating rule curves in terms of much higher degree of satisfaction and lower GSI values. These
results provide further evidence to support that the hybrid GA–SA method is a robust and efficient way for opti-mizing variables in a highly non-linear problem, such as water resources management.
ACKNOWLEDGEMENTS
This paper is based on partial work supported by National Science Council, ROC (grant no. NSC91-2313-B-002-315).
REFERENCES
Archibald TW, McKinnon KIM, Thomas LC. 1996. An aggregate stochastic dynamic programming model of multireservoir systems.
Water Resources Research 33: 333– 340.
Bellman R. 1957. Dynamic Programming. Princeton University Press: Princeton, NJ.
Bellman RE, Zadeh LA. 1970. Decision making in a fuzzy environmen-tal. Management Science B 17: 141– 164.
Chang FJ, Chen L. 1991. Investigation of the objection functions for real-time reservoir operation. Taiwan Water Conservancy Quarterly 39(2): 50– 58.
Chang FJ, Chen L. 1998. Real-coded genetic algorithm for rule-based flood control reservoir management. Water Resources Management 12: 185– 198.
Chang FJ, Hui SC, Chen YC. 2002. Reservoir operation using grey fuzzy stochastic dynamic programming. Hydrological Processes 16: 2395– 2408.
Chang FJ, Chen L, Chang LC. 2005. Optimizing the reservoir operation rule curves by genetic algorithms. Hydrological Processes 19: 2277– 2289.
Chang LC, Chang FJ. 2001. Intelligent control for modelling of real-time reservoir operation. Hydrological Processes 15: 1621– 1634.
Chaves P, Kojiri T, Yamashiki Y. 2003. Optimization of storage reservoir considering water quantity and quality. Hydrological Processes 17: 2769– 2793.
Datta B, Houck MH. 1984. A stochastic optimization model for real-time operation of reservoirs using uncertain forecasts. Water Research 20: 1039– 1064.
Dougherty DE, Marryott RA. 1991. Optimal groundwater management -1. Simulated annealing. Water Resources Research 27: 2493– 2508. Dubrovin T, Jolma A, Turunen E. 2002. Fuzzy model for real-time
reservoir operation. Journal of Water Resources Planning and
Management 128: 66– 73.
Eshelman LJ, Schaffer JD. 1993. Real-coded genetic algorithms and interval-schemata. In Foundations of Genetic Algorithms 2 , Whitley LD (ed.). Morgan Kaufmann: San Mateo, CA; 187– 202. Goldberg DE. 1989. Genetic Algorithms in Search, Optimization, and
Machine Learning, 2nd edn. Addison-Wesley: Reading, MA.
Hall WA, Buras N. 1961. The dynamic programming approach to water resources development. Journal of Geophysical Research 66: 517– 520. Holland JH. 1992. Adaptation in Natural and Artificial Systems. MIT
Press: Cambridge, MA.
Hsu SK. 1995. Shortage indices for water-resources planning in Taiwan.
Journal of Water Resources Planning and Management 121: 119– 131.
Kirpatrick S, Gelatt CD, Vecchi MP. 1983. Optimization by simulated annealing. Science 220(4598): 671– 680.
Lin FT, Kao CY, Hsu CC. 1993. Applying the Genetic approach to simulated annealing in solving some NP-hard problems. IEEE
Transactions on Systems, Man, and Cybernetics 23: 1752– 1767.
Maskey S, Guinot V, Price RK. 2004. Treatment of precipitation uncertainty in rainfall– runoff modelling: a fuzzy set approach.
Advances in Water Resources 27: 889– 898.
Mousavi SJ, Mahdizadeh K, Afshar A. 2004. A stochastic dynamic programming model with fuzzy storage states for reservoir operations.
Advances in Water Resources 27: 1105– 1110.
Oliveira R, Loucks DP. 1997. Operating rules for multireservoir systems.
Water Resources Research 33: 839– 852.
Pardo-Ig´uzquiza E. 1998. Optimal selection of number and location of rainfall gauges for areal rainfall estimation using geostatistics and simulated annealing. Journal of Hydrology 210: 206– 220.
Trezos T, Yeh WW-G. 1987. Use of stochastic dynamic programming for reservoir management. Water Resources Research 23: 983– 996. Wang DW. 1997. An inexact approach for linear programming problems
with fuzzy objective and resources. Fuzzy Sets and Systems 89: 61– 68. Wang PP, Zheng C. 1998. An efficient approach for successively perturbed groundwater models. Advances in Water Resources 21: 499– 508.
Wardlaw R, Sharif M. 1999. Evaluation of genetic algorithms for optimal reservoir system operation. Journal of Water Resources Planning and
Management 125: 25– 33.
Yu PS, Chen CJ, Chen SJ. 2000. Application of gray and fuzzy methods for rainfall forecasting. Journal of Hydrologic Engineering 5: 339– 345. Zadeh LA. 1965. Fuzzy sets. Information and Control 8: 338– 353. Zimmermann HJ. 1976. Description and optimization of fuzzy systems.
