Applying Simulated Annealing For Optimal Operation of Multi-Reservoir Systems

The need to optimize the operation of water reservoirs is an issue that is becoming increasingly a concern for water resources planners in developing countries. This issue particularly becomes more significant in large systems with multiple reservoirs where operation of one reservoir has an impact on the others. In other word, the set of reservoirs in these systems act like a united series and require specific methods to handle various modeling issues. Problem statement: Sirvan River basin in west of Iran, standing in fifth order in respect of discharge, is an example of such a complex system. The project of water transfer from western tropical regions which is one of the large-scale projects in water resources management in Middle East, consists of a number of reservoirs and transfer systems. The matter of optimum operation of such a collection is one of the complicated and outstanding issues in water resources management. Approach: It was found that, due to increasing decision-making variables, conventional models used in optimizing water resources systems were not any longer capable of obtaining a desired solution, either because of low precision or time constraints. Hence, the intelligent random research approach "Simulated Annealing" has been used which in recent, decades has revealed appropriate results in solving major problems. Results: The results of this research indicate that the annealing approach is capable of solving such complex problems in water resources management with good precision in a reasonable period of time. Conclusions/Recommendations: The results were also compared with outputs of MODSIM which is a widely known model for solving complex water resources systems problems and benefits from “out of kilter” algorithm. The results indicate SA as a very robust and effective model in optimization of large real multi-reservoir systems.


INTRODUCTION
For a multi reservoir system, where the number of reservoirs is large, the conventional modeling by Dynamic Programming (DP) or classical Stochastic Dynamic Programming (SDP) presents difficulty, due to the curse of dimensionality inherent in the model solution. It takes a long time to obtain a steady state policy; also, it requires large amount of computer storage space, which form drawbacks in application. So an attempt is made to explore the concept of "local search" to provide a viable alternative in this context. Simulated Annealing (SA), proposed by Kirkpatrick et al., is a randomized search method for optimization. It tries to improve a solution by walking randomly in the space of possible solutions and gradually adjusting a parameter called "temperature".
At high temperature, the random walk is almost unbiased and it converges to essentially the uniform distribution over the whole space of solutions. As the temperature drops, each step of the random walk is more likely to move towards solutions with a better objective value and the distribution is more and more biased towards the optimal solutions [1] . The sequence of temperatures and lengths of time for which they are maintained is called the annealing schedule in analogy with statistical mechanics.
The simulated annealing algorithm, though by itself it is a local search algorithm, avoids getting trapped in a local minimum by also accepting cost increasing neighbors with some probability. In SA, first an initial solution is randomly generated and a neighbor is found and is accepted with a probability between 0 and 1 [2] . Although "SA" can be implemented quite easily with the degree of coding quite minimal relative to other nonlinear optimization algorithms; but it is not without its critics. One negative feature of SA is that it can be quite time-consuming to find an optimal fit, especially when using the "standard" Boltzmann technique.
Cunha and Sousa [3] utilized this method to obtain the least-cost design of a looped water distribution network. The implemented simulated annealing for optimum installation of sampling sites in the river. Cunha used for groundwater management. Teegavarapu and Simonovic [4] used that for the best exploitation from four-reservoir system with irrigation and hydroelectro purposes. In order to utilized this approach to prepare a program for distribution of water in irrigating channels. By applying this method to prepare a mathematics model of optimizing the hydraulic performance of irrigating channels.
In this study we demonstrate how simulated annealing, a Heuristic algorithm, can be used to solve high dimensional, linear and non-linear optimization problems for multi-reservoirs allocation system. The main objective pursued in the research is to obtain the optimal operating policies by minimizing the total irrigation deficits during a critical drought year in Sirvan river basin.
It should be noted that all of the computations in this study accomplished by a personal computer with Pentium IV, 2800 MHz processor and 256 MB RAM and computation times are corresponding to this computer.

