Cross Entropy multiobjective optimization for water distribution ...

Click Here

WATER RESOURCES RESEARCH, VOL. 44, W09413, doi:10.1029/2007WR006248, 2008

for

Full Article

Cross Entropy multiobjective optimization for water distribution systems design Lina Perelman,1 Avi Ostfeld,1 and Elad Salomons2 Received 11 June 2007; revised 13 June 2008; accepted 23 June 2008; published 10 September 2008.

[1] A methodology extending the Cross Entropy combinatorial optimization method

originating from an adaptive algorithm for rare events simulation estimation, to multiobjective optimization of water distribution systems design is developed and demonstrated. The single objective optimal design problem of a water distribution system is commonly to find the water distribution system component characteristics that minimize the system capital and operational costs such that the system hydraulics is maintained and constraints on quantities and pressures at the consumer nodes are fulfilled. The multiobjective design goals considered herein are the minimization of the network capital and operational costs versus the minimization of the maximum pressure deficit of the network demand nodes. The proposed methodology is demonstrated using two sample applications from the research literature and is compared to the NSGA-II multiobjective scheme. The method was found to be robust in that it produced very similar Pareto fronts in almost all runs. The suggested methodology provided improved results in all trails compared to the NSGA-II algorithm. Citation: Perelman, L., A. Ostfeld, and E. Salomons (2008), Cross Entropy multiobjective optimization for water distribution systems design, Water Resour. Res., 44, W09413, doi:10.1029/2007WR006248.

1. Introduction [2] The optimal design problem of a water distribution system is commonly defined as a single objective optimization problem of finding the water distribution system component characteristics which minimize the system capital and operational costs, such that the system hydraulic laws are maintained, and constraints on quantities and pressures at the consumer nodes are fulfilled. [3] However, in reality, the design problem of a water distribution system involves competing objectives, such as: minimizing cost, maximizing reliability, minimizing risks, minimizing deviations from specific targets of quantity, pressure, and quality, etc. The design problem is thus inherently of a multiobjective nature. In a multiobjective optimization problem there is not a single optimal solution but a set of compromised solutions that form the so-called Pareto optimal solution set. Incorporating multiple objectives in the optimal design of water distribution systems provides an improvement compared to using a single design objective as a broader range of alternatives is explored, thus making the design outcome more realistic. [4] The problem of water distribution systems optimal design has attracted numerous studies over the last four decades. The studies concentrated mainly on minimizing the construction and operational costs of the system, but also considered redundancy [e.g., Jacobs and Goulter, 1989],

1 Faculty of Civil and Environmental Engineering, Technion-Israel Institute of Technology, Haifa, Israel. 2 OptiWater, Haifa, Israel.

Copyright 2008 by the American Geophysical Union. 0043-1397/08/2007WR006248$09.00

reliability [e.g., Su et al., 1987; Duan et al., 1990; Walski, 1993], and water quality [e.g., Ostfeld and Shamir, 1996]. The use of multiobjective optimization for water distribution systems has started recently. Halhal et al. [1999] were the first to introduce a multiobjective procedure to solve a water distribution systems management problem. Minimizing network cost versus maximizing the hydraulic benefit served as the two conflicting objectives, with the total hydraulic benefit evaluated as a weighted sum of pressures, maintenance cost, flexibility, and a measure of water quality benefits. Kapelan et al. [2003] used a multiobjective genetic algorithm to find sampling locations for optimal calibration. Keedwell and Khu [2003] applied a hybrid multiobjective evolutionary algorithm to the optimal design problem of a water distribution system. The hybrid approach employed a nondominated sorting genetic algorithm coupled with a neighborhood search technique. Prasad and Park [2004] applied a nondominated sorting genetic algorithm for minimizing the network cost versus maximizing a reliability index. Babayan et al. [2005] used a multiobjective genetic algorithm to solve the design problem of a water distribution system under uncertainty. [5] In this manuscript a new method for multiobjective optimization is developed and demonstrated. The proposed methodology is based on the Cross Entropy (CE) algorithm [Rubinstein, 1999; Rubinstein and Kroese, 2004] coupled with selected features from existing multiobjective evolutionary algorithm techniques [Fonseca and Fleming, 1995]. [6] The proposed method extends the study of Perelman and Ostfeld [2007] on CE to multiobjective optimization of water distribution systems design. It is demonstrated using the New York tunnels (NYT) water system [Schaake and Lai, 1969], and an extension of the two loop network (TLN)

W09413

1 of 15

W09413

PERELMAN ET AL.: CROSS ENTROPY MULTIOBJECTIVE OPTIMIZATION

of Alperovits and Shamir [1977] to multiple loadings, pumping, and storage [Salomons, 2001].

2. Problem Formulation [7] A water distribution system can be modeled as a graph G (V, E) with the set of nodes V representing the sources, consumers and intermediate nodes, and the set of edges E representing the control elements (e.g., pumps, valves) and connecting pipes. Let K (V, E) be the directed incidence matrix of G defined as: kij ¼ þ1 If edge i is directed away from node j kij ¼ 1 If edge i is directed toward node j kij ¼ 0 If edge i is not adjacent to node j

ð1Þ

For ease of representation the methodology presented herein is for gravitational networks considering two competing objectives, but the method is capable of handling general systems. The optimal design problem of a gravity driven water distribution system, using the Hazen – Williams headloss equation, can be defined as: Minimize : F1 ¼ cðd; LÞ Maximize : F2 ¼ f h; hmin

minimization [Rubinstein, 1997]. Rubinstein [1999] suggested the CE algorithm for combinatorial optimization through modifying the rare events simulation estimation algorithm by: (1) employing the Kullback – Leibler [Kullback and Leibler, 1951] distance measure for solving the rare event simulation estimation problem, and (2) by utilizing a supplementary random mechanism which translates the optimization problem considered into an associated stochastic problem (ASP). The rationale and theory underlying the CE optimization method are briefly described below, starting from the problem of rare event estimations, then through invoking the Kullback – Leibler distance measure, and finally through utilizing an associated stochastic problem which needs to be solved. A complete outline of the CE method theory can be found in Rubinstein and Kroese [2004]. [10] Consider estimating the probability ‘(l): ‘ðlÞ ¼ PðZ lÞ

ð7Þ

where: Z = a random variable, and l = a real number. A straightforward attempt to estimate ‘(l) is through Monte Carlo sampling:

ð2Þ

RT 1 X \ ðZ i l Þ ‘^ðlÞ ¼ RT i¼1 i

ð8Þ

ð3Þ

Subject to: Kq þ qnode ¼ 0

ð4Þ

KT h DhðqÞ ¼ 0

ð5Þ

where: F1, F2 = objective functions, c = cost function, f = pressure deficiency function, d = vector of pipe diameters, L = vector of pipe lengths, h = vector of nodal heads, hmin = vector of minimum required nodal heads, q = vector of pipe flows, qnode = vector of nodal demands, and Dh (q) = vector of headlosses along the pipes, whose ith component equals to: Dhi ðqi Þ ¼ RPi ðdi ; CHWi ; Li Þqi jqi j0:852

W09413

ð6Þ

