A Set of New Benchmark Optimization Problems for Water Resources Management

Water Resour Manage (2013) 27:3333–3348 DOI 10.1007/s11269-013-0350-z A Set of New Benchmark Optimization Problems for Water Resources Management Dimitrios K. Karpouzos & Konstantinos L. Katsifarakis Received: 18 August 2011 / Accepted: 31 March 2013 / Published online: 16 April 2013 # Springer Science+Business Media Dordrecht 2013 Abstract In this paper, we introduce four new benchmark problems, which are based on rather common optimization issues of water resources management. These problems have the following features: a) adjustable difficulty, to cover a wide range of common engineering problems b) physical background familiar to scientists working on water resources management c) known global optimal solution and known range of values of the objective function d) easy application and e) low computational volume (analytical solution of the respective groundwater flow model). First we calculate the optimal solutions of these problems and then we evaluate their difficulty and their suitability as benchmarking tools, based on theoretical considerations and on the performance of a genetic algorithm and a simulated annealing code in finding their optimal solutions. Results show that the proposed set of benchmark problems is useful for evaluating heuristic optimization codes in the field of water resources management. Keywords Optimization . Benchmark problem . Genetic algorithms . Simulated annealing . Water resources management 1 Introduction Optimization problems are encountered more or less frequently, in almost every scientific field (e.g. Floudas & Pardalos 2008). There are already many optimization methods, starting from simple differentiation (if the features of the problem permit its application) to sophisticated D. K. Karpouzos Department of Hydraulics, Soil Science and Agricultural Engineering, Faculty of Agriculture, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece e-mail: dimkarp@agro.auth.gr K. L. Katsifarakis (*) Division of Hydraulics and Environmental Engineering, Department of Civil Engineering, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece e-mail: klkats@civil.auth.gr 3334 D.K. Karpouzos, K.L. Katsifarakis heuristic techniques. No method is superior to all others for every problem (Wolpert & Macready 1997); or, as Reeves & Raw (2003) put it, “there is no royal road to optimization”. Moreover, there is no general rule for selecting the most suitable method, not even the best features of a certain method (e.g. Wu & Walski 2005). Experience, intuition and use of more than one optimization approaches may help with this process. Undoubtedly, though, there is a need for benchmark problems to compare optimization techniques (e.g. Matott et al. 2006, Younis & Dong 2010). There are already large sets of benchmark problems in the literature, with different features and degree of difficulty (e.g. Koon & Sebald 1995; Floudas et al. 1999; Gaviano et al. 2003; Shcherbina et al. 2003; Ali et al. 2005). There are also several approaches to the presentation and evaluation of comparative results (e.g. Youssef et al. 2001; Dolan & Moré 2002). In any case, to select the proper method for a real world application, one should use a benchmark problem with similar degree of difficulty (Tang et al. 2007). This paper introduces four new benchmark problems of adjustable difficulty, based on rather common aspects of water resources management problems, namely groundwater pumping cost (e.g. Fragoso et al. 2009; Chiu et al. 2010) and network amortization cost. First, we calculate their optimal solutions and then we evaluate their range of difficulty and their suitability as benchmarks for heuristic optimization techniques, based on theoretical considerations, but also on the performance of genetic algorithms and simulated annealing in their solution. These two methods are extensively used in most scientific fields, including water resources management (e.g. Teegavarapur & Simonovic 2002; Tsai et al. 2009). We opted for simple versions of the methods, not adapted to the known optimal solutions of the investigated problems, to arrive at impartial estimates. 2 Desirable Features of Benchmark Problems Our aim is to produce benchmark problems having the following features: (a) Adjustable difficulty, to cover a wide range of common engineering problems. (b) Physical background familiar to scientists working in the field of water resources management. (c) Known global optimal solution and known range of values of the objective function. (d) Easy application. (e) Low computational volume. In many cases, the objective function consists mainly of the solution of the respective physical problem. In heuristic methods, such as genetic algorithms and simulated annealing, the objective function must be evaluated hundreds or even thousands of times. It is very important then that the respective computational volume is as low as possible. 3 The Proposed Benchmark Problems The proposed benchmark problems have most of the aforementioned assets. Although they have common features with the benchmark problems described by Mayer et al. (2002), their computational load is much lower, since the respective groundwater flow problems are solved analytically. Moreover, the optimal solutions are known exactly for the three of them. Problem 1. Minimize the cost of pumping a given total well flow rate QT from an “infinite” aquifer with two zones of different transmissivities T1 and T2 (Fig. 1). Four A Set of New Benchmark Optimization Problems 3335 C(600,1200) D(-600,1200) T2 T1 W3 W4 c W1 b W2 Zone 1 Zone 2 B(600,0) A(-600,0) Fig. 1 The two-zone aquifer of the benchmark problems existing and two new wells will be used for this task. As shown in Fig. 1, the existing wells (W1 to W4) are located at the vertices of a square, whose side is equal to c, while the new ones should be constructed inside the square area ABCD, (AB=BC=b), shown with broken line. The solid line in Fig. 1 represents the interface between the two aquifer zones and coincides with the y-axis, while AB lies on x-axis. To standardize the problem, we set b=1,200 m and c=300 m. Then the coordinates of the existing wells W1 to W4 are (−150, 450), (150, 450), (150, 750) and (−150, 750), respectively. For the new wells, x can vary between −600 and 600, while y ranges from 0 to 1,200. Regarding the features of the aquifer, T1 =0.002 m2/s and T2 =0.001 m2/s, while the required total flow rate QT =200 l/s. The radius r0 of each well and the radius of influence of the system of wells R, which are used in subsequent equations, are equal to 0.25 m and 2,000 m, respectively (the latter is definitely larger than the maximum possible distance between wells). Pumping cost is defined as: C1 ¼ A N X QI h I ð1Þ I¼1 where N is the total number of wells (in our case N=6), QI is the flow rate of well I, hI the distance between water level at well I and a predefined level (e.g. highest ground elevation) and Α is a constant, depending on energy cost and duration of pumping. The “cost” function C1 that should be minimized can be reduced to: C1 ¼ A N X QI sI ð2Þ I¼1 where sI is the hydraulic head level drawdown at well I, which is due to pumping. In Eq. (2), QI are expressed in m3/s, sI in m, while A is set to 1,000 for convenience. At any point of the flow field, s depends on: a) its distance from the wells, b) the well flow rates and c) the aquifer features and boundaries. Given the set of QI and the coordinates of the wells, the 3336 D.K. Karpouzos, K.L. Katsifarakis respective sI can be calculated analytically, using: a) the method of images (e.g. Bear 1979), which introduces imaginary wells, symmetrical to the real ones with respect to the zone interface, in order to fulfil the boundary condition there and b) the superposition principle. Let us suppose that wells 1 to K are in zone 1, while the rest of them (K+1 to N) lie in zone 2. Then, sI for the wells 1 to K is given as: K P 1 2pT1 sI ¼ 1 pðT1 þT2 Þ J ¼1 N P QJ ln rRIJ QJ ln J ¼Kþ1 T1 T2 2pT1 ðT1 þT2 Þ K P QJ ln J ¼1 rIj R ð3Þ rIJ R while for wells, K+1 to N, which lie in zone 2, sI is given as: 1 2pT2 sI ¼ J ¼Kþ1 K P 1 pðT1 þT2 Þ N P QJ ln rRIJ T2 T1 2pT2 ðT1 þT2 Þ N P QJ ln J ¼Kþ1 rIj R ð4Þ rIJ R QJ ln J ¼1 In Eqs. (3) and (4) the real wells are denoted with capital letters, while their images with lower case ones. Moreover, R denotes the radius of influence of the system of wells, rIJ the distance between wells I and J (therefore rIJ =rJI) and rIj the distance between well I and imaginary well j (therefore rIj =rjI). The value of rII in particular, is taken equal to the radius of well I and is denoted by r0 for all I. It should be noted that summation is done for Q only. Using Eqs. (3) and (4) to express sI in Eq. (2), the following expression for C1 is obtained (Katsifarakis & Tselepidou 2009): C1 ¼ þ A 1 pðT1 þT2 Þ A N P I¼Kþ1 þ 1 pðT1 þT2 Þ K P QI I¼1 N P  J ¼Kþ1  QI K P J ¼1 1 2pT1 QJ ln 1 2pT2 K P 2 QJ ln rRIJ þ 2pTT11ðT1TþT 2Þ J ¼1  QJ ln J ¼1 rIj R rIJ R N P J ¼Kþ1  QJ ln K P 1 QJ ln rRIJ þ 2pTT22ðT1TþT 2Þ N P J ¼Kþ1 ð5Þ QJ ln rIj R rIJ R So, the objective function C1 of the optimization problem includes ten variables of two different kinds, namely six well flow rates and four coordinates (x, y values for the new wells, used to calculate the respective distances r). Moreover, the following constraints should be respected: 600  xI  600 0  yI  1200 N X I¼1 QI ¼ QT ð6Þ ð7Þ ð8Þ A Set of New Benchmark Optimization Problems 3337 Observance of constraints (6) and (7) should be checked for I=5 and 6, namely for the coordinates of the new wells only. If the coordinates of all wells are known, Katsifarakis & Tselepidou (2009) proved that the well flow rate values QI, which minimize C1, can be found analytically by solving a linear system of N equations and N unknowns. Equations 1 to K have the following form: K P QJ J ¼1 N P  QJ J ¼Kþ1 1 2pT1  ln rMJ R T1 T2 2pT1 ðT1 þT2 Þ ln rMj R rNJ 1 ln rMJ R þ 2pT2 ln R 1 p ðT1 þT2 Þ  þ pðT11þT2 Þ ln rRNJ þ  rNj 1 þ 2pTT22ðT1TþT ln ¼0 R 2Þ ð9Þ while equations K+1 to N-1 read: K P QJ J ¼1 N P  QJ  rNJ 1 ln rMJ R þ p ðT1 þT2 Þ ln R þ 1 p ðT1 þT2 Þ  J ¼Kþ1 T2 T1 2pT2 ðT1 þT2 Þ 1 2pT2 ln ð10Þ rNJ 1 ln rMJ R þ 2pT2 ln R rMj R 1 þ 2pTT22ðT1TþT ln 2Þ rNj R  ¼0 Finally, the N-th equation, which completes the system, is the constraint that well flow rate values should fulfil, namely Eq. (8). Calculation of C1 is then straightforward. In our test problem, the coordinates of the two new wells (W5 and W6) are not known a priori. It is known, though, from groundwater hydraulics, that: a) Given the well flow rate QI, the hydraulic head level drawdown sI at well I is smaller when it is placed in the zone with the larger transmissivity and it generally decreases with the distance from the interface. This can be inferred from Eqs. (3) and (4) for sI, where the respective TJ is in the denominator and b) Given the distribution of QI, the respective sI decrease with the distances among the wells. Then, most probably, the best locations for the two new wells are A and D, since T1 >T2. Setting x5 =−600, y5 =0, x6 =−600 and y6 =1,200, the aforementioned system of six equations and six unknowns results in the set of optimal well flow rate values which appears in Table 1 together with well coordinates. The respective “cost” is C1 =9486.81. If the aquifer were uniform, namely if T1 =T2, the optimal locations of the new wells would be A and C (or B and D), which maximize the distances among wells. For the given transmissivity values, this combination of well locations results in C1 =10329.1>9486.8. This confirms that the optimal solution coincides with our initial choice, given in Table 1. So, there is one global optimum, with C1 =9486.1 and two local optima, with C1 =10329.1. A third, weaker, local optimum arises when the new wells are placed at B and C, namely at the corners of zone 2. The respective “cost” value is C1 =11628.1. Moreover, a pair of local optima has been identified with both wells in zone 1, namely either y5 =y6 =0 or y5 =y6 = Table 1 Optimal set of well flow rates and coordinates of Problem 1 (6 wells) Well 1 xI (m) yI (m) QI (l/s) 2 3 4 −150 150 150 −150 −600 −600 450 450 750 750 0 1,200 18 18 32.5 5 32.5 6 49.5 49.5 3338 D.K. Karpouzos, K.L. Katsifarakis 1,200, while one of the x values is −600 and the other approximately −216.5. The respective “cost” value is C1 =10302.24. It should be mentioned that the difference between the C1 values at the global and the local optima depends on the ratio T1/T2. We have chosen a ratio that leads to a quite small difference. The complexity of the objective function can be increased, by increasing the number of existing wells. It is safe to add them along the line y=600 (with -300≤x≤600), in order not to affect the optimal location of the new wells. If the additional well is located at (−150,600) or at (−300,600), the optimal flow rate distributions, shown in Table 2, result in minimum cost values of 9158.11 and 9021.07 respectively. Problem 2. Find the optimal solution for pumping a given total well flow rate QT from the aquifer of Fig. 1, as in Problem 1. The number of new wells can be either one or two, though, and the cost function reads: C2 ¼ 6 X QI sI þ AMORT ð11Þ I¼1 where AMORT represents the annual amortization cost of the construction of the second new well. Hence AMORT=0, if Q5 =0 or Q6 =0, namely when only one new well is constructed and used, while it takes a given positive value, if both Q5 and Q6 are larger than 0. We expect that problem 2 is harder than problem 1, since the cost function exhibits a step. If only one new well is constructed and used, its best location is A (or D). Solving the pertinent system of Eqs. (8), (9) and (10), we end up with C2 =11594.25, larger than 9486.8, namely the respective cost for six wells, as expected. Thus the optimal solution depends on the value of AMORT. If AMORT<11594.25−9486.8=2107.45, two new wells should be constructed and used and the global minimum of C2 is 9486.8+AMORT, while 11594.25 is a local minimum. If, on the contrary, AMORT>2107.45, then only one new well should be constructed and the global minimum of C2 is 11594.25. It is obvious that the optimization problem becomes harder, as the value of AMORT approaches 2107.45, since the local minimum approaches the value of the global one. Thus, the difficulty of the benchmark problem can be adjusted through the selection of AMORT value. Problem 3. Minimize the annual amortization cost of the pipe network carrying water from a set of wells to a central water tank TC, located at (600, 600), namely on the boundary of the available area and shown as a small square in Fig. 2a. Six existing and two new wells will be used. As shown in Fig. 2, the locations of the existing wells (W1 to W6) are those of the optimal solution of Problem 1. Table 2 Optimal set of well flow rates and coordinates of Problem 1 (7 wells) Well 1 xI (m) yI (m) QI (l/s) 2 3 4 5 6 −150 150 150 −150 −150 −600 −600 450 450 750 750 600 0 1,200 26 16 16 26 23 xI (m) −150 150 150 −150 −300 yI (m) 450 450 750 750 600 QI (l/s) 25.84 15.9 15.9 25.84 26.42 7 46.5 46.5 −600 −600 0 1,200 45.05 45.05 A Set of New Benchmark Optimization Problems 3339 (600,1200) (-600,1200) W6 T1 T2 W3 W4 TC W7 W8 W1 W2 Zone 2 Zone 1 W5 (-600,0) (600,0) (a) (600,1200) (-600,1200) (600,1200) (-600,1200) W6 W6 T1 T1 T2 W4 T2 W4 W3 W3 TC TC W7 W1 W1 W2 Zone 2 Zone 1 W2 Zone 2 Zone 1 W5 W5 (-600,0) (600,0) (-600,0) (b) (600,0) (c) Fig. 2 Well layouts for the shortest (a), second shortest (b) and third shortest (c) pipe network (Problem 3) Since the amortization cost is proportional to the length of the pipe network, the “cost” function is reduced to: C3 ¼ 8 X LI ð12Þ I¼1 where LI is the length of the pipe that carries water away from well I. The decision variables are the coordinates of the 2 new wells only, while the global optimum can be calculated analytically, based on the property of the Fermat point of a triangle to minimize the sum of the distances from its three vertices. So, in our case, one new well should be placed at the Fermat point of triangle W2TCW3, while the other at the Fermat point of triangle W1W2W4 (or W1W3W4), as shown in Fig. 2a. The coordinates of the aforementioned Fermat points appear in last section of Table 3, while the respective smallest network length is equal to 3340 D.K. Karpouzos, K.L. Katsifarakis 2562.15 m. One should mention that the Fermat points of the triangles W1W2W5 and W3W4W6 coincide with vertices W1 and W4 respectively, since the respective angles exceed 120o. The interest in using Problem 3 as benchmark stems from the number and value of local optima. If one new well is placed at the Fermat point of the triangle W2TCW3, while the other at any point of the network shown in Fig. 2b, the “cost” value is equal to 2582.6, just 0.8 % larger than the global optimum. Moreover, the number of these strong local optima is theoretically infinite (practically very large, depending on the network length and the step in coordinate values). A second set of “infinite” local optima results when both new wells lie on any segment of the shortest network which connects the 6 existing wells with the water tank (shown in Fig. 2c, mention that pipe W2TC could be substituted by W3TC). The respective “cost” value is 2647.14, namely 3.3 % larger than the global minimum. Problem 4. Minimize the sum of two cost items: a) the annual cost for pumping a given total well flow rate QT from the aquifer considered in Problem 1 and b) the annual amortization of the construction cost of the pipe network, carrying water from the wells to a central water tank TC. Six existing and two new wells will be used for this task, as in Problem 3 and as shown in Fig. 2. The “cost” function in this case has the following form: C 4 ¼ A1 8 X I¼1 QI sI þ A2 8 X LI ð13Þ I¼1 where A1 and A2 are cost coefficients. We have set A1 =1,000 and A2 =6, after a number of trials, in order to balance the two terms of the “cost” function. The exact optimal solution of this problem is not known a priori. As a starting point we have used the best solution of Problem 3. The respective smallest network length is equal to 2562.15 m, which, multiplied by A2, gives the respective smallest “network cost” Cn =15372.9. The optimal flow rate distribution for this well layout can be calculated by solving a system of eight equations and eight unknowns, as in Problem 1. The results appear in Table 3, together with well coordinates. The respective “pumping cost” is Cp =8,961. Then C4 =Cn +Cp =24333.9. So we know that the global optimum does not exceed this value. Further investigation, using the methods of genetic algorithms and simulated annealing, led to the approximate identification of fourteen local optima (twelve of them being symmetric pairs), with C4 values smaller than 24333.9. Possible local optima have been initially identified through a large number of algorithm runs, and characterized as such if the respective solutions correspond to: a) well layouts with distinct geometrical features, deriving from different compromises between the optimal pumping and network layouts and b) combination of flow rates resulting in almost equal hydraulic head level drawdowns at the wells. Finally, we have conducted local search around each of these points, to verify whether they are local minima. We believe that the value C4 =24030.73, which corresponds to the best pair, can be considered as the global optimum. It represents the best compromise between pumping and network cost. One additional well is placed close to the Fermat point of the triangle W2TCW3, substantially reducing network cost, and the other on the pipe connecting existing wells W5 and W1 (or W6 and W4), namely in the higher transmissivity zone, comparatively away from existing wells and with zero additional network cost. The C4 value, which corresponds to the next local optimum, is 24055.2, namely just 0.1 % larger than the global optimum, while the well layout is quite different, as one of the new wells is placed on the pipe connecting wells W5 and W1 and the other on the pipe connecting W6 and W4. Well coordinates and flowrates for the global and local optima appear in Table 3, while the network layouts are shown in Fig. 3. A Set of New Benchmark Optimization Problems 3341 Table 3 Well coordinates and flow rates for the global and local optima of Problem 4 Well 1 2 3 4 5 6 7 8 fitness Global optimuma xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 −482.5 253.4 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 115.6 600.8 24.3 QI (l/s) Local Optimum 1 12.8 13.3 26.1 36.0 43.1 30.9 13.6 xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 −484.0 −484.1 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 1088.5 111.5 QI (l/s) 23.0 13.9 13.9 23.0 34.0 34.0 29.1 29.1 24030.73 24055.20 Local Optimum 2a xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 600.0 −486.2 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 600.0 1086.6 QI (l/s) 25.3 Local Optimum 3 13.2 12.7 23.6 41.7 34.7 19.0 29.8 xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 600.0 245.7 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 600.0 600.2 QI (l/s) 27.4 12.6 12.6 27.3 44.5 44.6 19.1 11.9 24121.19 24180.88 Local Optimum 4a xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 −383.2 −529.0 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 985.9 1130.4 25.5 QI (l/s) Local Optimum 5a 14.8 14.0 22.3 42.0 30.7 24.8 26.1 xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 −83.0 −486.6 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 515.3 1088.9 QI (l/s) 21.8 13.9 13.9 22.9 42.5 35.6 18.8 30.6 24255.18 24296.89 Local Optimum 6a xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 −488.8 −82.4 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 1091.0 682.0 24.7 QI (l/s) Local Optimum 7a 14.5 13.6 20.5 42.8 35.4 30.3 18.3 xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 255.9 −85.3 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 599.4 684.3 QI (l/s) 26.2 13.4 12.9 23.3 45.4 45.1 13.9 19.9 24309.49 24325.25 Shortest pipe networka a xI (m) −150.0 150.0 150.0 −150.0 −600.0 −600.0 236.6 −86.6 yI (m) 450.0 450.0 750.0 750.0 0.0 1200.0 600.0 513.4 QI (l/s) 23.2 13.0 13.5 26.3 45.2 45.5 13.5 19.9 24333.90 The symmetrical solution with respect to the axe y=600 is also an optimum with the same fitness value 4 Evaluation of the Relevant Difficulty of the Benchmark Problems Assessing (and even defining) the difficulty of optimization problems is not an easy task and has attracted the interest of many scientists (e.g. Goldberg 2002), who have developed sets of evaluation and comparison criteria (e.g. Kollat & Reed 2006). An important source of 3342 D.K. Karpouzos, K.L. Katsifarakis W1 Fig. 3 Well layouts for the global and local optima (Problem 4) difficulty is multimodality (Singh & Deb 2006; Deb & Saha 2010). All of the proposed benchmark problems share the feature of multimodality. Problem 3 exhibits a very large set of local optima, much larger than the other three. Problem 2 exhibits a step in its evaluation function. Problem 4 has a more complicated evaluation function, while it exhibits more local optima than Problems 1 and 2. So, we believed that Problem 1 is the easiest one and we anticipated that Problem 3 would be the hardest. To further assess the difficulty of the four benchmark problems, we have used simple versions of two of the most common heuristic optimization methods, namely genetic A Set of New Benchmark Optimization Problems 3343 algorithms and simulated annealing. The first method deals with a population of solutions, while the second one processes a single solution only. In our simple binary genetic algorithm code the evaluation process includes the flow simulation model for each benchmark problem, while selection is performed by the tournament method. Moreover, the best chromosome of each generation is separately preserved through the selection process (elitist approach). The code uses also the following operators: typical one-point crossover (recombination of the features of pairs of chromosomes), mutation (change of the value of a gene from 1 to 0 or vice versa) and antimetathesis, which has been introduced by Katsifarakis et al. (1999) and is used interchangeably with mutation. Simulated annealing (SA), first proposed by Kirkpatrick et al. (1983) and Cerny (1985), is an efficient method for finding the global optimal (or a suboptimal) solution to multidimensional optimization problems. Simulated annealing algorithm begins by generating a single initial solution (instead of a population of solutions used in GA) at random, which is evaluated by the problem-specific objective function. In our case, a solution is a real vector consisting of ten parameters (six well flow rates and four coordinate values), while the objective function requires the simulation of the groundwater flow. Although there is a number of variants of simulated annealing algorithm found in relative literature (Ingber 1993; Cunha 2003; Ji et al. 2006; Suman & Kumar 2006; Liu et al. 2008), the main structure is almost preserved and comprises the three following operators: a temperature cooling schedule T, a function g for generating a perturbation and a state transition with an acceptance probability p (Nam et al. 2004). The SA, implemented in this paper, uses the following annealing scheme (Yang et al. 2005): T¼ g¼  kmax k q ð1 þ μÞjbj μ p ¼ exp  ; q>0 1 signðbÞ Δf T jf jef 1  ð14Þ ð15Þ ð16Þ where k and kmax are the iterations index and the maximum number of iterations respectively, q a quenching factor, b a random vector having uniform distribution U[-1,1] and μ the parameter of the inverse μ-law function g that is increased according to a rule μ ¼ 1 10ð100T Þ : Finally, Δf, f, and ef are the change, the value, and the relative tolerance of the objective function, respectively. 5 Results By virtue of the low computational volume per evaluation of the objective function, the genetic algorithm and simulated annealing codes were run hundreds of times, to get meaningful statistical results. 3344 D.K. Karpouzos, K.L. Katsifarakis 5.1 Performance of Genetic Algorithms We have used the simple genetic algorithm, outlined in 4, since our aim is to evaluate the suitability of the proposed problems as benchmarks. Seven genes are needed for each QI, which is allowed to vary between 0 and 127, and 11 genes for each coordinate. Thus the chromosome length SL for the first problem equals 86. The coordinate values, which result from chromosome decoding, are scaled to observe constraints (5) and (6), namely to range between −600 and 600. Moreover, since the sum SQ of the QI is not equal to QT, their values are multiplied by QT/SQ. Finally, we have used mutation/antimetathesis probabilities slightly larger or equal to 1/SL and we have set crossover probability to 0.6 in all problems, based on a number of trials with the first problem. In the first problem, we have set the population size (PS) equal to 25 and the number of generations (NG) equal to 80. Thus the code evaluated the objective function, namely it solved the groundwater flow problem, 2,000 times. The code has been successful in more than 96 % of the runs, namely it hits exactly or approximates satisfactorily the global optimum (the difference being less than 1 %). Very few times it is trapped around one of the two local optima of zone 1. Addition of one well affects performance of the genetic algorithm, but rather slightly. If the well is placed at (−150,600) the code approximates the global optimum satisfactorily in around 91 % of the runs, but very rarely hits it exactly. If the well is placed at (−300,600), success rate falls to 89 %. In cases of failure the code is either trapped at a local optimum of zone 1 or it places one well at x=600. To solve the second problem, one more gene has been added to each chromosome, to define the number of new wells (gene values of 0 and 1 indicate one and two new wells respectively). As the problem is considered more difficult, we have set PS=40 and NG=100, namely the code solved the hydraulic problem 4,000 times in each run. We tried different values of AMORT, to check whether the code ends up with the correct solution, namely with two and one additional wells, if AMORT is smaller and larger than 2107.45, respectively. Results are shown in Table 4. It can be seen that when AMORT is well below 2107.45, the genetic algorithm code tends correctly to the 2-well solution. Its efficiency is reduced, though, as AMORT approaches the limiting value. But if AMORT exceeds 2107.45, the code always tends correctly to the 1-well solution. The chromosome length for the third problem is 44, namely substanstially smaller than that of the first problem. Nevertheless, the rate of success fell to 55 %. In the rest of the runs the genetic algorithms was trapped at a local optimum with C3 =2582.6 or, very rarely (around 1 % of the runs), at a local optimum with C3 =2647.14 Table 4 Genetic algorithm and Simulated annealing performance for different AMORT values of Problem 2 Genetic algorithm performance AMORT Correct number of wells Success (%) 1,000 2 1,500 2 2,000 2 2,215 1 2,715 1 100 65 8 100 100 1,000 1,500 2,000 2,215 2,715 2 2 2 1 1 94 72 51 75 95 Simulated annealing performance AMORT Correct number of wells Success (%) A Set of New Benchmark Optimization Problems 3345 The chromosome length SL for the fourth problem equals 100. We have tried different combinations of population size and number of generations, their product ranging from 2,000 to 5,000. The algorithm managed to hit exactly or to reach approximately the global optimum in around 28 % of the runs, while in around 65 % of them it was trapped at the strongest local optimum. Interestingly, the increase of the product PSxNG improved the accuracy of the successful runs, but its effect on the rate of success was not significant. 5.2 Performance of Simulated Annealing The simulated annealing algorithm follows the scheme given by Eqs. (14)–(16). Each solution in its general form is represented by a state vector of real parameters z=[Qi, xj,, yj], where index i covers the total set of wells and j the new wells only, in a domain defined by un upper bound u=[200, 600, 1200] and a lower bound l=[0,−600, 0]. The size of the random perturbation Δz, defined in the domain of z, is the element-wise product of the function g returned vector and the vector (u-l). In case that the new generated solution z΄= z+Δz, exceeds the upper or/and lower bound, it is confined inside the admissible region. Regarding the constraint on the sum SQ of the QI, it is handled in the same way as it is treated in the genetic algorithm. The quenching factor q in the following applications is taken equal to one. Finally, in all applications the maximum number of iterations kmax (objective function evaluations) is set equal to those of the respective GA code. In the first problem, the algorithm had an easy task, since it managed to find the optimum in all runs. Furthermore, the addition of one more well did not alter the algorithm’s convergence behaviour, since the optimal was also found in almost all trials. The code that has been employed in the first problem is also used to solve the second problem. The only change concerns the adjustment of the form of the objective function, in order to incorporate the annual amortization cost in case a solution indicates construction of two new wells (Q5 and Q6 values are different from zero) instead of one (Q5 or Q6 values are equal to zero). The investigation of different values of AMORT is performed in a similar way to the one adopted in Section 5.1. The results of this analysis are presented in Table 4. It can be observed that for all AMORT values the performance of simulated annealing is quite stable and satisfactory (more than 50 %). The more the AMORT values diverge from the critical value 2107.45 the more the performance of simulated annealing algorithm is improved. The worst score of the code corresponds to an AMORT value close to the critical one, which is the most challenging case. The third problem concerning the optimal network layout of an eight wells system, requires only four real bounded variables – the four coordinates of the two new wells. Nevertheless, the performance of SA in reaching the global minimum (2562.15) is drastically reduced. Successful runs are just a small subset (11 %) of a large number (>1,000) of algorithm trials. Local optima corresponding to values 2582.6 and 2647.14 dominate the 60 % and 29 % of SA runs respectively. The infinite number of variable combinations leading to these two local optima situated in the vicinity of the global minimum in conjunction with the single starting point property of SA can justify the above low performance score. In order to further investigate the difficulty of this benchmark problem, a local non-linear optimization algorithm is also tested. For this purpose, Sequential Quadratic Programming (SQP), as implemented in Matlab (R2010a) optimization function “fmincon”, is properly coded and applied to solve this problem. The SQP succeeded to find the global minimum only in 5 % of runs, connoting the domination of local minima in fitness landscape. 3346 D.K. Karpouzos, K.L. Katsifarakis The fourth benchmark problem is properly formulated in order to find a real state vector z of twelve variables (i.e. eight flow rates and four well coordinates) that minimizes the weighted sum of pumping and network costs. Due to the various local optima that characterize this specific problem the success in reaching the global optimum is rather limited, namely around 27 %. Local optima form a condensed fitness landscape that can significantly augment the possibility an optimization algorithm is trapped in. Besides that, most of the local optima correspond to two different combinations of well positions, symmetrical to the central horizontal axis y=600, thus increasing their number and subsequently the hardness of the optimization task. The local optimum of 24055.20 was found to be the strongest one, as SA was stacked in it for more than half of the runs. This could be explained by the fact that this local optimum is the nearest to the global one and seems to have a quite large attraction region with a rather mild slope, compared to the other local optima. 6 Conclusions This study proposed and evaluated a set of four benchmark problems through the application of two heuristic optimization algorithms - genetic algorithms and simulated annealing. Both methods performed quite satisfactorily in the first problem, even when the number of existing wells was increased. Simulated annealing performed slightly better, since it predicted the global optimum almost exactly in every run. In any case, we can classify problem 1 as rather easy for heuristic techniques. In the second problem, the existence of the step in the evaluation function plays the dominant role. Success of both methods depends heavily on the value of AMORT, which determines the difficulty of the benchmark problem. Performance of both techniques in the third problem verified our initial estimate that it is rather difficult. While the performance of GA remained above 50 %, that of SA sunk to 11 %. Moreover, the advanced local search technique, which has been additionally used, had a rate of success of 5 % only. Both methods had a rather low success in the fourth problem, identifying the global optimum in less than 30 % of the runs. The high multimodality, the intense local optimal density of fitness landscape and the increased number of dimensions affect the effectiveness of the optimization algorithms. While the difficulty of the aforementioned optimization problems varies, they share an important feature: low computational volume of the solution evaluation processes. This feature renders these problems particularly attractive for comparing and tuning heuristic optimization techniques, which involve repetitive calculations of the objective function. If, for instance, the objective function includes numerical simulation of a groundwater flow field, it is crucial to know the lowest values for the population size and the number of generations that will allow an accurate estimate of the optimal solution. This can be achieved by applying the investigated optimization technique to one or more of the proposed benchmark problems. Moreover, these proposed benchmark problems can serve in evaluating additional features of genetic algorithm or simulated annealing codes, which can be rated in comparison with the simple codes used in this paper. Problem 1 can serve for an initial evaluation of code efficiency. Problem 2 allows comparative rating in the presence of steps in the objective function, while Problem 3 can be used when the solution field exhibits a very large number of local optima. Finally, Problem 4 would be useful with multimodal fields combined with more complex evaluation functions. A Set of New Benchmark Optimization Problems 3347 Finally, scientists working in the field of water resources management are very familiar with the physical background of the proposed benchmark problems. This feature is expected to facilitate their use in the related research area. References Ali MM, Khompatraporn C, Zabinsky ZB (2005) A numerical evaluation of several stochastic algorithms on selected continuous global optimization test problems. J Global Optim 31(4):635–672 Bear J (1979) Hydraulics of Groundwater. McGraw-Hill. Cerny V (1985) Thermodynamical approach to the travelling salesman problem: An efficient simulation algorithm. J Optim Theory Appl 45:41–51 Chiu Y-C, Nisshikawa T, Yeh WW-G (2010) Optimal Pump and Recharge Management Model for Nitrate Removal in the Warren Groundwater Basin. California J Water Resour Plann Manage 136(3):299–308 Cunha MC (2003) On Solving Aquifer Management Problems with Simulated Annealing Algorithms. Water Resour Manage 13:153–170 Deb K, Saha A (2010) Finding multiple solutions for multimodal optimization problems using multiobjective evolutionary approach. GECCO 2010:447–454 Dolan ED, Moré JJ (2002) Benchmarking optimization software with performance profiles. Math Program 91:201–213 Floudas CA, Pardalos PM (eds.) 2008 Encyclopedia of Optimization. 2nd edn. Kluwer Academic Publishers. Floudas CA, Pardalos PM, Adjiman C, Esposito W, Gümüs Z, Harding S, Klepeis J, Meyer C, Schweiger C (1999) Handbook of Test Problems in Local and Global Optimization. Kluwer Academic Publishers, Dordrecht Fragoso T, Cunha MC, Lobo-Ferreira JP (2009) Optimal pumping from Palmela water supply wells (Portugal) using simulated annealing. Hydrogeol J 17(8):1935–1948 Gaviano M, Kvasov DE, Lera D, Sergeyev YD (2003) Algorithm 829: Software for generation of classes of test functionswith known of local and global minima for global optimization. ACM Trans Math Softw 29(4):469–480 Goldberg DE (2002) Design of innovation: lessons from and for competent genetic algorithms. Springer. Ingber L (1993) Simulated annealing: practice versus theory. Math Comput Modell 18:29–57 Ji M, Jin Z, Tung H (2006) An improved simulated annealing for solving the linear constrained optimization problems. Appl Math Comput 183:251–259 Katsifarakis KL, Tselepidou K (2009) Pumping cost minimization in aquifers with regional flow and two zones of different transmissivities. J Hydrol 377(1–2):106–111 Katsifarakis KL, Karpouzos DK, Theodossiou N (1999) Combined use of BEM and genetic algorithms in groundwater flow and mass transport problems. Engin Anal Bound Elem 23(7):555–565 Kirkpatrick S, Gelatt CD, Vecchi MP (1983) Optimization by Simulated Annealing. Sci 4598:671–680 Kollat JB, Reed PM (2006) Comparison of multi-objective evolutionary algorithms for long-term monitoring design. Adv Water Resour 29(6):792–807 Koon GH, Sebald AV (1995) Some interesting test functions for evaluating evolutionary programming strategies. Proc fourth Annual Conf on Evolut. Program, San Diego, California, 1-3 March 1995, Evolutionary Programming IV, MIT Press, 479-499. Liu JS, Caley AJ, Waddie AJ, Taghidech MR (2008) Comparison of simulated quenching algorithms for design of diffractive optical elements. Appl Opt 47(6):807–816 Matott LS, Bartelt-Hunt SL, Rabideau AJ, Fowler KR (2006) Application of heuristic optimization techniques and algorithm tuning to multilayered sorptive barrier design. Environ Sci Technol 40:6354–6360 Mayer AS, Kelley CT, Miller CT (2002) Optimal design for problems involving flow and transport in saturated porous media. Adv Water Resour 12:1233–1256 Nam D, Lee JS, Park CH (2004) N-dimensional Cauchy neighbour generation for the Fast Simulated Annealing. IEICE Trans Inf Syst E87-D(11):2499–2502 Reeves CR, Raw JE (2003) Genetic algorthms-Principles and perspectives. Kluwer Academic Publishers. Shcherbina O, Neumaier A, Sam-Haroud D, Vu XH, Nguyen TV (2003) Benchmarking global optimization and constraint satisfaction codes. In: Global Optimization and Constraint Satisfaction. Lecture Notes in Computer Science 2861:211–222 Singh G, Deb K (2006) Comparison of multi-modal optimization algorithms based on evolutionary algorithms. GECCO 2006:1305–1312 3348 D.K. Karpouzos, K.L. Katsifarakis Suman B, Kumar P (2006) A survey of simulated annealing as a tool for single and multiobjective optimization. J Oper Res Soc 57:1143–1160 Tang Y, Reed PM, Kollat JB (2007) Parallelization strategies for rapid and robust evolutionary multiobjective optimization in water resources applications. Adv Water Res 30:335–353 Teegavarapur RSV, Simonovic SP (2002) Optimal Operation of Reservoir Systems using Simulated Annealing. Water Resour Manage 16:401–428 Tsai FTC, Katiyar V, Toy D, Goff RA (2009) Conjunctive Management of Large-Scale Pressurized Water Distribution and Groundwater Systems in Semi-Arid Area with Parallel Genetic Algorithm. Water Resour Manage 23:1497–1517 Wolpert DH, Macready WG (1997) No free lunch theorems for optimization. IEEE Trans Evol Comput 1(1):67–82 Wu ZY, Walski T (2005) Self-adaptive penalty approach compared with other constraint-handling techniques for pipeline optimization. J Water Resour Plann Manage 131(3):181–192 Yang WY, Cao W, Chung TS, Morris J (2005) Applied numerical methods using Matlab. Hoboken, New Jersey, John Wiley & Sons, Inc Younis A, Dong Z (2010) Trends, features, and tests of common and recently introduced global optimization methods. Eng Optim 42(8):691–718 Youssef H, Sait SM, Adiche H (2001) Evolutionary algorithms, simulated annealing and tabu search: a comparative study. Eng Appl Artif Intell 14(2):167–181