MATERIALS AND METHODS
In order to choose a proper approach to model the main problem (herein Sirvan river basin), the computation was initially performed over a single reservoir model. For this purpose Dez reservoir was considered, the highest double arched concrete reservoir in southwest of Iran with the total capacity of 3460 MCM.
Dez reservoir is one of the first multi objective reservoirs in Iran. The main goals of this structure and relative establishments are as follows: • Controlling the vernal flood water and prevention of damages of flood • Providing the water for irrigation and agriculture uses • Producing hydro-power energy and supplying the water for massive industry of state of Khuzestan To better illustrate the efficiency of various methods, two models have been defined: "short-term and long term model". The Short-term model was performed only for one year by means of the average inflows over a 42 year history record. On the other hand, the long-term model was computed for all the 42year historic records. Either of models was considered for both instances that the final storage of water in the reservoir is equal to its initial storage at the beginning of the period and also unequal to its initial storage.
For simplicity, water supply is considered as the only operating aim. Evaporation and direct rainfall on the reservoir also are ignored and monthly demands are fixed. All of the simulations begin in October 1956-September 1998. Figure 1 shows Annual Reservoir inflows and related moving averages. The objective of these models is to minimize the annual sum of squared deviations from target demands. In Eq. 1, belongs to short-term model, TD t is the total of municipal and agricultural Demands needed in each particular period t and R t represents the allocated flow from the reservoir in period t: These optimization methods are based on massbalance equations for routing flows through the reservoir. The mass balance or continuity equations which have been shown in Eq. 2, explicitly define storage volume at the beginning of each period t.
where, S t , Q t and R t are the storage volume in the reservoir, the inflow to and the release from the reservoir at the beginning of period t, respectively. It was assumed that the initial and final reservoir storages in the whole period are unknown but equal. However, In the second scenario this constraint was relaxed.
The objective function of long-term model, similar to the short-term model, is to minimize the total longterm deficits as explained by Eq. 3: The model was developed on monthly basis for operation. In the model an initial storage was considered as initial volume in the first period and then amount of release in each period was selected according to request approach. Then with regard to mass balance equation shown in Eq. 4 the final storage of each period was computed.
The performance of Simulated Annealing is also compared with the Genetic Algorithm (GA), Ant Colony Optimization (ACO) and DP programming. After Dez-model the proposed algorithm was applied to a complex system with 22 reservoirs and several water demand sections. In all, the system contains 52 nodes which are basically points between links. This resource allocation system requires input from extensive spatial databases and involves a complex decision making problem. Figure 2 shows the schematic topology of this system. The availability of high-speed digital computers of large capacity makes it possible to simulate the performance of relatively complex river-basin systems for periods of any desired length, whereas conventional methods using only desk computers permit operating only selected parts of the system during a limited period of hydrologic record.
It is obvious from Fig. 2 that an interconnected reservoir system can be represented as a network of nodes and links. Nodes depict storage or non storage points of confluence or diversion and links represent reservoir release, river or channel flows. So the attention must be paid to operational effectiveness and the efficiency of existing reservoir system for maximizing the beneficial uses of this project. The statistics and information that have been used are as follows: 46 years statistics of historic yields, minimum and maximum storage volumes in different months, total water demand in different months and surface-volume-height curve data.
Like previous model, here, we have also considered two different states. In the state I, the optimum operating policies of reservoirs are obtained by one by one reservoirs' optimization, but in state II, the system is considered as an integrated model.
As the performance and the accuracy of "simulated annealing" has introduced it as an elite method, in this large scale model this approach has been considered and also the results of that have been compared with the result of MODSIM model.
MODSIM is a general-purpose simulationoptimization model, incorporated into a Decision Support System (DSS). It was originally developed by Dr. John Labadie of Colorado State University (CSU) in the mid-1970's and is currently being used by the U.S. Bureau of Reclamation for operational planning in the Upper Snake River Basin. The Windows-based Graphical User Interface (GUI) in MODSIM allows the user to create any reservoir system topology by simply clicking on various icons and placing system objects in any desired configuration on the screen. Data structures embodied in each model object on the screen are controlled by a database management system, with formatted data files prepared interactively and a network flow optimization model automatically executed from the interface. Results of the optimization are presented in useful graphical plots, or even customized reports available through a scripting language included with MODSIM complex, nonnetwork constraints on the optimization in MODSIM are incorporated through an iterative procedure using the embedded PERL scripting language [5] .
For adjusting the algorithm parameters of SA, we should specify the Initial Temperature (T 0 ), the Final Value of Temperature, the Length of Markov Chains (epoch length), the Temperature Schedule and the Acceptance Criteria. Also we must declare the randomly tweak solution.
The initial value of T, T 0 , is determined in such a way that virtually all transitions are accepted, i.e. T 0 is such that exp (-∆ E/T 0 ) ≈ 1 for almost all i and j.
Skiscim and Golden [6] and Lundy and Mees [7] determine T 0 by calculating the average increase in cost, where, x 0 or acceptance ratio defined as a number of accepted transitions divided by a number of proposed transitions and E represents the energy or objective function. Equation 6 leads to the following choice for T 0 : In current optimization models x 0 = 0.9. A stop criterion has been considered by terminating execution of the algorithm if the last configurations of consecutive Markov chains are identical for a number of chains.
The simplest choice for L k , the length of k-th Markov chain, is a value depending on the size of the problem. Thus, L k is independent of k. Here L k has been selected among some various values in sensitivity analysis test. Sensitivity analysis is a technique that allows the identification of the parameters and variables which have a greater influence on system performance and allows evaluation of the scale of this influence.
First we chose a geometric schedule using the equation below: where, α is a constant smaller than but close to 1and [.8 .99] α ∈ − . In SA algorithm, we have two solutions. The first is our original solution called the current solution and the second is the tweaked version called the working solution. Each has an associated energy which is the strength of the solution (let's say that the lower the energy, the better the solution).Our working solution is then compared to the current solution. If the working solution has less energy than the current solution (i.e., is better solution), then we copy the working solution to the current solution and move on to temperature reduction [8] .
However, if the working solution is worse than the current solution, we evaluate the acceptance criteria to figure out what to do with the current working solution. The probability of acceptance is based on Eq. 6 (which is based upon the law of thermodynamics) [1] : To randomly modify the working solution, the following equation is used: Gam(a) + = + β (10) Let i k a be a parameter in dimension i generated at annealing-time k, β be a random variable [ 1,1] β ∈ − and Gam(a) be the pace of motion in the neighborhood space of possible solutions and is dependent to the size of problem and is selected through a sensitivity analysis of various values of this parameter.

RESULTS
The results of Dez short-term model are tabulated in Table 1. The quality of the solution obtained by SA algorithm is quite much better than the other methods. Because the parameters and the resulted objective function may vary in different runs of the model, here average results of 10 runs are considered for evaluation and comparison purposes. Figure 3 shows the values of objective function which have been made by accepted transitions.   Table 3 makes a comparison between the results of SA via long-term model and those of different methods. The comparison shows that simulated annealing provides optimal solutions, converging into the fitness values within a short time. The SA is able to produce more satisfactory results as compared to those by GA or DP in much shorter run times.
There are still other procedures that facilitate the SA and hence making it possible to obtain similar results in less time. Furthermore in course of the algorithm's execution, the pace of motion values is gradually lowered, eventually to 0, in which case only improvements are accepted. This change made the results much better, whereas using logarithmic schedules [9] , shown in (11), had no more effective improvement: The results of above tables confirm the advantage of SA model in comparison with other methods both in terms of objective values and computational times. In other words, not only the objective value of SA is less than by an amount or at least equal to the objective value of DP or GA, but also the running time of this method is less than the others. Figure 4 and 5 show the comparison between SA and GA and SA using different logarithmic temperature schedules, respectively.
The first model of Sirvan large scale system has been executed for initial theoretical release values of 0, 0.7, 0.8, 0.9 and 1.0 times each node's demands. Fig. 4: Comparison between GA and SA is given for function f 0 ; the abscissa indicates the number of function calls, while the ordinate shows the best function evaluation found so far. Although the obtained optimal solutions are the same, but the run time of GA is nearly 16 times more than SA      Table 4. The best result from different initial theoretical release values has been obtained for initial release of 0.7 times demand. Thus, accrued results from this model were compared with MODSIM model. In order to assimilate the results of two models, the values of losses caused by deviations from target operating conditions in MODSIM model were being to the power of 2 and compared with SA model. Accrued results are shown in Table 5. In Fig. 6 quantitative reliability of both models have been compared. As it is implied, in particular, in downstream nodes of Sirvan, MODSIM model was better capable to provide region demands than simulated annealing model. Of course it is due to the linear objective function of MODSIM model which optimize the linear deviations from desired demands. Despite the previous model, in this model system is considered to be integrated. The results of this model were also compared with those of MODSIM model. In Table 6 the results are shown. As results indicate, in this situation, annealing model has done better than MODSIM model, in both supplying the downstream nodes' demands and obtaining a good result for both linear and squared deviations from target demands.    Fig. 11: The average of optimal storage volume in lower SIRVAN and MODSIM model. Figure 9 shows the comparison between maximum monthly shortages of SA and MODSIM model and finally Fig. 10 and 11 show the average of optimal storage volume for all 22 reservoirs.

DISCUSSION
Application of Simulated Annealing algorithm is evaluated in this study for optimization of simple to multi reservoir system. As the first problem, a singlereservoir system is considered and the results of SA and some heuristic methods have been compared. As far as the performance of SA is concerned, it can be inferred that for many applications the quality of the solution produced by the algorithm is much better than the other heuristic methods. Therefore SA can be considered as a new method in optimization of multi-reservoir systems, especially in real large system, even with nonlinear complex objective function.

CONCLUSION
Optimization and exploitation management of a multi-reservoir system is such an important issue that must be considered by planning staff. It well can be elicited that small improvements in operation of these systems will result in significant profits. Obtained results from simulated annealing method show that to achieve the best answer, proper determination of the pace of motion in the possible answer space is of great importance and if these variations do not fall within an appropriate range, the obtained answer may be locally optimum. The demonstrated model in this study reveals the high potential of annealing method in solving the water resources problems and in particular the exploitation of reservoir. Hence with regard to easy coding and acceptable results of this approach, it can be selected and used as the first option in solving the optimization problems.