where: RPi = the resistance of pipe i (a function of pipe diameter di, Hazen – Williams headloss coefficient CHWi, and length Li), qi = the flow along pipe i which can become negative if flow reverses, and 0.852 = an empirical constant. [8] The objective function F1 in equation (2) is the cost of the water distribution system design, which is a nonlinear function of the pipe diameters, the link lengths, and the installation cost of the pipes. The objective function F2 in equation (3) is a measure of the hydraulic performance of the system, which is a function of the actual heads at the consumer nodes versus the minimal head requirements. In this study, the hydraulic performance is measured as the maximum pressure deficit below a minimum required pressure head at the demand nodes.

3. Cross Entropy for Combinatorial Optimization [9] The CE method originates from an adaptive algorithm for rare events simulation estimation based on variance

^ = the estimator of ‘(l), Zi = the ith (Z1,. . ., ZRT) where: ‘(l) generated sample from the probability density function (PDF) of Z, and \i = the ith Indicator function outcome, receiving the value of one if at the ith sample Zi l, and zero otherwise. [11] For a very low value of l, ‘(l) is very small and the estimation of ‘(l) is of a rare event. Applying crude Monte Carlo sampling to estimate ‘(l) [i.e., equation (8)] is practically impossible (e.g., about 1010 samples are required to estimate ‘(l) with a relative error of 0.01, if ‘(l) = 106 [Rubinstein and Kroese, 2004]). [12] Rubinstein [1997] proposed Importance Sampling (IS) for rare events probability estimations. IS is a variance reduction technique in which sampling is performed using a probability distribution which increases the occurrence likelihood of rare events, and thus estimates their probabilities more efficiently. The major deficiency of employing an IS scheme is its unknown associated parameters, which are complex to calibrate [Rubinstein and Melamed, 1998]. An alternative approach for estimating the IS reference parameters [Rubinstein, 1999] is to minimize the Kullback – Leibler distance measure which quantifies the distance D [g(x), w(x)] between two density functions g (x) and w (x) as: Z D½gðxÞ; wðxÞ ¼

Z gðxÞLN½gðxÞ dx

gðxÞLN½wðxÞ dx

ð9Þ

where LN denotes the natural logarithm. [13] The minimization of the Kullback – Leibler distance as in equation (9) for estimating the IS reference parameters is particularly appealing as its solution can often be obtained analytically [Rubinstein, 1999]. [14] Minimization of the Kullback – Leibler distance measure within the context of IS yields the main CE algorithm

2 of 15

W09413


for rare events simulation estimation, which further forms the CE algorithm for combinatorial optimization. [15] To apply the CE method, the problem in hand needs to be transformed into an associated stochastic graph problem representation of one of two types: (1) a stochastic edge network (SEN): introducing randomness to the graph links; or (2) a stochastic node network (SNN): introducing randomness to the graph nodes. An example of a SNN representation is the max-cut problem (i.e., given a graph with a weight associated to each link, find the maximum cut which partitions the graph nodes into two sets); an example of a SEN representation is the traveling salesman problem (i.e., find the shortest tour length of an agent required to visit a predefined number of cities exactly once, given that all cities are directly connected and that their associated distance length are known). A detailed description of the associated stochastic graph problem construction and the SNN and SEN examples for the max-cut and traveling salesman problems, respectively, can be found in De Boer et al. [2005]. [16] Once an associated stochastic problem is formulated the CE iterates between two main steps: (1) generation of random sample data solutions, and (2) updating of the parameters of the associated stochastic problem on the basis of the sampled data in the direction of solution improvements. [17] As a water distribution system is intuitively modeled as a graph with the links representing the pipes and the nodes representing connections between pipes and hydraulic control elements, each random parameter of the associated stochastic problem can uniquely define a node or an edge of the system.

4. Multiobjective Optimization [18] Real-world problems often require the simultaneous optimization of multiple, possibly conflicting objectives. The most common goal of an analysis is to select the best trade-offs among competing objectives. Frequently, optimization problems are structured as single objective problems by lumping all goals into a single minimization or maximization framework. This type of analysis is beneficial for gaining insights to the considered problem. However, it does not allow for a clear trade-off among the employed objectives. Contrary to that, in a multiobjective formulation there is no single optimal solution, but sets of nondominated solutions which form the Pareto trade-off curve among all goals. [19] In recent years several methods have been developed which extend single objective evolutionary schemes to multiobjective algorithms. Three of the more utilized algorithms are the multiobjective genetic algorithm (MOGA) [Fonseca and Fleming, 1995], the nondominated sorting genetic algorithm II (NSGA II) [Deb et al., 2002], and the strength Pareto evolutionary algorithm II (SPEA II) [Zitzler et al., 2001].

5. Multiobjective CE for Water Distribution Systems Optimization [20] In this study the maximum pressure deficit at a demand node, as in Farmani et al. [2005], serves as an objective for system performance quantification which is

W09413

negatively correlated with the system’s design cost (i.e., a zero pressure deficit implies a full attainment of pressure design constraints which corresponds to maximum cost). The decision maker is provided with a trade-off curve between minimum system cost and maximum pressure deficit, and thus can quantitatively evaluate the cost of pressure constraints attenuation which implies a reduction in the system service to its consumers. Other objectives such as maximizing reliability versus design cost minimization can also be employed. [21] The extension of the CE method [Perelman and Ostfeld, 2007] for multiobjective water distribution systems optimization derives some of its features from the MOGA algorithm suggested by Fonseca and Fleming [1995]. [22] Fonseca and Fleming [1995] used nondominated ranking and selection to move solutions toward the Pareto front. The ranking method assembles the density information into a rank value, where an individual solution rank corresponds to how many individuals in a given population dominate it. All the nondominated individuals are assigned a rank value of one, whereas dominated solutions are penalized according to the population density of the corresponding region of the trade-off surface.Consider u(i) t as the ith solution at iteration t. Each u(i) t solution is assigned a rank R in the following manner: ð iÞ ð iÞ R ut ¼ 1 þ Ut

ð10Þ

(i) where U(i) t is the number of solutions which dominate ut in iteration t. [23] Fitness assignment is performed by: (1) sorting the population according to the assigned ranks, (2) assigning fitness to solutions by interpolating from the best rank to the worst, and (3) averaging the fitness of the solutions with thesame rank, so that all of them will be further sampled at the same rate. [24] The main limitation in current multiobjective evolutionary algorithms is in finding a suitable fitness assignment approach to search for the near-optimal and near-complete diverse Pareto front. Because of the inherent nature of the CE algorithm the fitness assignment procedure is not required at the developed optimization scheme. [25] The proposed multiobjective CE algorithm, based on the underling theory of the CE method (i.e., equations (7) – (9)), involves the following steps: [26] 1. Iteration counter: set an iteration counter t = 0. [27] 2. Initialization: choose an initial probability vector ^ p0 with components p0,i (i = 1,. . ., m) where p0,i is the probability of success of a diameter i (i.e., selection) at t = 0, and m is the total number of available diameters. [28] 3. Sample solutions: randomly generate N sample vectors Xi (i = 1,. . .,N), at iteration t using the probability vector ^ pt (i.e., generate N ‘‘zero – one’’ solution vectors each of size m, where a ‘‘one’’ implies a diameter selection and a ‘‘zero’’ otherwise). [29] 4. Performance assignment: assign each generated solution vector XI a vector of performance objective values S(Xi) (e.g., for the two objectives employed: the cost of design Xi and the maximum pressure deficit of the system corresponding to design Xi).

3 of 15


W09413

[30] 5. Rank assignment: assign each sampled vector Xi a rank R[S(Xi)]: R½SðXi Þ ¼ 1 þ

