- Email: [email protected]
An innovative simulated annealing approach to the long-term hydroscheduling problem
Electrical Power and Energy Systems 25 (2003) 41±46
www.elsevier.com/locate/ijepes
An innovative simulated annealing approach to the long-term hydroscheduling problem A.H. Mantawy a, S.A. Soliman b,*, M.E. El-Hawary c a
Department of Electrical Power and Machines, Ain Shams University, Abbassia, Cairo, Egypt b Department of Electrical Engineering, University of Qatar, P.O. Box 2713, Doha, Qatar c Department of Electrical and Computer Engineering, Dalhousie University, P.O. Box 1000, Halifax, NS, Canada Received 21 May 2001; accepted 13 February 2002
Abstract This paper presents a new simulated annealing algorithm (SAA) to solve the long-term hydroscheduling problem. A new algorithm for randomly generating feasible trial solutions is introduced. The problem is a hard nonlinear optimization problem in continuous variables. An adaptive cooling schedule and a new method for variables discretization are implemented to enhance the speed and convergence of the original SAA. A signi®cant reduction in the number of the objective function evaluations, and consequently less iteration are required to reach the optimal solution. The proposed algorithm has been applied successfully to solve a system with four series cascaded reservoirs. Numerical results show an improvement in the solutions compared to previously obtained results. q 2002 Published by Elsevier Science Ltd. Keywords: Simulated annealing; Simulated annealing algorithm; Hydroelectric plants
1. Introduction In this paper the long-term hydroscheduling problem (LTHSP) of a multi-reservoir hydropower plant connected in series on a river is considered [1±17]. The LTHSP optimization problem involves the use of a limited resource over a period of time. The resource is the water available for hydro generation. Most hydroelectric plants are used as multipurpose plants. In such cases, it is necessary to meet certain obligations other than power generation. These may include a maximum fore bay elevation, not to be exceeded because of danger of ¯ooding, and a minimum plant discharge and spillage to meet irrigational and navigational commitments. Thus, the optimal operation of the hydro system depends on the conditions that exist over the optimization interval [9]. Other distinctions among power systems are the number of hydro stations, their location, and special unit operating characteristics. The problem is quite different if the hydro stations are located on the same stream or on different ones. An upstream station will highly in¯uence the operation of the next downstream station. The latter, however, also * Corresponding author. Tel.: 1974-479-4381; fax: 1974-485-2118. E-mail addresses: [email protected] (S.A. Soliman), [email protected] (M.E. El-Hawary). 0142-0615/03/$ - see front matter q 2002 Published by Elsevier Science Ltd. PII: S 0142- 061 5( 02) 00019-4
in¯uences the upstream plant by its effect on the tail water elevation and effective head. Close coupling of stations by such a phenomenon is a complicating factor [1±9]. The evaluation of the optimum monthly operating policy of a multi-reservoir hydroelectric power system is a stochastic nonlinear optimization problem with continuous variables. The problem of determining the optimal long-term operation of multi-reservoir power system has been the subject of numerous publications over the past 50 years, and yet no completely satisfactory solution has been obtained, since in every proposed method the problem is simpli®ed. Many classical methods [1±17] have been used to solve the LTHSP including; dynamic programming, stochastic, aggregation linear programming and decomposition methods. These methods need complex mathematical manipulations in addition to requiring excessive computing time and storage requirements. Simulated annealing (SA) is a powerful technique to solve many optimization problems [18±23]. It has the ability of escaping local minima by incorporating a probability function in accepting or rejecting new solutions. It was ®rst proposed in the area of combinatorial optimization, that is, when the cost function is de®ned in a discrete domain. Therefore, many important problems are de®ned as functions of continuous variables, as the problem under study. The
42
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
3. The water conservation equation for each reservoir is adequately described by the following continuity difference equations: xi;t xi;t21 1 Ij;t 2 ui;t 2 si;t 1 s i21;t
1
1 # i # N; 1 # t # K Fig. 1. The system under study.
application of SA method to such type of problems requires an ef®cient strategy to generate feasible trial solutions. The other important point is the selection of an ef®cient cooling schedule that adapt itself according to statistics of the trial moves (acceptance/rejection) during the search. In this paper, we propose a simulated annealing algorithm (SAA) to solve the LTHSP of a multi-reservoir hydropower plant connected in series on a river. In the proposed algorithm an adaptive strategy to control the step sizes of varying the variables between coarse and ®ne, according to the accepted and rejected trials. Moreover, a polynomial time cooling schedule is used to control the temperature scheme based on the statistics of trial solutions generated during the search process [18,23]. A major contribution of this work is to introduce new rules for randomly generating feasible solutions satisfying all the system constraints including the well-known continuity equation of all reservoirs at all time periods. An example system includes four series reservoirs, is taken from literatures [2,9]. Comparing the results obtained with those of other classical optimization methods show an improvement in the quality of solutions obtained with the proposed SAA. In Section 2, the notation will be given along with a mathematical formulation of the problem. In Section 3, a general description of the SA method is presented, followed, in Section 4, by the general steps of the proposed SAA to solve the LTHSP. Section 5 presents the new proposed rules for generating feasible trial solutions. Section 6 is intended to offer the detailed description of the proposed SAA. In Section 7 the computational results along with a comparison with previously published work are presented. Section 8 outlines the conclusions. 2. Problem statement The long-term optimization problem of the hydroelectric power system shown in Fig. 1 is to determine the discharge, ui,t, the storage, xi,t, and the spillage, si,t for each reservoir, i 1; ¼; n; at all scheduling time periods, t 1; ¼; K; under the following conditions [9]: 1. The expected value of the water left in storage at the end of the last period studied is a maximum. 2. The expected value of the MWH generated during the optimization interval is a maximum.
where xi,t is the storage of reservoir i at the end of period t, Ii,t the in¯ow to the reservoir i during the period t and si,t is the spillage from reservoir i during the period t. Spillage occurs when the discharge exceeds the maximum discharge value and the storage exceeds the maximum allowed storage of the reservoir. 4. To satisfy the multipurpose stream use requirements, the following operational constraints should be satis®ed: max xmin i; # xi;t # xi; ; max umin i;t # ui;t # ui;t ;
1 # i # N; 1 # t # K
2
1 # i # N; 1 # t # K
3
Where xmax is the capacity of the reservoir, xmin i; i; the minimum storage, umin the minimum discharge through the i;t turbine i during the period t, and umax is the maximum i;t discharge through the turbine i during the period t. If ui;t $ umax and xi,t is equal to xmax then ui;t 2 umax is i;t i i;t discharged through the spillways. In mathematical terms, the long-term optimization problem for the power system shown in Fig. 1 is to determine the discharge ui,t and the storage xi,t (consequently the spillage si,t) that maximize " JE
N X i1
Vi Xi;t 1
N X K X i1 t1
# ct Gi Xi;t21 ; ui;t
4
Subject to satisfying the equality constraints given by Eq. (1) and the inequality constraints given by Eqs. (2) and (3). Vi(xi,t) is the value of water left in storage in reservoir i at the end of the last period studied, Gi(xi,t21,ui,t) the generation of plant i during period t in MWh and ct is the value (in dollars) of the water on the river in month t. As shown above the LTHSP is a nonlinear optimization problem with continuous variables. The problem size can be determined as the product of three times the number of reservoirs and the number of time periods in the required schedule. In the LTHSP under consideration, one is interested in a solution that maximizes the total power generated during the scheduling time horizon, Eq. (4), while several constraints are satis®ed. 3. SA method 3.1. Physical concepts of SA SA is a Monte Carlo technique for ®nding solutions to
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
optimization problems [18±23]. Annealing, physically, refers to the process of heating up a solid to a high temperature followed by slow cooling achieved by decreasing the temperature of the environment in steps. At each step the temperature is maintained constant for a period of time suf®cient for the solid to reach thermal equilibrium. At equilibrium, the solid could have many con®gurations, each corresponding to different spins of the electrons and to speci®c energy level. At equilibrium the probability of a given con®guration, Pconfg ; is given by Boltzman distribution; Pconfg K exp 2Econfg =Cp; where Econfg is the energy of the given con®guration and K is a constant [21]. Ref. [22] proposed a Monte Carlo method to simulate the process of reaching thermal equilibrium at a ®xed temperature Cp. In this method, a randomly generated perturbation of the current con®guration of the solid is applied so that a trial con®guration is obtained. Let Ec and Et denote the energy level of the current and trial con®gurations, respectively. If Ec . Et ; then a lower energy level has been reached, the trial con®guration is accepted, and becomes the current con®guration. On the other hand, if Ec # Et then the trial con®guration is accepted as the current con®guration with probability exp Ec 2 Et =Cp: The process continues where a transition to a con®guration of higher energy level is not necessarily rejected. Eventually thermal equilibrium is achieved after a large number of perturbations, where the probability distribution of a con®guration approaches Boltzman distribution. By gradually decreasing Cp and repeating Metropolis simulation, new lower energy levels become achievable. As Cp approaches zero, the least energy con®gurations will have a positive probability of occurring. 3.2. Application of SA to optimization problems In applying the SAA, to solve optimization problems, the basic idea is to choose a feasible solution at random and then get a neighbor to this solution. A move to this neighbor is performed if either it has a better (lower) objective value or, in case the neighbor has a higher objective function value, if exp 2DE=Cp $ R 0; 1; where DE is the increase in objective value if we move to the neighbor and R(0,1) is a random number between 0 and 1. The effect of decreasing Cp is that the probability of accepting an increase in the objective function value is decreased during the search. 4. The proposed SAA for the LTHSP Implementation details of the SAA to solve the LTHSP are given in Sections 5 and 6. The major steps of the algorithm are summarized as follows: Step (0): Set iteration counter K 0. Set the initial temperature of the cooling schedule that results in high probability of accepting new solutions.
43
Initialize all step size vectors (Ustepo(i,t), Xstepo(i,t), and Sstepo(i,t)) for the variables (U(i,t), X(i,t), and S(i,t)), respectively, for all units i 1; 2; ¼; N at all time periods t 1; 2; ¼; K: Step (1): Find, randomly, an initial feasible solution (see Section 6.1). Step (2): Calculate the objective function at the initial solution (Eq. (4)). Step (3): Generate an integer random number J IR 1; 3: If J 1 then generate the trial solution by randomly varying the reservoir discharge value. If J 2 then generate the trial solution by randomly varying the reservoir storage value. If J 3 then generate the trial solution by randomly varying the reservoir spillage discharge value (see Section 5). Step (4): Calculate the objective function at the generated trial solution (Eq. (4)). Step (5): Perform the SA acceptance test; to accept or reject the trial solution (see Section 6.2). Step (6): Check for equilibrium at this temperature (see Section 6.5). If equilibrium reached go to step (7), Else go to step (3). Step (7): If the prespeci®ed maximum number of iterations is reached then stop, Else go to step (8). Step (8): If the step size vector values for all variables (Ustep, Xstep, Sstep) are less than a prespeci®ed value then stop, Else go to step (9). Step (9): Update the step size vector values (see Section 6.3). Decrease the temperature according to the polynomial time cooling schedule (see Section 6.4). Go to step (3). 5. Generating trial solution (neighbor) The most important part in the SAA is to have a good rule for ®nding a diversi®ed and intensi®ed neighborhood (trial solution) so that a large amount of the solution space can be explored. Neighbors should be randomly generated, feasible, and span as much as possible the problem solution space. Because of the dif®culty of satisfying the continuity equation (Eq. (1)) of all reservoirs at all time periods generating random trial feasible solution is not a simple matter. A major contribution of this work is the implementation of new rules to obtain randomly feasible solutions faster. Three different routines are implemented to generate trial solutions based on the random variation of one of the three variables (U, X, and S) at a time instant. The following steps describe the rules for generating trial solutions by randomly
44
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
varying the reservoir discharge (U): Step (1): Generate randomly a unit i IR 1; N; and a time instant t IR 1; K IR : integer random: Step (2): Generate R(0,1) Let DU R 0; 1 £ Ustep i; t; and U i; t U i; t 1 or 2 randomlyDU Step (3): Check for the feasibility of U(i,t): If Umin(i,t) # U(i,t) # Umax(i,t) Then go to step (4), Else go to step (1). Step (4): Use the continuity Eq. (1) to calculate the storage value of unit i at time t, X(i,t). Step (5): Check for the feasibility of X(i,t): If Xmin(i) # X(i,t) # Xmax(i) Then go to step (6), Else go to step (1). Step (6): Since the variation of U(i,t) will affect U, X and S at units i; i 1 1; ¼; N at periods t; t 1 1; ¼; K: The next step is then used to check the feasibility of all those values. Use the continuity equation to calculate X(i1,t1) for i1 i; i 1 1; ¼; N and t1 t; t 1 1; ¼; K At each value perform the following substeps: (6a) If Xmin(i1) # X(i1,t1) # Xmax(i1) Then go to step (7), Else go to (6b) (6b) If X(i1,t1) $ Xmax(i1) Then Let Utemp U(i1,t1) 1 X(i1,t1) 2 Xmax(i1) and X(i1,t1) Xmax(i1) Check If Utemp # Umax(i1,t1) Then U(i1,t1) Utemp Else U(i1,t1) Umax(i1,t1) and S(i1,t1) S(i1,t1) 1 Utemp 2 Umax(i1,t1) (6c) If X(i1,t1) # Xmin(i1) Then Let Utemp U(i1,t1) 1 X(i1,t1) 2 Xmin(i1) and X i1; t1 Xmin i1 Check If Utemp $ Umin(i1,t1) Then U(i1,t1) Utemp, Else U(i1,t1) Umin(i1,t1) and S(i1,t1) S(i1,t1) 1 Utemp 2 Umin(i1,t1) Step (7): If all reservoir and time periods have been treated then go to step (8), otherwise go to step (6a) Step (8): Check for the continuity equation at all units and all time periods. If satis®ed, then stop; otherwise go to step (1). 6. Details of the SAA for the LTHSP 6.1. Generating an initial feasible solution The proposed SAA requires a starting feasible solution, which satis®es all system constraints. This solution is randomly generated. The following steps are used in ®nding this starting solution: Step (1): (1a) Let i 1 and t 1
(1b) Let U i; t Umax i; t; X i; t Xmax i and S i; t 0 Step (2): Check if S(i,t) $ 0, then go to Step (4), Else let U1 Uin £ Umax i; t and X1 Xin £ Xmax i: Where Uin and Xin are a prespeci®ed parameters for ®nding initial feasible starting values for the variables (assumed 0.25) Step (3): U i; t R U1; Umax i; t X i; t R X1; Xmax i Go to step (2) Step (4): Check: If i N and t K then go to step (5), Else increase i and t and go to Step (1b). Step (5): Check the continuity equation for all units at all time periods. If satis®ed, then stop; Else go to step (1). 6.2. SA test The implementation steps of the SA test, as applied to each iteration in the algorithm, are described as follows [18]: Step (1): At the same calculated temperature, ckp ; apply the following acceptance test for the new trial solution. Step (2): Acceptance test: If Ej 2 Ei ; or If exp Ei 2 Ej =Cp $ R 0; 1; then accept the trial solution, set Xi Xj and Ei Ej : Otherwise reject the trial solution. Where Xi ; Xj ; Ei ; Ej is the SA current solution, the trial solution and their corresponding cost, respectively. Step (3): Go to the next step in the algorithm. 6.3. Step size vector adjustment In this work the step vector is updated jointly with temperature, according to the acceptance rate of the attempted moves at the previous temperature stage [23]. All its components are updates simultaneously. It must be mentioned here that a good step vector adjustment scheme must not yield steps that are too large with too many uncorrelated and unaccepted moves, whereas steps that are too small result in an incomplete exploration of the variation domain, with small and frequent increasing of the objective function. The following explains how the step vector is adjusted mathematically: Step (1): Calculate P i NTRACP i=NTR i; i 1; ¼; N: Where NTRACP(i) is the number of trial accepted when the variable i is changed. NTR(i) number of trials attempted by changing variable i. Step (2): If P(i) . PMAX, then STEP i STEP i £ STEPMAX
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46 Table 1 Comparison with DEC-DP and M-NORM
DEC-DP [2,9] M-NORM [9] Our SAA % Increase of SAA
Table 3 The storage (mm 3)
Total pro®t ($)
VWE ($)
28,145,330 28,227,174 28,799,568 2.32
10,403,170 10,414,598 10,307,287 22.02
If P(i) , PMIN, then STEP i STEP i £ STEPMIN Where PMAX, PMIN, STEPMAX and STEPMIN are parameters taken in our implementation as 0.2, 0.8, 0.8 and 1.2, respectively [23]. 6.4. Cooling schedule A ®nite-time implementation of the SA algorithm can be realized by generating homogenous Markov chains of ®nite length for a ®nite sequence of descending values of the control parameter. To achieve this, one must specify a set of parameters that governs the convergence of the algorithm. These parameters form a cooling schedule. The parameters of the cooling schedules are: an initial value of the control parameter decrement function for decreasing the control parameter and a ®nal value of the control parameter speci®ed by the stopping criterion, and a ®nite length of each homogenous Markov chain. In this work a polynomial-time cooling schedules is used in which the temperature is decreased based on the statistics of the trial solutions acceptance or rejection during the search. Details of the implemented cooling schedule are described in details in Ref. [18]. 6.5. Equilibrium test The sequence of trial solutions generated in the SAA at a ®xed temperature is stopped as soon as thermodynamic equilibrium, detected by some adequate condition, is Table 2 The discharge (mm 3) Month
1 2 3 4 5 6 7 8 9 10 11 12
45
Reservoir number 1
2
3
4
43.87 1037 1071.08 1068.41 968 1064.36 1037 0 0 0 0.13 0
411.97 1418 1290.33 1464.95 1322.77 1160.65 1105 127 510.31 314.82 491.05 363.09
572.7 1500.39 1330.8 1546.35 1338.41 1178.01 1109.42 405.17 695.25 470.73 636.48 495.24
2298.13 2704.04 2143.63 2450.56 2854.51 3160.28 1363.13 711.3 1695.27 3030.41 2220.36 1948.48
Month
1 2 3 4 5 6 7 8 9 10 11 12
Reservoir number 1
2
3
4
7473.13 7265.13 6772.05 6097.64 5394.64 4563.28 3719.28 3938.28 5039.28 6926.28 8076.15 8900.15
569.9 519.9 524.65 274.11 14.35 0.06 0.06 0.05 103.75 569.93 570 569.91
49.27 48.88 49.4 0 2.36 0 2.37 4.06 0.12 49.21 49.79 49.63
3419.58 3416.92 3414.09 2996.52 1782.42 58.82 0 1179.7 3419.68 3420 3419.12 3419.89
reached. Then the temperature and step vectors are suitably adjusted [23]. The test of equilibrium is done as follows. If the NTRACP(T ) , N1 £ n and NTR(T ) , N2 £ n then continue at the same temperature, otherwise end of temperature stage. Where NTRACP(T ), NTR(T ) are the number of trial accepted and attempted at temperature T, respectively, n the number of variables in the problem, N1 and N2 are the end temperature stage parameters (taken as 12 and 100) [23].
7. Numerical examples A computer model has been implemented based on the SAA described earlier and using the new proposed rules for generating random trial solutions. To test the model, an example system with four reservoirs connected in series on a river is solved. The scheduling time horizon is 1 year divided into 12 periods. The full data of this example are given in Refs. [2,9]. This Example was solved by combined decomposition-dynamic programming (DEC-DP) method [2,9] and by the minimum norm (MNORM) approach [9]. Table 1 shows a comparison of the results obtained by the DEC-DP approach, the M-NORM method and our proposed SAA. As shown in the table, our SAA achieves a higher total pro®t with a percentage increase of 2.32% relative to the DEC-DP approach and a value of 2.02% relative to the MNORM method. The value of water left in storage in the reservoirs at the end of the year (VWE) is also given in $ for the 3 methods. Detailed results for the Example test system are given in Tables 2 and 3. Table 2 shows the monthly discharge (releases) from the four reservoirs. Table 3 gives monthly storage in the four reservoirs. The solution obtained requires zero spillage from all reservoirs during the 12 months
46
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
except for the ®rst month, which requires 0.5, 0.6, 0.7 and 0.3 for the four reservoirs, respectively. 8. Conclusions In this paper we consider the LTHSP. The problem is a hard nonlinear optimization problem in continuous variables. The paper presents a new SAA for solving the LTHSP. The proposed algorithm introduces three new innovations; new rules for randomly generating feasible trial solutions, an adaptive cooling schedule and a method for variable discretization. This results in signi®cant reduction in the number of the objective function evaluations, and consequently less iterations are required to reach the optimal solution is obtained. The proposed SAA is applied successfully to solve a problem involving four series cascaded reservoirs hydroelectric power system. The results are compared with those obtained by DEC-DP [2,9] and MNORM [9]. Numerical results show a signi®cant increase of the total pro®t per year for the solved test system compared to previously obtained results using classical techniques. References [1] El-Hawary ME, Christensen GS. Optimal operation of electric power systems. New York: Academic Press, 1979. [2] Turgeon A. A decomposition method for the long-term scheduling of reservoirs in series. Water Resour Res 1981;17(6):1565±70. [3] Arvantidies NV, Rosing J. Composite representation of a multi reservoir hydroelectric power system. IEEE Trans PAS 1970; PAS-89(2):319±26. [4] Arvantidies NV, Rosing J. Optimal operation of multi reservoir systems using a composite representation. IEEE Trans PAS 1970; PAS-89(2):3327±35. [5] Turgeon A. Optimal operation of multi reservoir power system with stochastic in¯ows. Water Resour Res 1980;16(2):275±83. [6] Grygier JC, Stendinger J. Algorithms for optimization of hydropower system operation. Water Resour Res 1985;21(1):1±10. [7] Olcer S, Harsa C, Roch A. Application of linear and dynamic programming to the optimization of the production of hydroelectric power. Optimal Control Appl Meth 1985;6:43±56.
[8] Marino MA, Loaiciga HA. Dynamic model for multi reservoir operation. Water Resour Res 1985;21(5):619±30. [9] Christensen GS, Soliman SA. Optimal long-term operation of electric power systems. New York: Plenum Press, 1988. [10] Maceira MEP, Pereira MVF. Analytical modeling of chronological reservoir operation in probabilistic production costing. IEEE Trans Power Syst 1996;11(1):171±80. [11] Christoforidis M, et al. Long-term/mid-term resource optimization of a hydro-dominant power system using the interior point method. IEEE Trans Power Syst 1996;11(1):287±94. [12] Escudero LF, et al. Hydropower generation management under uncertainty via scenario analysis and parallel computation. IEEE Trans Power Syst 1996;11(2):683±9. [13] da Cruz Jr G, Soares S. Non-uniform composite representation of hydroelectric systems for long-term hydrothermal scheduling. IEEE Trans Power Syst 1996;11(2):702±7. [14] Wong KP, Wong SYW. Hybrid-genetic/simulated annealing approach to short-term multi-fuel constrained generation scheduling. IEEE Trans Power Syst 1997;12(2):776±84. [15] Guan X, et al. An optimization-based algorithm for scheduling hydrothermal power systems with cascaded reservoirs and discrete hydro constraints. IEEE Trans Power Syst 1997;12(4):1775±80. [16] Ponrajah RA, Galiana FD. Systems to optimize conversion ef®ciencies at Ontario hydro's hydroelectric plants. IEEE Trans Power Syst 1998;13(3):11044±50. [17] Yu Z, Sparrow FT, Bowen BH. A new long-term hydro production scheduling method for maximizing the pro®t of hydroelectric systems. IEEE Trans Power Syst 1998;13(1):66±71. [18] Mantawy AH, Abdel-Magid YL, Selim SZ. A simulated annealing algorithm for unit commitment. IEEE Trans Power Syst 1997; 13(1):197±204. [19] Aarts E, Korst J. Simulated annealing and Boltzman machines: a stochastic approach to combinatorial optimization and neural computing. New York: Wiley, 1989. [20] Cerny V. Thermodynamical approach to the traveling salesman problem: an ef®cient simulation algorithm. J Optim Theor Appl 1985;45(1). [21] Shokri Z, Selim K, Alsultan A. Simulated annealing algorithm for the clustering problem. Pattern Recog 1991;24(10):1003±8. [22] Metropolis N, Rosenbluth A, Rosenbluth M, Teller A, Teller E. Equations of state calculations by fast computing machines. J Chem Phys 1953;21:1087±982. [23] Siarry P, Berthiau G, Haussy J. Enhanced simulated annealing for globally minimizing functions of many continuous variables. ACM Trans Math Software 1997;23(2):209±28.
www.elsevier.com/locate/ijepes
An innovative simulated annealing approach to the long-term hydroscheduling problem A.H. Mantawy a, S.A. Soliman b,*, M.E. El-Hawary c a
Department of Electrical Power and Machines, Ain Shams University, Abbassia, Cairo, Egypt b Department of Electrical Engineering, University of Qatar, P.O. Box 2713, Doha, Qatar c Department of Electrical and Computer Engineering, Dalhousie University, P.O. Box 1000, Halifax, NS, Canada Received 21 May 2001; accepted 13 February 2002
Abstract This paper presents a new simulated annealing algorithm (SAA) to solve the long-term hydroscheduling problem. A new algorithm for randomly generating feasible trial solutions is introduced. The problem is a hard nonlinear optimization problem in continuous variables. An adaptive cooling schedule and a new method for variables discretization are implemented to enhance the speed and convergence of the original SAA. A signi®cant reduction in the number of the objective function evaluations, and consequently less iteration are required to reach the optimal solution. The proposed algorithm has been applied successfully to solve a system with four series cascaded reservoirs. Numerical results show an improvement in the solutions compared to previously obtained results. q 2002 Published by Elsevier Science Ltd. Keywords: Simulated annealing; Simulated annealing algorithm; Hydroelectric plants
1. Introduction In this paper the long-term hydroscheduling problem (LTHSP) of a multi-reservoir hydropower plant connected in series on a river is considered [1±17]. The LTHSP optimization problem involves the use of a limited resource over a period of time. The resource is the water available for hydro generation. Most hydroelectric plants are used as multipurpose plants. In such cases, it is necessary to meet certain obligations other than power generation. These may include a maximum fore bay elevation, not to be exceeded because of danger of ¯ooding, and a minimum plant discharge and spillage to meet irrigational and navigational commitments. Thus, the optimal operation of the hydro system depends on the conditions that exist over the optimization interval [9]. Other distinctions among power systems are the number of hydro stations, their location, and special unit operating characteristics. The problem is quite different if the hydro stations are located on the same stream or on different ones. An upstream station will highly in¯uence the operation of the next downstream station. The latter, however, also * Corresponding author. Tel.: 1974-479-4381; fax: 1974-485-2118. E-mail addresses: [email protected] (S.A. Soliman), [email protected] (M.E. El-Hawary). 0142-0615/03/$ - see front matter q 2002 Published by Elsevier Science Ltd. PII: S 0142- 061 5( 02) 00019-4
in¯uences the upstream plant by its effect on the tail water elevation and effective head. Close coupling of stations by such a phenomenon is a complicating factor [1±9]. The evaluation of the optimum monthly operating policy of a multi-reservoir hydroelectric power system is a stochastic nonlinear optimization problem with continuous variables. The problem of determining the optimal long-term operation of multi-reservoir power system has been the subject of numerous publications over the past 50 years, and yet no completely satisfactory solution has been obtained, since in every proposed method the problem is simpli®ed. Many classical methods [1±17] have been used to solve the LTHSP including; dynamic programming, stochastic, aggregation linear programming and decomposition methods. These methods need complex mathematical manipulations in addition to requiring excessive computing time and storage requirements. Simulated annealing (SA) is a powerful technique to solve many optimization problems [18±23]. It has the ability of escaping local minima by incorporating a probability function in accepting or rejecting new solutions. It was ®rst proposed in the area of combinatorial optimization, that is, when the cost function is de®ned in a discrete domain. Therefore, many important problems are de®ned as functions of continuous variables, as the problem under study. The
42
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
3. The water conservation equation for each reservoir is adequately described by the following continuity difference equations: xi;t xi;t21 1 Ij;t 2 ui;t 2 si;t 1 s i21;t
1
1 # i # N; 1 # t # K Fig. 1. The system under study.
application of SA method to such type of problems requires an ef®cient strategy to generate feasible trial solutions. The other important point is the selection of an ef®cient cooling schedule that adapt itself according to statistics of the trial moves (acceptance/rejection) during the search. In this paper, we propose a simulated annealing algorithm (SAA) to solve the LTHSP of a multi-reservoir hydropower plant connected in series on a river. In the proposed algorithm an adaptive strategy to control the step sizes of varying the variables between coarse and ®ne, according to the accepted and rejected trials. Moreover, a polynomial time cooling schedule is used to control the temperature scheme based on the statistics of trial solutions generated during the search process [18,23]. A major contribution of this work is to introduce new rules for randomly generating feasible solutions satisfying all the system constraints including the well-known continuity equation of all reservoirs at all time periods. An example system includes four series reservoirs, is taken from literatures [2,9]. Comparing the results obtained with those of other classical optimization methods show an improvement in the quality of solutions obtained with the proposed SAA. In Section 2, the notation will be given along with a mathematical formulation of the problem. In Section 3, a general description of the SA method is presented, followed, in Section 4, by the general steps of the proposed SAA to solve the LTHSP. Section 5 presents the new proposed rules for generating feasible trial solutions. Section 6 is intended to offer the detailed description of the proposed SAA. In Section 7 the computational results along with a comparison with previously published work are presented. Section 8 outlines the conclusions. 2. Problem statement The long-term optimization problem of the hydroelectric power system shown in Fig. 1 is to determine the discharge, ui,t, the storage, xi,t, and the spillage, si,t for each reservoir, i 1; ¼; n; at all scheduling time periods, t 1; ¼; K; under the following conditions [9]: 1. The expected value of the water left in storage at the end of the last period studied is a maximum. 2. The expected value of the MWH generated during the optimization interval is a maximum.
where xi,t is the storage of reservoir i at the end of period t, Ii,t the in¯ow to the reservoir i during the period t and si,t is the spillage from reservoir i during the period t. Spillage occurs when the discharge exceeds the maximum discharge value and the storage exceeds the maximum allowed storage of the reservoir. 4. To satisfy the multipurpose stream use requirements, the following operational constraints should be satis®ed: max xmin i; # xi;t # xi; ; max umin i;t # ui;t # ui;t ;
1 # i # N; 1 # t # K
2
1 # i # N; 1 # t # K
3
Where xmax is the capacity of the reservoir, xmin i; i; the minimum storage, umin the minimum discharge through the i;t turbine i during the period t, and umax is the maximum i;t discharge through the turbine i during the period t. If ui;t $ umax and xi,t is equal to xmax then ui;t 2 umax is i;t i i;t discharged through the spillways. In mathematical terms, the long-term optimization problem for the power system shown in Fig. 1 is to determine the discharge ui,t and the storage xi,t (consequently the spillage si,t) that maximize " JE
N X i1
Vi Xi;t 1
N X K X i1 t1
# ct Gi Xi;t21 ; ui;t
4
Subject to satisfying the equality constraints given by Eq. (1) and the inequality constraints given by Eqs. (2) and (3). Vi(xi,t) is the value of water left in storage in reservoir i at the end of the last period studied, Gi(xi,t21,ui,t) the generation of plant i during period t in MWh and ct is the value (in dollars) of the water on the river in month t. As shown above the LTHSP is a nonlinear optimization problem with continuous variables. The problem size can be determined as the product of three times the number of reservoirs and the number of time periods in the required schedule. In the LTHSP under consideration, one is interested in a solution that maximizes the total power generated during the scheduling time horizon, Eq. (4), while several constraints are satis®ed. 3. SA method 3.1. Physical concepts of SA SA is a Monte Carlo technique for ®nding solutions to
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
optimization problems [18±23]. Annealing, physically, refers to the process of heating up a solid to a high temperature followed by slow cooling achieved by decreasing the temperature of the environment in steps. At each step the temperature is maintained constant for a period of time suf®cient for the solid to reach thermal equilibrium. At equilibrium, the solid could have many con®gurations, each corresponding to different spins of the electrons and to speci®c energy level. At equilibrium the probability of a given con®guration, Pconfg ; is given by Boltzman distribution; Pconfg K exp 2Econfg =Cp; where Econfg is the energy of the given con®guration and K is a constant [21]. Ref. [22] proposed a Monte Carlo method to simulate the process of reaching thermal equilibrium at a ®xed temperature Cp. In this method, a randomly generated perturbation of the current con®guration of the solid is applied so that a trial con®guration is obtained. Let Ec and Et denote the energy level of the current and trial con®gurations, respectively. If Ec . Et ; then a lower energy level has been reached, the trial con®guration is accepted, and becomes the current con®guration. On the other hand, if Ec # Et then the trial con®guration is accepted as the current con®guration with probability exp Ec 2 Et =Cp: The process continues where a transition to a con®guration of higher energy level is not necessarily rejected. Eventually thermal equilibrium is achieved after a large number of perturbations, where the probability distribution of a con®guration approaches Boltzman distribution. By gradually decreasing Cp and repeating Metropolis simulation, new lower energy levels become achievable. As Cp approaches zero, the least energy con®gurations will have a positive probability of occurring. 3.2. Application of SA to optimization problems In applying the SAA, to solve optimization problems, the basic idea is to choose a feasible solution at random and then get a neighbor to this solution. A move to this neighbor is performed if either it has a better (lower) objective value or, in case the neighbor has a higher objective function value, if exp 2DE=Cp $ R 0; 1; where DE is the increase in objective value if we move to the neighbor and R(0,1) is a random number between 0 and 1. The effect of decreasing Cp is that the probability of accepting an increase in the objective function value is decreased during the search. 4. The proposed SAA for the LTHSP Implementation details of the SAA to solve the LTHSP are given in Sections 5 and 6. The major steps of the algorithm are summarized as follows: Step (0): Set iteration counter K 0. Set the initial temperature of the cooling schedule that results in high probability of accepting new solutions.
43
Initialize all step size vectors (Ustepo(i,t), Xstepo(i,t), and Sstepo(i,t)) for the variables (U(i,t), X(i,t), and S(i,t)), respectively, for all units i 1; 2; ¼; N at all time periods t 1; 2; ¼; K: Step (1): Find, randomly, an initial feasible solution (see Section 6.1). Step (2): Calculate the objective function at the initial solution (Eq. (4)). Step (3): Generate an integer random number J IR 1; 3: If J 1 then generate the trial solution by randomly varying the reservoir discharge value. If J 2 then generate the trial solution by randomly varying the reservoir storage value. If J 3 then generate the trial solution by randomly varying the reservoir spillage discharge value (see Section 5). Step (4): Calculate the objective function at the generated trial solution (Eq. (4)). Step (5): Perform the SA acceptance test; to accept or reject the trial solution (see Section 6.2). Step (6): Check for equilibrium at this temperature (see Section 6.5). If equilibrium reached go to step (7), Else go to step (3). Step (7): If the prespeci®ed maximum number of iterations is reached then stop, Else go to step (8). Step (8): If the step size vector values for all variables (Ustep, Xstep, Sstep) are less than a prespeci®ed value then stop, Else go to step (9). Step (9): Update the step size vector values (see Section 6.3). Decrease the temperature according to the polynomial time cooling schedule (see Section 6.4). Go to step (3). 5. Generating trial solution (neighbor) The most important part in the SAA is to have a good rule for ®nding a diversi®ed and intensi®ed neighborhood (trial solution) so that a large amount of the solution space can be explored. Neighbors should be randomly generated, feasible, and span as much as possible the problem solution space. Because of the dif®culty of satisfying the continuity equation (Eq. (1)) of all reservoirs at all time periods generating random trial feasible solution is not a simple matter. A major contribution of this work is the implementation of new rules to obtain randomly feasible solutions faster. Three different routines are implemented to generate trial solutions based on the random variation of one of the three variables (U, X, and S) at a time instant. The following steps describe the rules for generating trial solutions by randomly
44
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
varying the reservoir discharge (U): Step (1): Generate randomly a unit i IR 1; N; and a time instant t IR 1; K IR : integer random: Step (2): Generate R(0,1) Let DU R 0; 1 £ Ustep i; t; and U i; t U i; t 1 or 2 randomlyDU Step (3): Check for the feasibility of U(i,t): If Umin(i,t) # U(i,t) # Umax(i,t) Then go to step (4), Else go to step (1). Step (4): Use the continuity Eq. (1) to calculate the storage value of unit i at time t, X(i,t). Step (5): Check for the feasibility of X(i,t): If Xmin(i) # X(i,t) # Xmax(i) Then go to step (6), Else go to step (1). Step (6): Since the variation of U(i,t) will affect U, X and S at units i; i 1 1; ¼; N at periods t; t 1 1; ¼; K: The next step is then used to check the feasibility of all those values. Use the continuity equation to calculate X(i1,t1) for i1 i; i 1 1; ¼; N and t1 t; t 1 1; ¼; K At each value perform the following substeps: (6a) If Xmin(i1) # X(i1,t1) # Xmax(i1) Then go to step (7), Else go to (6b) (6b) If X(i1,t1) $ Xmax(i1) Then Let Utemp U(i1,t1) 1 X(i1,t1) 2 Xmax(i1) and X(i1,t1) Xmax(i1) Check If Utemp # Umax(i1,t1) Then U(i1,t1) Utemp Else U(i1,t1) Umax(i1,t1) and S(i1,t1) S(i1,t1) 1 Utemp 2 Umax(i1,t1) (6c) If X(i1,t1) # Xmin(i1) Then Let Utemp U(i1,t1) 1 X(i1,t1) 2 Xmin(i1) and X i1; t1 Xmin i1 Check If Utemp $ Umin(i1,t1) Then U(i1,t1) Utemp, Else U(i1,t1) Umin(i1,t1) and S(i1,t1) S(i1,t1) 1 Utemp 2 Umin(i1,t1) Step (7): If all reservoir and time periods have been treated then go to step (8), otherwise go to step (6a) Step (8): Check for the continuity equation at all units and all time periods. If satis®ed, then stop; otherwise go to step (1). 6. Details of the SAA for the LTHSP 6.1. Generating an initial feasible solution The proposed SAA requires a starting feasible solution, which satis®es all system constraints. This solution is randomly generated. The following steps are used in ®nding this starting solution: Step (1): (1a) Let i 1 and t 1
(1b) Let U i; t Umax i; t; X i; t Xmax i and S i; t 0 Step (2): Check if S(i,t) $ 0, then go to Step (4), Else let U1 Uin £ Umax i; t and X1 Xin £ Xmax i: Where Uin and Xin are a prespeci®ed parameters for ®nding initial feasible starting values for the variables (assumed 0.25) Step (3): U i; t R U1; Umax i; t X i; t R X1; Xmax i Go to step (2) Step (4): Check: If i N and t K then go to step (5), Else increase i and t and go to Step (1b). Step (5): Check the continuity equation for all units at all time periods. If satis®ed, then stop; Else go to step (1). 6.2. SA test The implementation steps of the SA test, as applied to each iteration in the algorithm, are described as follows [18]: Step (1): At the same calculated temperature, ckp ; apply the following acceptance test for the new trial solution. Step (2): Acceptance test: If Ej 2 Ei ; or If exp Ei 2 Ej =Cp $ R 0; 1; then accept the trial solution, set Xi Xj and Ei Ej : Otherwise reject the trial solution. Where Xi ; Xj ; Ei ; Ej is the SA current solution, the trial solution and their corresponding cost, respectively. Step (3): Go to the next step in the algorithm. 6.3. Step size vector adjustment In this work the step vector is updated jointly with temperature, according to the acceptance rate of the attempted moves at the previous temperature stage [23]. All its components are updates simultaneously. It must be mentioned here that a good step vector adjustment scheme must not yield steps that are too large with too many uncorrelated and unaccepted moves, whereas steps that are too small result in an incomplete exploration of the variation domain, with small and frequent increasing of the objective function. The following explains how the step vector is adjusted mathematically: Step (1): Calculate P i NTRACP i=NTR i; i 1; ¼; N: Where NTRACP(i) is the number of trial accepted when the variable i is changed. NTR(i) number of trials attempted by changing variable i. Step (2): If P(i) . PMAX, then STEP i STEP i £ STEPMAX
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46 Table 1 Comparison with DEC-DP and M-NORM
DEC-DP [2,9] M-NORM [9] Our SAA % Increase of SAA
Table 3 The storage (mm 3)
Total pro®t ($)
VWE ($)
28,145,330 28,227,174 28,799,568 2.32
10,403,170 10,414,598 10,307,287 22.02
If P(i) , PMIN, then STEP i STEP i £ STEPMIN Where PMAX, PMIN, STEPMAX and STEPMIN are parameters taken in our implementation as 0.2, 0.8, 0.8 and 1.2, respectively [23]. 6.4. Cooling schedule A ®nite-time implementation of the SA algorithm can be realized by generating homogenous Markov chains of ®nite length for a ®nite sequence of descending values of the control parameter. To achieve this, one must specify a set of parameters that governs the convergence of the algorithm. These parameters form a cooling schedule. The parameters of the cooling schedules are: an initial value of the control parameter decrement function for decreasing the control parameter and a ®nal value of the control parameter speci®ed by the stopping criterion, and a ®nite length of each homogenous Markov chain. In this work a polynomial-time cooling schedules is used in which the temperature is decreased based on the statistics of the trial solutions acceptance or rejection during the search. Details of the implemented cooling schedule are described in details in Ref. [18]. 6.5. Equilibrium test The sequence of trial solutions generated in the SAA at a ®xed temperature is stopped as soon as thermodynamic equilibrium, detected by some adequate condition, is Table 2 The discharge (mm 3) Month
1 2 3 4 5 6 7 8 9 10 11 12
45
Reservoir number 1
2
3
4
43.87 1037 1071.08 1068.41 968 1064.36 1037 0 0 0 0.13 0
411.97 1418 1290.33 1464.95 1322.77 1160.65 1105 127 510.31 314.82 491.05 363.09
572.7 1500.39 1330.8 1546.35 1338.41 1178.01 1109.42 405.17 695.25 470.73 636.48 495.24
2298.13 2704.04 2143.63 2450.56 2854.51 3160.28 1363.13 711.3 1695.27 3030.41 2220.36 1948.48
Month
1 2 3 4 5 6 7 8 9 10 11 12
Reservoir number 1
2
3
4
7473.13 7265.13 6772.05 6097.64 5394.64 4563.28 3719.28 3938.28 5039.28 6926.28 8076.15 8900.15
569.9 519.9 524.65 274.11 14.35 0.06 0.06 0.05 103.75 569.93 570 569.91
49.27 48.88 49.4 0 2.36 0 2.37 4.06 0.12 49.21 49.79 49.63
3419.58 3416.92 3414.09 2996.52 1782.42 58.82 0 1179.7 3419.68 3420 3419.12 3419.89
reached. Then the temperature and step vectors are suitably adjusted [23]. The test of equilibrium is done as follows. If the NTRACP(T ) , N1 £ n and NTR(T ) , N2 £ n then continue at the same temperature, otherwise end of temperature stage. Where NTRACP(T ), NTR(T ) are the number of trial accepted and attempted at temperature T, respectively, n the number of variables in the problem, N1 and N2 are the end temperature stage parameters (taken as 12 and 100) [23].
7. Numerical examples A computer model has been implemented based on the SAA described earlier and using the new proposed rules for generating random trial solutions. To test the model, an example system with four reservoirs connected in series on a river is solved. The scheduling time horizon is 1 year divided into 12 periods. The full data of this example are given in Refs. [2,9]. This Example was solved by combined decomposition-dynamic programming (DEC-DP) method [2,9] and by the minimum norm (MNORM) approach [9]. Table 1 shows a comparison of the results obtained by the DEC-DP approach, the M-NORM method and our proposed SAA. As shown in the table, our SAA achieves a higher total pro®t with a percentage increase of 2.32% relative to the DEC-DP approach and a value of 2.02% relative to the MNORM method. The value of water left in storage in the reservoirs at the end of the year (VWE) is also given in $ for the 3 methods. Detailed results for the Example test system are given in Tables 2 and 3. Table 2 shows the monthly discharge (releases) from the four reservoirs. Table 3 gives monthly storage in the four reservoirs. The solution obtained requires zero spillage from all reservoirs during the 12 months
46
A.H. Mantawy et al. / Electrical Power and Energy Systems 25 (2003) 41±46
except for the ®rst month, which requires 0.5, 0.6, 0.7 and 0.3 for the four reservoirs, respectively. 8. Conclusions In this paper we consider the LTHSP. The problem is a hard nonlinear optimization problem in continuous variables. The paper presents a new SAA for solving the LTHSP. The proposed algorithm introduces three new innovations; new rules for randomly generating feasible trial solutions, an adaptive cooling schedule and a method for variable discretization. This results in signi®cant reduction in the number of the objective function evaluations, and consequently less iterations are required to reach the optimal solution is obtained. The proposed SAA is applied successfully to solve a problem involving four series cascaded reservoirs hydroelectric power system. The results are compared with those obtained by DEC-DP [2,9] and MNORM [9]. Numerical results show a signi®cant increase of the total pro®t per year for the solved test system compared to previously obtained results using classical techniques. References [1] El-Hawary ME, Christensen GS. Optimal operation of electric power systems. New York: Academic Press, 1979. [2] Turgeon A. A decomposition method for the long-term scheduling of reservoirs in series. Water Resour Res 1981;17(6):1565±70. [3] Arvantidies NV, Rosing J. Composite representation of a multi reservoir hydroelectric power system. IEEE Trans PAS 1970; PAS-89(2):319±26. [4] Arvantidies NV, Rosing J. Optimal operation of multi reservoir systems using a composite representation. IEEE Trans PAS 1970; PAS-89(2):3327±35. [5] Turgeon A. Optimal operation of multi reservoir power system with stochastic in¯ows. Water Resour Res 1980;16(2):275±83. [6] Grygier JC, Stendinger J. Algorithms for optimization of hydropower system operation. Water Resour Res 1985;21(1):1±10. [7] Olcer S, Harsa C, Roch A. Application of linear and dynamic programming to the optimization of the production of hydroelectric power. Optimal Control Appl Meth 1985;6:43±56.
[8] Marino MA, Loaiciga HA. Dynamic model for multi reservoir operation. Water Resour Res 1985;21(5):619±30. [9] Christensen GS, Soliman SA. Optimal long-term operation of electric power systems. New York: Plenum Press, 1988. [10] Maceira MEP, Pereira MVF. Analytical modeling of chronological reservoir operation in probabilistic production costing. IEEE Trans Power Syst 1996;11(1):171±80. [11] Christoforidis M, et al. Long-term/mid-term resource optimization of a hydro-dominant power system using the interior point method. IEEE Trans Power Syst 1996;11(1):287±94. [12] Escudero LF, et al. Hydropower generation management under uncertainty via scenario analysis and parallel computation. IEEE Trans Power Syst 1996;11(2):683±9. [13] da Cruz Jr G, Soares S. Non-uniform composite representation of hydroelectric systems for long-term hydrothermal scheduling. IEEE Trans Power Syst 1996;11(2):702±7. [14] Wong KP, Wong SYW. Hybrid-genetic/simulated annealing approach to short-term multi-fuel constrained generation scheduling. IEEE Trans Power Syst 1997;12(2):776±84. [15] Guan X, et al. An optimization-based algorithm for scheduling hydrothermal power systems with cascaded reservoirs and discrete hydro constraints. IEEE Trans Power Syst 1997;12(4):1775±80. [16] Ponrajah RA, Galiana FD. Systems to optimize conversion ef®ciencies at Ontario hydro's hydroelectric plants. IEEE Trans Power Syst 1998;13(3):11044±50. [17] Yu Z, Sparrow FT, Bowen BH. A new long-term hydro production scheduling method for maximizing the pro®t of hydroelectric systems. IEEE Trans Power Syst 1998;13(1):66±71. [18] Mantawy AH, Abdel-Magid YL, Selim SZ. A simulated annealing algorithm for unit commitment. IEEE Trans Power Syst 1997; 13(1):197±204. [19] Aarts E, Korst J. Simulated annealing and Boltzman machines: a stochastic approach to combinatorial optimization and neural computing. New York: Wiley, 1989. [20] Cerny V. Thermodynamical approach to the traveling salesman problem: an ef®cient simulation algorithm. J Optim Theor Appl 1985;45(1). [21] Shokri Z, Selim K, Alsultan A. Simulated annealing algorithm for the clustering problem. Pattern Recog 1991;24(10):1003±8. [22] Metropolis N, Rosenbluth A, Rosenbluth M, Teller A, Teller E. Equations of state calculations by fast computing machines. J Chem Phys 1953;21:1087±982. [23] Siarry P, Berthiau G, Haussy J. Enhanced simulated annealing for globally minimizing functions of many continuous variables. ACM Trans Math Software 1997;23(2):209±28.