N X G S Xj

i ¼ 1; . . . ; N; i 6¼ j

ð11Þ

j¼1

1 if Xj dominates Xi G S Xj ¼ 0 otherwise

ð12Þ

[31] 6. Solution sorting: sort solutions according to their assigned ranks: R1 ½SðX1 Þ R2 ½SðX2 Þ . . . RN ½SðXN Þ

ð13Þ

where R1 [S (X1)] is the lowest rank value corresponding to design X1. The solution ensemble associated with the lowest rank values comprises the best Pareto front of the current iteration. ^ t: [32] 7. Updating of the probability vector p [33] 7.1. Elite sample: select a percentage (e.g., 1%) of the top best performance vectors X (denoted also as the Elite sample) corresponding to the top best (i.e., lowest) rank solutions. ^t to [34] 7.2. Updating: update the probability vector p ^t+1: p ^ ptþ1;i ¼

Ft;i rN

i ¼ 1; . . . ; m

ð14Þ

where: ^pt + 1,i = the ith component of the probability vector ^t + 1 at iteration t + 1; r = the Elite sample percentage (e.g., p 1%); and Ft,i = the number of times diameter i is selected at iteration t within the Elite sample [e.g., if N = 1000, r = 0.01, and diameter i is selected 6 times within the best 1% ranked performance vectors X, then: Ft,i = 6, rN = 10, and ^pt+1,i = 0.6]. [35] 7.3. Smoothing: to reduce the likelihood that the probability vector will get ‘‘stuck’’ on a ‘‘zero – one’’ solution which will prevent the exploration of new domains, a smoothing parameter a [a 2 (0, 1)] is employed: ^ptþ1;i

a^ptþ1;i þ ð1 aÞ^pt;i

8i ¼ 1; :::; m

ð15Þ

For example, the probability value of 0.6 of diameter i above at iteration t = 1 with ^ p0,i = 0.5, and a = 0.7, is smoothed to: ^p1,i = 0.7 0.6 + (1 0.7) 0.5 = 0.57. [36] 8. Check convergence criteria: assemble the Pareto front using the solutions of all previous iterations. Convergence is declared if for a predefined number of subsequent iterations no additional nondominated solutions are added to the assembled Pareto front. Other convergence conditions such as a minimum generational distance (explored further below) between Pareto fronts in subsequent iterations could also be employed [Hanne, 2007]. If convergence is attained, STOP and define all current nondominated solutions as the best approximated Pareto front; otherwise: t t + 1 and return to step (3). [37] Table 1 and Figure 1 summarize an illustrative example for implementing the proposed multiobjective CE algorithm. The system layout, data, and an illustration of a

W09413

possible cut solution [Perelman and Ostfeld, 2007] are shown in Figure 1. The distribution system is comprised of a gravitational one loading one loop network of five pipes which is to be constructed using one of nine possible candidate diameters for each pipe, thus the design space incorporates 95 = 59049 possible solutions. The CE parameters employed are a sample size of 450 (i.e., 10 times the number of decision variables 5 9 = 45), an Elite r percentage of 0.02 (i.e., the size of the Elite sample is 9), and a smoothing parameter a of 0.7. Table 1 presents the solutions obtained at the first four iterations. [38] Note that at iteration two the global optimal single objective least cost solution (i.e., the minimum cost corresponding to zero pressure deficit) was attained (i.e., the least cost result out of 95 = 59049 possible solutions found by enumeration). Following iteration four the Pareto front was not altered any further thus iteration four holds the best Pareto front. [39] The CE algorithm parameters as described above are the sample size N, the Elite sample percentage r, the smoothing parameter a, and the convergence criteria been a predefined number of iterations for which no additional nondominated solutions are added to the Pareto front. All four parameter values are generally to be set through model calibration. For most applications [Perelman and Ostfeld, 2007; Rubinstein, 1999; Rubinstein and Kroese, 2004] N is within 1 to 20 times the total number of nodes of the SNN problem formulation; 0.005 r 0.03; and 0.3 a 0.9. [40] It should be emphasized that the maximum pressure deficit is not a spatially distributed criterion, thus one can not tell from the analysis the number of households, for example, which will experience a reduction/shortage in flow rates. Inclusion of spatially distributed objectives is possible and is not expected to pose conceptual difficulties to the proposed methodology. [41] The algorithm performance capabilities, as well as the parameters selection and influence on the proposed scheme, are further explored in the sample applications section below.

6. Sample Applications [42] Two sample applications are presented in this section for demonstrating the proposed methodology capabilities: (1) the one loading gravitational expansion of the New York tunnels (NYT) water system [Schaake and Lai, 1969], and (2) an extension of the two loop network (TLN) of Alperovits and Shamir [1977] to multiple loadings, pumping, and storage [Salomons, 2001]. The multiple loadings, pumping, and storage example [Salomons, 2001] was solved using CE in Perelman and Ostfeld [2007] for minimizing cost alone. System costs for both examples are expressed in terms of present values (i.e., compound costs). 6.1. Example 1: New York Tunnels Water Distribution System [43] The problem of the New York tunnels (NYT) water system dates back to Schaake and Lai [1969]. The city was interested to increase its capacity for meeting future demands by adding parallel pipes to its existing 21 tunnels. Since the early linear programming model of Schaake and Lai [1969], the NYT system was the subject of extensive research for minimizing its expansion piping cost [Perelman

4 of 15

W09413


W09413

Figure 1. Illustrative example layout, SNN representation, and a cut solution. and Ostfeld, 2007; Eusuff and Lansey, 2003; Maier et al., 2003; Savic and Walters, 1997; Loganathan et al., 1995; Fujiwara and Khang, 1990; Morgan and Goulter, 1985; Quindry et al., 1981]. Farmani et al. [2005] compared the nondominated sorting genetic algorithm II (NSGA-II) [Deb et al., 2002] to the strength Pareto evolutionary algorithm II (SPEA-II) [Zitzler et al., 2001] for the NYT system and for the Hanoi [Fujiwara and Khang, 1990] water system. In this study the same objectives as in Farmani et al. [2005] of cost minimization versus minimization of the maximum pressure deficit, are employed. The proposed multiobjective CE algorithm is compared to NSGA-II. [44] The layout of the NYT system is shown in Figure 2. The system is subject to a one demand loading condition, and consists of 21 links and 20 demand nodes supplied by a single reservoir at a constant head of 91.4 m. All the 21 links are candidates for reinforcement by additional parallel pipes, with each pipe having one of 15 possible pipe diameters: 914, 1219, 1524,. . ., and 5182 mm, and a ‘‘zero’’ alternative of ‘‘do nothing’’ (i.e., a total of 16 possible options for each link). Thus the number of decision variables is: ndv = 336 [21 (number of links) 16 (15 candidate diameters for each link plus a ‘‘do nothing’’ option)], and the design space is 1621 = 1.934 1025. [45] All nodes are at zero elevation. The minimum pressure head requirements for the most downstream nodes 16 and 17 are 79.2 and 83.1 m, respectively. For all other nodes the minimum pressure head requirement is 77.7 m. The cost Ci ($) of installing a new parallel pipe of diameter di (mm) and length Li (m) is: Ci ¼ 6:537 102 d1:24 Li i

i ¼ 1; . . . ; np

ð16Þ

The Hazen –Williams roughness coefficient for the existing and considered new parallel pipes is assumed to be 100. The demands, links, length, and existing diameters data are as in Schaake and Lai [1969]. [46] Tables 2, 3, and 4 and Figures 3, 4, 5, 6, and 7 describe the analysis for the NYT system. Two types of analyses were performed: evaluation of the proposed multiobjective CE methodology capabilities, and comparison to NSGA-II. 6.1.1. Evaluation of the Proposed Multiobjective CE Methodology Capabilities [47] There are two basic distinct goals in multiobjective optimization: (1) to find solutions as close as possible to the true Pareto optimal front, and (2) to disclose solutions as diverse as possible. The first objective requires a search toward the Pareto optimal region, while the second entails a search along the Pareto optimal front. Solutions are sought to be close to the true Pareto optimal front while also well distributed. In some sense thus these two objectives are conflicting. [48] Various performance metrics were proposed for Pareto optimal fronts comparison. Of the more commonly utilized metrics are: diversity [Deb et al., 2000], attainment surface [Fonseca and Fleming, 1996], attainment surface sampling [Knowles and Corne, 2000], generational distance [Van Veldhuizen, 1999], and, distance and distribution [Ang et al., 2002]. Further detailed discussion of various other multiobjective optimization comparison metrics can be found in Sarker et al. [2002]. [49] To evaluate the proposed CE methodology capabilities, three of the above metrics were employed: generational distance [Van Veldhuizen, 1999], distance measure

5 of 15


W09413

W09413

Table 1. Results for the CE Multiobjective Illustrative Examplea Pareto Front Solutions at First Four Iterations

Pipe Diameters, mm Solution Number

1

2

Sampled Solutions at Iteration 1 1 356 305 2 254 356 3. 254 356 . .

.. .. .

.. .. .

.. .. .

450 356 203 Elite Solutions at Iteration 1 1 356 203 2 305 254 3 305 203 4 356 254 5 356 203 6 356 203 7 356 254 8 356 51 9 254 356 Elite Solutions at Iteration 2 1 305 254 2 305 254 3 356 203 4 305 254 5 305 254 6 356 254 7 305 203 8 356 203 9 254 203 Elite Solutions at Iteration 3 1 356 203 2 356 254 3 305 254 4 305 254 5 356 203 6 356 254 7 305 203 8 356 203 9 254 203 Elite Solutions at Iteration 4 1 356 254 2 305 203 3 356 203 4 305 254 5 305 254 6 356 254 7 305 254 8 305 254 9 356 254 Number of nondominated solutions

3

4

5

Cost ($103)

Maximum Pressure Deficit (m)

Rank

254 254 305 .

356 152 25 .

254 203 356 .

234 163 204 .

0.0 21.3 19.0 .

7 3 12 .

203

305

356

216

8.2

9

305 152 254 254 305 152 305 254 254

203 51 254 254 254 76 102 254 152

76 254 51 152 76 305 254 203 203

164 135 142 172 173 157 185 152 163

5.2 12.9 9.3 0.0 0.0 28.6 0.0 29.5 21.3

1 1 1 1 2 3 3 3 3

152 152 254 152 254 152 254 254 305

51 51 254 25 152 51 203 203 254

254 203 76 203 203 254 76 102 51

135 126 *155 123 153 145 136 149 142

12.9 23.9 0.0 25.0 5.0 5.1 21.2 8.8 21.1

1 1 1 1 1 1 2 2 2

254 203 152 152 254 152 305 254 254

254 152 25 51 254 25 254 254 254

76 203 254 203 51 254 51 152 102

*155 154 132 126 152 142 160 163 130

0.0 0.6 12.7 23.9 1.4 4.8 0.0 0.0 27.1

1 1 1 1 1 1 2 3 3

25 254 254 25 51 25 25 25 152 front

254 51 51 254 203 254 203 254 254

149 160 152 139 126 142 123 132 163

2.1 0.0 1.4 10.0 23.9 4.8 25.0 12.7 0.0

1 1 1 1 1 1 1 1 2

.. .. .

203 305 254 203 152 152 152 152 203 at Pareto

.. .. .

.. .. .

.. .. .

.. .. .

1

2

3

4

.. .. .

+ + + +

+

+ + + + + +

+

+ + + + + +

+ +

+ + + + + + + 4

7

7

9

a + = Solution present at Pareto front; *155 = global optimal single objective least cost solution (i.e., for zero pressure deficit) out of 95 = 59049 possible results.

[Ang et al., 2002] for assessing the proximity of individual solutions of a Pareto front to the best approximated Pareto front, and distribution measure [Ang et al., 2002] for evaluating the diversity of the solutions along the Pareto frontier. [50] Generational distance is the averaged normalized Euclidean distance of the solutions of a Pareto front to the nearest best approximated Pareto front solution. In Figure 3, for example, d1, d2, and d3 are the normalized Euclidean distances of each solution to the nearest best approximated Pareto front solution, and the generational distance is the sum of the normalized Euclidean distances d1, d2, and d3

divided by three. The distance measure is similar to the generational distance metric except that it is measuring individual distances rather than an overall average distance. A zero value distance measure indicates that a solution coincides with a best approximated Pareto front solution, and any value above zero indicates a deviation from the best Pareto front. The distribution measure is the normalized Euclidean distance between each two solutions of the Pareto front, with the boundary solutions measuring the gap to the boundary solutions of the best approximated Pareto front. In Figure 3, g2 and g3 represent distribution distances between solutions, and g1 and g4 the distribution distances between

6 of 15

W09413


Figure 2. Layout of the New York tunnels water distribution system.

the boundaries of the evaluated Pareto front to the approximated best Pareto frontier. Low distribution values, in particular, for the boundary solutions, characterize a Pareto front with a good distribution. [51] Table 2 describes the outcome of 50 runs of the proposed multiobjective CE methodology as applied to the NYT problem. The parameters used are a sample size of 20 ndv = 6720 (ndv = 336); r = 0.02, and a = 0.7. Initially, all diameters and the option of ‘‘do nothing’’ are assigned the same likelihood to be selected, hence the initial probability vector was set to be uniform with: p0,i = 1/nd = 1/16, 8i = 1,. . .,336. [52] Convergence was declared if at five subsequent iterations no additional nondominated solutions were added to the Pareto front. The number of five subsequent iterations was selected after testing the model for several runs and observing its number of nondominated solutions as a function of the number of iterations. Other stopping conditions such as a minimum generational distance between subsequent Pareto fronts could also be applied [Hanne, 2007]. [53] It can be seen from Table 2 that the average number of iterations to convergence was 29.1 with a standard deviation of 3.89; and that the average number of nondominated solutions at a Pareto front was 117.74 with a standard deviation of 3.32. The number of nondominated solutions at the approximated best Pareto front (i.e., the Pareto front utilized by merging all 50 Pareto front solutions) was 123. At a single run a Pareto front had on average 90.52 solutions with a standard deviation of 9.77 at the approximated best Pareto front, and an average generational distance to the approximated best Pareto front of 9.33 104 with a standard deviation of 3.77 104.

W09413

[54] Figure 4 presents plots of the best approximated Pareto front, and the Pareto front corresponding to the maximum generational distance of 2.14 103 (Table 2, run 22). It can be seen from Figure 4 that the two Pareto fronts are apart up to cost solutions of about 14 106 $ and almost coincide otherwise. Since run 22 corresponds to the maximum generational distance, the Pareto fronts of all other runs will lay in between the two plots shown in Figure 4. It should also be noted that for zero pressure deficit (i.e., minimizing cost alone) the best known solution of 38.64 106 $ [Maier et al., 2003] is attained, and that for zero cost (i.e., the ‘‘do nothing’’ option) the maximum pressure deficit is 47.6 m. [55] Figure 5 describes the distance measure [Ang et al., 2002] for run 22 and for run 50 (Table 2), corresponding to the maximum and minimum generational distances, respectively. It can be seen from Figure 5 that solutions differ mostly for the first twenty nondominated low cost solutions for run 22 (see also Figure 4), and have a zero distance (i.e., coincide with the best approximated Pareto front) for most of the solutions (i.e., for 110 solutions out of 123, see Table 2) for run 50. [56] Figure 6 presents distribution measure plots [Ang et al., 2002] for run 12 and for run 45 which correspond to the maximum and minimum accumulated distributions, respectively. Figure 7 describes the Pareto front of the maximum accumulated distribution (run 12) as compared to the best approximated Pareto front. As with the distance measure (Figure 5), it can be seen from Figure 6 that high distribution values occurred for about the first 20 nondominated solutions (see also Figure 7 up to solutions corresponding to costs of 10 106 $), and less for the others. 6.1.2. Comparison to NSGA-II [57] In this section the proposed multiobjective CE algorithm is compared to the nondominated sorting genetic algorithm II (NSGA-II) [Deb et al., 2002], which is one of the most employed multiobjective evolutionary algorithms. [58] NSGA-II employs a nondominated sorting approach using a selection operator, which creates a mating pool by combining the parent and offspring populations, and by selecting solutions with respect to fitness and spread. Generations are populated starting with the best nondominated front and succeeding until the specified population size is reached. If at the final stage there are more individuals in the nondominated front than the available space, the crowded distance-based niching strategy is invoked to choose which individuals of that front will enter into the next population. A detailed description of the NSGA-II scheme can be found in Deb et al. [2002], and thus is not repeated herein. Tables 3 and 4, and Figure 8 summarize the comparison of NSGA-II to the proposed multiobjective CE algorithm. [59] In Table 3, each of the CE run outcomes, as given in Table 2, is compared to NSGA-II for the same CE computational effort evaluated as the total number of function evaluations (i.e., system cost designs) required until convergence. For example, run number 1 converged after 168000 system design evaluations [i.e., 25 (iterations) 6720 (sample size of a single iteration)]. This was translated to 840 NSGA-II generations with an assumption of 200 strings at each population. The NSGA-II mutation probability was 0.075. Crossover was performed by ran-

7 of 15


W09413

Table 2. CE Results for the New York Tunnels Problem

Run Number

Number of Iterations to Convergence

Number of Nondominated Solutions

Number of Nondominated Solutions at Best Pareto Fronta

Euclidean Normalized Generational Distance to Best Pareto Frontb

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 Average Standard deviation

25 33 29 35 29 30 24 38 30 31 26 34 31 33 23 32 25 28 22 24 30 25 26 33 25 33 30 29 35 25 34 29 28 31 24 28 31 32 26 25 25 37 25 30 27 28 25 33 32 32 29.10 3.89

114 119 117 120 122 121 116 118 118 111 125 117 122 126 116 121 117 116 120 116 108 117 119 117 117 115 116 120 121 116 116 120 116 120 121 115 118 116 121 111 117 117 119 120 115 119 114 115 117 122 117.74 3.32

92 93 89 98 90 92 87 103 95 71 93 72 95 93 86 91 93 98 86 102 83 67 96 93 79 102 90 85 101 102 84 81 99 100 101 100 95 77 90 81 75 69 101 96 90 97 96 78 89 110 90.52 9.77

6.42E04 8.74E04 4.26E04 8.17E04 9.52E04 1.03E03 1.00E03 6.15E04 7.82E04 1.07E03 8.96E04 1.02E03 9.45E04 9.10E04 9.74E04 8.21E04 7.56E04 7.34E04 1.20E03 6.20E04 1.21E03 2.14E03 8.33E04 6.91E04 1.65E03 6.69E04 9.23E04 1.16E03 8.18E04 6.00E04 9.05E04 1.10E03 6.83E04 5.94E04 6.40E04 8.42E04 5.65E04 2.05E03 8.79E04 1.69E03 1.59E03 1.58E03 6.28E04 7.55E04 7.17E04 7.00E04 6.29E04 1.09E03 8.85E04 3.14E04 9.33E04 3.77E04

a

Out of 123 non-dominated solutions at best Pareto front. Euclidean normalized generational distance to best Pareto front [Van Veldhuizen, 1999].

b

8 of 15

W09413

W09413


W09413

Table 3. CE Comparison to NSGA-II – CE Convergence Criteria Multiobjective Cross Entropy

NSGA-II [Deb et al., 2002]

Number of Nondominated Solutions at Merged Pareto Front

Run Number

Number of Iterations to Convergence



Total

Multiobjective Cross Entropy



25 33 29 35 29 30 24 38 30 31 26 34 31 33 23 32 25 28 22 24 30 25 26 33 25 33 30 29 35 25 34 29 28 31 24 28 31 32 26 25 25 37 25 30 27 28 25 33 32 32 29.10 3.89

114 119 117 120 122 121 116 118 118 111 125 117 122 126 116 121 117 116 120 116 108 117 119 117 117 115 116 120 121 116 116 120 116 120 121 115 118 116 121 111 117 117 119 120 115 119 114 115 117 122 117.74 3.32

52 59 60 53 68 66 54 63 54 68 51 57 61 59 60 65 66 58 52 60 56 50 55 59 64 60 59 60 57 65 49 63 63 57 54 49 60 62 64 64 47 62 63 59 54 49 66 75 57 64 59.04 5.88

114 119 117 120 124 121 116 118 118 105 125 116 122 126 116 122 120 116 120 116 111 116 119 117 117 115 116 120 121 116 114 120 117 120 121 115 118 116 121 109 108 118 119 120 114 119 114 116 118 122 117.56 3.99

114 119 117 120 121 121 116 118 118 95 125 115 122 126 116 121 114 116 120 116 105 115 118 116 109 115 116 120 120 116 107 120 116 120 121 115 117 116 120 103 100 117 119 120 114 119 114 115 117 122 116.24 5.82

9 of 15

1

5

10

1

1 6

7 1 1 1 8

1

8

1

1

1 7 8 1 1

1

1 1 2 1.52 2.69

W09413


W09413

Table 4. CE Comparison to NSGA-II – NSGA-II Convergence Criteria NSGA-II [Deb et al., 2002]


Number of Nondominated Solutions at Merged Pareto Front

Run Number

Number of Generations to Convergence



Total


1043 1093 619 733 725 1406 819 581 863 815 945 832 671 828 1051 778 384 955 660 942 633 828 700 909 862 940 1011 894 747 946 824 1357 1314 870 1150 965 1078 698 643 1134 816 1096 883 564 946 875 899 1109 1107 884 888.50 204.28

62 74 50 59 61 62 60 58 50 50 68 68 56 53 56 64 43 61 68 63 53 64 59 58 60 56 63 61 55 54 53 56 54 60 64 53 56 56 51 56 55 65 58 49 67 62 67 75 61 63 59.00 6.41

114 119 109 115 118 123 116 102 118 112 125 114 116 124 116 117 88 116 120 117 104 117 116 116 117 114 116 120 117 116 112 124 119 119 121 115 118 116 118 112 117 117 119 108 115 119 115 115 117 121 115.78 5.90

114 119 109 115 118 123 116 102 119 112 125 114 116 124 117 119 88 118 120 117 104 118 117 116 114 114 121 122 117 118 113 122 119 119 120 115 119 115 118 106 106 110 119 108 117 119 115 116 117 120 115.58 6.30

10 of 15



114 119 109 115 118 122 115 102 118 112 124 114 116 124 116 116 88 114 120 116 98 117 116 115 112 113 116 120 117 116 112 113 119 119 117 115 118 105 118 91 98 106 119 108 114 119 114 115 117 119 113.76 7.42

1 1

1

1

1 3

4

1 6 1 1 1 2 1 5 2

2 1 9

3

1 10

15 8 4

3

1 1

1 1.82 3.04

W09413


W09413

Figure 3. Distribution and distance measures. Figure 5. Distance measures for the New York tunnels problem. domly selecting and mating pairs of strings from the best 50% strings. For this run, the number of nondominated solutions was 114 and 52 for the CE and NSGA-II algorithms, respectively. At the merged Pareto front the 114 CE solutions dominated all 52 NSGA-II results. This same procedure was repeated for the rest 49 runs. It can be seen from Table 3 that from the average 117.56 nondominated solutions at a single merged optimal Pareto front only 1.52 solutions were contributed by the NSGA-II algorithm. [60] In Table 4, 50 NSGA-II runs are summarized with a population of 200 strings in each generation and a mutation probability of 0.075. Convergence of an NSGA-II run was declared if at 50 subsequent generations no-additional nondominated solutions were added to the Pareto front. The value of 50 subsequent generations was selected after testing the model for several runs and observing its number of nondominated solutions as a function of the number of generations (i.e., likewise the five iterations selected for the CE algorithm). Similarly to Table 3 the same computational effort of each of the NSGA-II runs was imposed on the CE algorithm. For example, in run number 1 the total number of design evaluations was 208600 [i.e., 1043 (generations) 200 (strings at a single generation)]. This was translated to 31 iterations with a sample size of 6720 for the CE algorithm. For this run, the number of nondominated solutions using NSGA-II and the CE algorithms was 62 and 114, respectively. At the merged Pareto front, the 114 CE solutions dominated all 62 NSGA-II results. This same procedure was repeated for the rest 49 runs. It can be seen

Figure 4. Best Pareto and maximum generational distance fronts for the NYT problem.

from Table 4 that from the average 115.58 nondominated solutions at a single merged optimal Pareto front only 1.82 solutions were contributed by NSGA-II. [61] Figure 8 describes the Pareto fronts of the multiobjective CE algorithm and NSGA-II for 3360000 design evaluations, which is a substantial increase over the number of evaluations utilized with the adopted stopping conditions as in Tables 3 and 4. The CE parameters used were 100 iterations; a sample size of 100 ndv = 336000; r = 0.005; and a = 0.7. Note that r was reduced from 0.02 to 0.005 as the sampling size was increased substantially, thus a coefficient of 0.02 might be too high for generating good elite samples for probability updating. The NSGA-II parameters were a population size of 1000 strings, 3360 generations, and a mutation probability of 0.075. The merged optimal Pareto front yielded a total of 123 solutions of which 108 were CE results and 21 NSGA-II solutions (six CE and NSGA-II solutions were identical). The generational distances to the merged optimal Pareto front were 6.18 104 and 1.40 103 for the CE and NSGA-II algorithms, respectively. The generational distance between the CE and the NSGA-II Pareto fronts was 1.83 103.

Figure 6. Distribution measures for the New York tunnels problem.

11 of 15

W09413


Figure 7. Best Pareto and maximum accumulated distribution fronts for the NYT problem.

6.2. Example 2: Multiple Loadings, Pumping, and Storage [62] Example 2 incorporates multiple loadings, pumping, and storage, and thus is a general implementation of the proposed multiobjective CE algorithm. [63] The system layout is shown in Figure 9. The unit costs of the candidate pipe diameters are as in Alperovits and Shamir [1977]. The system is subjected to four loading conditions and to four variable energy tariffs which span a typical demand day: 00 – 06 (hours), 0.2 (base demand multiplier), 0.25 (energy tariff in $/kWh); 06 – 12, 0.8, 0.80; 12– 18, 1.2. 1.50; 18– 24, 0.6, 0.60. All pipes have a length of 1000 m, excluding pipe 9 whose length is 100 m. The Hazen – Williams coefficient for all pipes is 130, the tank diameter is 36 m (assumed constant), the pump efficiency 0.75 (assumed constant), the minimum required pressure head at all consumer nodes is 30 m, and the initial (and final) tank water level is constrained to 2 m at midnight. [64] The tank construction cost is assumed to be 40000 $/ 1000 m3, the cost of constructing the pump is as defined in equation (17), and the cost of the pump operation and maintenance is as given in equation (18). The relationships defined in equations (17) and (18) follow Kessler and Shamir [1991]. The reader is referred to Kessler and Shamir [1991] for further explanations on the underlying assumptions for the development of equations (17) and (18) which are employed herein. [65] The cost of the pump construction is: 0:325 PCC ¼ 194680 Q0:65 m Hm

W09413

provided by the pump at time period i; ECi ($/kWh) is the energy cost; and 7200 is an empirical coefficient value. [67] The water distribution system design cost is thus the sum of constructing the pipes, the pump (i.e., PCC), the storage, and the pump operational cost (i.e., POC). [68] There are a total of 530 decision variables corresponding to the number of available pipes and the operating power of the pump, resulting in a 530 dimensional random vector associated with the optimization problem. The random vector has an initial probability distribution of p0,i = 1/14 for i = 1,. . .,126 and p0,i = 1/101 for i = 127,. . .,530. The first 126 elements of the probability vector correspond to the 14 available pipe diameters and the rest probability elements correspond to the operating power of the pump during each time pattern of the simulation. [69] The CE parameters used were a sample size of 10 ndv = 5300, r = 0.01, and a = 0.7. Convergence was declared if at five subsequent iterations the number of nondominated solutions at the Pareto front remained unchanged. [70] Figure 10 describes the approximated best optimal Pareto front received using 30 runs and the Pareto front corresponding to the run with the maximum generational distance to the approximated best optimal Pareto front. The average number of iterations until convergence (standard deviation) was 34.03 (5.95); the average number of nondominated solutions at a single run (standard deviation) was 132.27 (26.49); the average Euclidean normalized generational distance of a single run to the approximated best Pareto front (standard deviation) was 2.73 103 (1.32 103); and the number of nondominated solutions at the approximated best Pareto front was 189. [71] The approximated best optimal Pareto front shown in Figure 10 presents a diverse and well spread set of solutions. The proximity of the Pareto front corresponding to the maximum generational distance to the approximated best optimal Pareto front established by merging the obtained solutions from all 30 runs, demonstrates the algorithm stability. For zero maximum pressure deficit the best cost solution obtained was 2.28 106 $ which is slightly higher than the single objective best cost solution of 2.204 106 $ received by Perelman and Ostfeld [2007].

7. Conclusions [72] Numerous studies were published over the last four decades mainly concentrating on the single objective opti-

ð17Þ

where Qm (m3/s) and Hm (m) are the maximum flow and head, respectively, provided by the pump during the simulation period, which determine the pump construction cost PCC ($). The values 194680, 0.65, and (0.325) are empirical coefficients. The pump is assumed to be of variable speed. [66] The cost of the pump operation and maintenance is: POC ¼ 7200

24 X

Qi Hi ECi

ð18Þ

i¼1

where POC ($) is the cost of operating the pump; Qi (m3/s) and Hi (m) are the flow and pressure head, respectively

Figure 8. Best Pareto fronts of the CE and NSGA-II algorithms for the NYT problem.

12 of 15

W09413


W09413

Figure 9. Layout and data for Example 2. mal design of water distribution systems, ranging from the early gradient type schemes [Alperovits and Shamir, 1977] to the more recently evolutionary applications of genetic algorithms [Savic and Walters, 1997] and ant colony [Maier et al., 2003]. The many number of publications are due to the challenging highly nonlinear and nonsmooth properties of the optimal design problem [Eiger et al., 1994], whose search space grows exponentially with system size [Savic and Walters, 1997; Maier et al., 2003]. [73] Extensions of the optimal design problem to multiobjective formulations started recently [Halhal et al., 1999; Babayan et al., 2005; Kapelan et al., 2003; Keedwell and Khu, 2003; Farmani et al., 2005]. [74] This manuscript extends the CE method for combinatorial optimization to multiobjective optimization by incorporating elements from multiobjective evolutionary algorithms, and in particular by using the rank of the generated elite solutions to update the CE probabilities instead of using fitness function values. This modification resulted in almost no changes to the original CE method and thus preserved the CE original qualities of fast convergence to optimal or near-optimal solutions [Rubinstein, 1999]. The developed methodology was demonstrated using two sample applications and compared to the NSGA-II algorithm [Deb et al., 2002]. [75] The results showed a high potential of receiving good solutions with a relatively small number of function evaluations and a small number of iterations until convergence. The methodology was found to be robust and reliable in that it produced a well spread solution front (i.e., small distribution measure) in the vicinity of the approximated optimal Pareto front (i.e., small distance measure) requiring the tuning of only three parameters: the sample size N, the elite sample size r, and the smoothing parameter a. In addition, the proposed CE method for multiobjective optimization demonstrated sound performances when its parameters were altered in a relatively wide range of values. In all trails the proposed methodology provided improved results compared to the NSGA-II algorithm. [76] As an outcome of this study, the following issues need further consideration and research efforts:

[77] 1. Exploring the usage of the proposed method for more complex water distribution systems and for the inclusion of additional objectives such as reliability and water quality considerations. [78] 2. The employment of the Hazen – Williams equation for pipe headloss computations at the explored examples follows the common engineering practice in water distribution systems analysis. Recent research [e.g., Koutsoyiannis, 2008] has pointed out deficiencies in using the Hazen – Williams formula and suggested an improved power-type headloss equation. The utilization of such equations should be explored in future water distribution systems research. [79] 3. To reduce the dimensionality and price of the optimization process, the incorporation of schemes which employ aggregation algorithms which construct ‘‘equivalent’’ reduced networks containing fewer nodes and links which match the hydraulics of the entire system, should be developed. [80] 4. Most of the current research, including this study, has focused on optimization refinement, rather that optimization synthesis/formulation. Recent optimization techniques tend to be computationally and algorithmically efficient at finding optimal solutions but less to fundamental understanding of the nature of the solution. From a practitioner perspective, model parameters (e.g., friction coefficients, nodal demands) are largely unknown with any

Figure 10. Best Pareto and maximum generational distance fronts for Example 2.

13 of 15


W09413

degree of certainty. The challenge facing practitioners is in making a decision in a framework of conflicting objectives. Development of tools for narrowing this research-practitioner gap should be a primary goal.

Acronyms and Notations Acronyms ASP CE IS MOGA NYT PDF SEN SNN TLN Notations CHWi Ci c (d, L) D [g(x),w(x)] d di E ECi F1 F2 f (h, hmin) G (V, E) g (x) Hi Hm h hmin IIi K (V, E) kij L LN Li ‘ (l) ^ ‘(l) m N nd ndv np P PCC POC ^0 p p0,

i

^t p

associated stochastic problem. Cross Entropy. Importance Sampling. multiobjective genetic algorithm. New York tunnels. probability density function. stochastic edge network. stochastic node network. two loop network. Hazen –Williams headloss coefficient of pipe i. the cost of installing a new pipe for the NYT and Hanoi problems, $. cost function. the Kullback – Leibler distance measure. vector of pipe diameters. diameter of pipe i. set of edges. pumping energy cost, $/kWh. objective function 1. objective function 2. pressure deficiency objective function. graph G. probability density function. pressure head provided by a pump at time period i, m. the maximum head provided by a pump, m. vector of nodal heads. vector of minimum required nodal heads. the ith Indicator function. directed incidence matrix of G. the (i, j) K (V, E) matrix component. vector of pipe lengths. natural logarithm. length of pipe i. estimator of l. estimator of ‘ (l). number of available diameters. number of sample vectors. number of candidate diameters. number of decision variables. number of pipes. probability. pump construction cost, ($). pump operational cost, ($). initial probability vector with components p0, i (i = 1,. . ., m). probability of success (selection) of the ith decision variable at t = 0. probability vector at iteration t.

W09413

Qi flow provided by a pump at time period i, m3/s. Qm maximum flow provided by the pump, m3/s. q vector of pipe flows. qnode vector of nodal demands. qi flow along pipe i. RT number of samples. (i) R(u(i) t ) rank of ut . RPi (di, CHWi, Li) ith pipe resistance. R[S(Xi)] rank of Xi. S(Xi) vector of performance objective values of sampled solution Xi. t iteration counter. number of solutions which dominate U(i) t u(i) t in iteration t. (i) ut the ith solution vector at iteration t. V set of nodes. w (x) probability density function. Xi sample vector (i = 1,. . .,N). Z random variable. Zi ith simulated sample from the probability density function (PDF) of Z. a smoothing parameter. Ft,i number of times node i is selected at iteration t within the Elite sample. G[S (Xj)] dominance function. l real number. Dh(q) vector of headlosses along pipes. Dhi(qi) ith component of Dh(q). r Elite sample percentage. [81] Acknowledgment. This study was supported by the Technion Grand Water Research Institute (GWRI).

References Alperovits, E., and U. Shamir (1977), Design of optimal water distribution systems, Water Resour. Res., 13(6), 885 – 900. Ang, K. H., G. Chong, and Y. Li (2002), Visualization technique for analyzing non-dominated set comparison, in Proceedings of the 4th Asia-Pacific Conference on Simulated Evolution and Learning, vol. 1, edited by L. Wang et al., pp. 36 – 40, Nanyang Technical University, Orchid Country Club, Singapore. Babayan, A., D. A. Savic, and G. A. Walters (2005), Multiobjective optimization for the least-cost design of water distribution systems under correlated uncertain parameters [CD-ROM], in Proceedings of the EWRI/ASCE World Water and Environmental Resources Congress, edited by R. Walton, American Society of Civil Engineers, Anchorage, Alaska. De Boer, P. T., D. P. Kroese, S. Mannor, and R. Y. Rubinstein (2005), A tutorial on the Cross-Entropy method, Ann. Oper. Res., 134(1), 19 – 67. Deb, K., A. Pratap, S. Agarwal, and T. Meyarivan (2000), A fast and elitist multi-objective genetic algorithm: NSGA-II, KanGAL Report 200001, Indian Institute of Technology, Kanpur, India. Deb, K., A. Pratap, S. Agarwal, and T. Meyarivan (2002), A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Trans. Evol. Comput., 6(2), 182 – 197. Duan, N., L. W. Mays, and K. E. Lansey (1990), Optimal reliability-based design of pumping and distribution systems, J. Hydraul. Eng., ASCE, 116(2), 249 – 268. Eiger, G., U. Shamir, and A. Ben-Tal (1994), Optimal design of water distribution networks, Water Resour. Res., 30(9), 2637 – 2646. Eusuff, M. M., and K. E. Lansey (2003), Optimization of water distribution network design using the shuffled frog leaping algorithm, J. Water Resour. Plan. Manage. Div., ASCE, 129(3), 210 – 225. Farmani, R., D. A. Savic, and G. A. Walters (2005), Evolutionary multiobjective optimization in water distribution network design, Eng. Optim., 37(2), 167 – 183.

14 of 15

W09413


Fonseca, C. M., and P. J. Fleming (1995), An overview of evolutionary algorithms in multiobjective optimization, Evol. Comput., 3(1), 1 – 16. Fonseca, C. M., and P. J. Fleming (1996), On the performance assessment and comparison of stochastic multiobjective optimizers, in Parallel Problem Solving From Nature-PPSN IV, Springer-Verlag Lecture Notes on Computer Science No. 1141, edited by H.-M. Voigt, pp. 584 – 593, Springer-Verlag, London, U.K. Fujiwara, O., and D. B. Khang (1990), A two-phase decomposition method for optimal design of looped water distribution networks, Water Resour. Res., 26(4), 539 – 549. Halhal, D., G. A. Walters, D. A. Savic, and D. Ouazar (1999), Scheduling of water distribution system rehabilitation using structured messy genetic algorithms, Evol. Comput., 7(3), 311 – 329. Hanne, T. (2007), A multiobjective evolutionary algorithm for approximating the efficient set, Eur. J. Oper. Res., 176, 1723 – 1734. Jacobs, P., and I. Goulter (1989), Optimization of redundancy in water distribution networks using graph theoretic principles, Eng. Optim., 15(1), 71 – 82. Kapelan, Z. S., D. A. Savic, and G. A. Walters (2003), Multiobjective sampling design for water distribution model calibration, J. Water Resour. Plan. Manage. Div., ASCE, 129(6), 466 – 479. Keedwell, E. C., and S. T. Khu (2003), More choices in water distribution system optimization, in Advances in Water Supply Management, Proceedings of Computers and Control in the Water Industry, London, edited by C. Maksimovic, D. Butler, and F. A. Memon, pp. 257 – 265, Taylor and Francis, London, U.K. Kessler, A., and U. Shamir (1991), Decomposition technique for optimal design of water supply networks, Eng. Optim., 17, 1 – 19. Knowles, J. D., and D. W. Corne (2000), Approximating the nondominated front using the Pareto archived evolution strategy, Evol. Comput., 8(2), 149 – 172. Koutsoyiannis, D. (2008), A power-law approximation of the turbulent flow friction factor useful for the design and simulation of urban water networks, Urban Water J., 5(2), 107 – 115. Kullback, S., and R. A. Leibler (1951), On information and sufficiency, Ann. Math. Stat., 22(1), 79 – 86. Loganathan, G. V., J. J. Greene, and T. J. Ahn (1995), Design heuristic for globally minimum cost water-distribution systems, J. Water Resour. Plan. Manage. Div., ASCE, 121(2), 182 – 192. Maier, H. R., A. R. Simpson, A. C. Zecchin, W. K. Foong, K. Y. Phang, H. Y. Seah, and C. L. Tan (2003), Ant colony optimization for design of water distribution systems, J. Water Resour. Plan. Manage. Div., ASCE, 129(3), 200 – 209. Morgan, G. R., and I. C. Goulter (1985), Optimal urban water distribution design, Water Resour. Res., 21(5), 642 – 652. Ostfeld, A., and U. Shamir (1996), Design of optimal reliable multiquality water supply systems, J. Water Resour. Plan. Manage. Div., ASCE, 122(5), 322 – 333. Perelman, L., and A. Ostfeld (2007), An adaptive heuristic cross entropy algorithm for optimal design of water distribution systems, Eng. Optim., 39(4), 413 – 428.

W09413

Prasad, T. D., and N.-S. Park (2004), Multiobjective genetic algorithms for design of water distribution networks, J. Water Resour. Plan. Manage. Div., ASCE, 130(1), 73 – 82. Quindry, G. E., E. D. Brill, and J. C. Liebman (1981), Optimization of looped water distribution systems, J. Environ. Eng., ASCE, 107(EE4), 665 – 679. Rubinstein, R. Y. (1997), Optimization of computer simulation models with rare events, Eur. J. Oper. Res., 99, 89 – 122. Rubinstein, R. Y. (1999), The cross-entropy method for combinatorial and continuous optimization, Methodol. Comput. Appl. Probab., 1(2), 127 – 190. Rubinstein, R. Y., and D. P. Kroese (2004), The Cross-Entropy Method: A Unified Approach to Combinatorial Optimization, Monte Carlo Simulation, and Machine Learning, Springer Science and Business Media Inc., New York. Rubinstein, R. Y., and B. Melamed (1998), Modern Simulation and Modeling, Wiley Series in Probability and Statistics, Wiley, New York. Salomons, E. (2001), Optimal design of water distribution systems facilities and operation, MS thesis, Technion, Haifa, Israel. Sarker, R., M. Mohammadian, and X. Yao (Eds.) (2002), Evolutionary Optimization, Kluwer Academic, New York. Savic, D. A., and G. A. Walters (1997), Genetic algorithms for least cost design of water distribution networks, J. Water Resour. Plan. Manage. Div., ASCE, 123(2), 67 – 77. Schaake, J. C., and D. Lai (1969), Linear programming and dynamic programming applications to water distribution network, Report No. 116, Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Mass. Su, Y. C., L. W. Mays, N. Duan, and K. E. Lansey (1987), Reliability-based optimization model for water distribution systems, J. Hydraul. Eng., ASCE, 114(12), 1539 – 1556. Van Veldhuizen, D. A. (1999), Multiobjective evolutionary algorithms: Classifications, analyses, and new innovations, Ph.D. thesis, Department of Electrical and Computer Engineering, Graduate School of Engineering, Air Force Institute of Technology, Wright-Patterson AFB, Ohio. Walski, T. M. (1993), Water distribution valve topology for reliability analysis, J. Reliab. Eng. Syst. Saf., 42(1), 21 – 27. Zitzler, E., M. Laumanns, and L. Thiele (2001), SPEA2: Improving the strength Pareto evolutionary algorithm for multiobjective optimization, Evolutionary Methods for Design, Optimisation and Control, in Proceedings of the EUROGEN2001 Conference, Athens, Greece, edited by K. C. Giannakoglou et al., pp. 95 – 100, CIMNE, Barcelona.

A. Ostfeld and L. Perelman, Faculty of Civil and Environmental Engineering, Technion-Israel Institute of Technology, Technion City, Haifa 32000, Israel. ([email protected]; [email protected]) E. Salomons, OptiWater, 6 Amikam Israel Street, Haifa 34385, Israel. ([email protected])

15 of 15