Water Distribution System Design Using MultiObjective
Particle Swarm Optimisation
Abstract
Application of the multiobjective particle swarm optimisation (MOPSO) algorithm to design of water distribution systems is described. An earlier MOPSO algorithm is augmented with (a) local search, (b) a modified strategy for assigning the leader, and (c) a modified mutation scheme. For one of the benchmark problems described in the literature, the effect of each of the above features on the algorithm performance is demonstrated. The augmented MOPSO algorithm (called MOPSO+) is applied to five benchmark problems, and in each case, it finds nondominated solutions not reported earlier. In addition, for the purpose of comparing Pareto fronts (sets of nondominated solutions) obtained by different algorithms, a new criterion is suggested, and its usefulness is pointed out with an example. Finally, some suggestions regarding future research directions are made.
1 Introduction
Design of water distribution systems (WDS) involves various types of challenging optimisation problems which have been investigated in several research papers to date [8][31]. An important subset of these problems is the twoobjective optimisation problem in which pipe diameters are to be determined with the objectives of minimising the total cost and maximising the network resilience. Recently, Wang et al. [30] have reported a systematic collection of twelve benchmark WDS design problems of various complexities to serve as an excellent starting point for researchers to try out new multiobjective optimisation algorithms and compare their performance with the results available in the literature. Furthermore, they have made the corresponding EPANET [27] source code and the bestknown Pareto front data available in the public domain [4].
Network  Search  

space size  
TLN (SP)  8  14  15 M  128  
NYT (MP)  21  16  90 M  627  
BLA (MP)  23  14  90 M  901  
HAN (MP)  34  6  90 M  575  
GOY (MP)  30  8  90 M  480  
PES (IP)  99  13  150 M  782 
Of the twelve benchmark problems presented in [30], six are summarised in Table 1. For a “small problem” (SP), the size of the search space is relatively small, e.g., for the TLN problem, it is =, and it is possible to compute the cost and resilience for each of the possible networks to obtain the true Pareto front. In larger problems, including medium and intermediate problems (marked as MP and IP, respectively, in the table), it is not possible to compute the objective function values for each possible network configuration, and a multiobjective evolutionary algorithm (MOEA) is employed to obtain the approximate Pareto front (i.e., set of nondominated solutions). In these cases, the true Pareto front is not known, but following [30], we will often use the term “Pareto front” (PF) to refer to the set of nondominated (ND) solutions obtained by the concerned algorithm.
Five MOEAs were used in [30] for each of the benchmark problems. The ND solutions obtained by the five algorithms were combined, and dominated solutions were removed to form the “bestknown” PF. The number of solutions () in the bestknown PF thus obtained is listed in Table 1 for each benchmark problem. The percentage of solutions contributed to for each problem by each of the five algorithms is given in Table 2. It can be observed that, except for TLN and HAN networks, no single algorithm was able to obtain all of the solutions in the bestknown PF.
Network  % contribution  

NSGAII  MOEA  NSGAII  AMALGAM  Borg  
TLN (SP)  100.0  83.1  83.1  98.7  84.4 
NYT (MP)  93.8  17.9  24.8  91.7  20.7 
BLA (MP)  77.0  30.3  26.3  90.1  28.3 
HAN (MP)  100.0  20.5  23.1  89.7  25.6 
GOY (MP)  43.3  3.0  58.2  85.1  22.4 
PES (IP)  38.1  27.0  10.7  18.6  23.3 
The computation effort involved in obtaining the bestknown PF can be described in terms of the total number of function evaluations by all five algorithms together since no single algorithm was found adequate in general. For the problems in the MP category, was restricted to 600,000 for each of the five algorithms [30], and 30 independent (randomly initiated) runs were performed for each benchmark problem. The total is therefore == million (90 M), as shown in Table 1.
The purpose of this paper is to present a modified multiobjective particle swarm optimisation algorithm (MOPSO+) for WDS design. In particular, we present results for some of the benchmark problems described in [30] and show that the proposed algorithm gives better solutions than the bestknown PFs of [30]. The paper is organised as follows. In Sec. 2, we present a new scheme for comparing PFs produced by two algorithms, to be used in subsequent sections to gauge the performance of the MOPSO+ algorithm. In Sec. 3, we describe the local search scheme implemented in MOPSO+. In Sec. 4, we describe the overall MOPSO+ algorithm and point out specifically the changes with respect to the original MOPSO algorithm presented in [5]. We compare the results obtained using MOPSO+ with the bestknown PFs for medium problems in Sec. 5. We then point out in Sec. 6 limitations of the local search scheme for larger problems, propose a suitable modification, and show its effectiveness for one of the benchmark problems in the IP category [30]. Finally, we present conclusion of this work in Sec. 7 and comment on related future work.
2 Comparison of Pareto fronts
In comparing the performance of an MOEA with another with respect to a given benchmark problem, we are mainly interested in (a) computation time, (b) the “quality” of the set of ND solutions. Computation time depends on several factors such as the hardware used, programming environment (compiled or interpreted), and whether parallelisation was employed. Because these factors would vary, a more objective metric, viz., the total number of function evaluations is used in [30] and also in this work.
To judge the quality of the ND set obtained by an algorithm, the following metrics are commonly used [5],[6] when the decision variables are continuous: (a) generational distance, a measure of how far the obtained ND solutions are from the true Pareto front, (b) spacing, a measure of the spread of solutions throughout the ND set, (c) error ratio, the percentage of the ND solutions which do not belong to the true Paretooptimal set.
In the context of the benchmark WDS problems considered here, the true Pareto front is not known. In any case, the above metrics are of limited value for the WDS problem because of the discrete nature of the decision variables. From a practical perspective, apart from cost and resilience, the decision maker (DM) may need to take into account other considerations for implementation [18], and it is important to make a large number of options available to the DM. A more meaningful metric is therefore the number of solutions in the bestknown PF which the concerned algorithm can contribute [30]. Based on this idea, we propose some new metrics for comparing the PFs obtained by two MOEAs for a given benchmark problem.
Consider a twoobjective problem in which the objective functions and are to be maximised and minimised, respectively. Fig. 1 shows an example with two sets of ND solutions referred to as PFA and PFB. The combined set of ND solutions (marked as “combined PF” in the figure) is obtained by combining solutions from PFA and PFB and then removing solutions dominated by any other solutions. The combined PF contains 14 solutions. PFA, which contains a total of 10 solutions ( to ), contributes 5 unique solutions (, , , , ), i.e., solutions which are not present in PFB. Similarly, PFB, which contains a total of 13 solutions ( to ), contributes 6 unique solutions (, , , , ) to the combined PF. The number of common solutions in the combined PF, i.e., solutions which are contributed by PFA as well as PFB, is 3 (/, /, /). Note that some of the solutions have been rejected from the original PFs, viz., 2 solutions (, ) from PFA and 4 solutions (, , ) from PFB.
We use the following nomenclature to describe the above scenario: =, =, =, = for PFA, and =, =, =, = for PFB. The superscripts t, u, c, r stand for total, unique, common, and rejected, respectively. Note that and must be equal by definition, and we will denote them by . Clearly, for this example, we cannot say that PFA is better than PFB (or vice versa) since each of them has contributed some unique solutions to the combined set.
Consider now the case =, =. In this case, PFB is certainly better than PFA because it contributes some unique solutions while PFA contributes none. In other words, we can do without the information provided by PFA since PFB provides all of the solutions which PFA would have contributed to the combined front ( of them), and additional solutions. We will use this concept to compare the PFs obtained by the proposed MOPSO+ algorithm with the bestknown PFs [4].
3 Local search implementation
Memetic algorithms (MA), which combine the exploration capability of an evolutionary algorithm with the exploitation capability of a local search method, have been extensively used for multiobjective optimisation [10][11]. Barlow and Tanyimboh [2] have presented a memetic algorithm for WDS optimisation. In their work, the genetic algorithm (GA) was used along with local improvement. The child population after every GA iterations was obtained with local and cultural improvement operators using the following procedure.

Choose a subset S of the current ND set.

Select one individual from S for local improvement. Define a scalar fitness function with linear weighting where the weights are obtained using an estimate of the gradient of the PF. For the selected individual,

Find the HookeJeeves pattern search direction using the above scalar fitness function.

Perform the cultural learning step by applying the same pattern search direction to a group of individuals in the current ND set.


Repeat (b) until a sufficient number of children are created.
The authors demonstrated that the introduction of local improvement led to improved convergence speed and solutions better than previously reported.
In this work, we use a simpler local search (LS) scheme based on the observations that (a) New ND solutions can be found in the neighbourhood of known ND solutions [15],[11], (b) Local search is effective in exploring the leastcrowded areas of the objective space (where a smaller number of ND solutions have been found so far) [22]. In the following, we illustrate the LS scheme used in this work.
Consider the twovariable, twoobjective test problem [13] in which
(1) 
are to be maximised, subject to the constraints specified in [13].
Suppose the current ND set contains three points , , in Fig. 2 (a). The corresponding objective function values are shown by , , in Fig. 2 (b). Our goal is to improve the ND set using local search. To this end, we generate four neighbouring points centred around each existing solution, as shown in Fig. 2 (a). The and positions of the neighbours are given by
(2) 
The corresponding neighbours in the objective space are shown in Fig. 2 (b). We now find the ND solutions from the entire set consisting of the old solutions and the newly generated solutions. The new ND set is made up of , , , in Fig. 2 (b). Note that the ND set has improved in two ways: (a) The number of solutions has increased from 3 to 4, (b) The ND points in the new set dominate the old ND points – in this example, all of the old ND points (, , ).
It is important to evaluate the neighbours of each solution in the ND set first and then update the ND set by removing dominated solutions [16]. For example, in Fig. 2 (a), suppose we first evaluate the neighbours of and obtain two new ND points and . If the ND set is updated at this stage, solution will get removed (since is dominated by ). As a consequence, the neighbours of will not get evaluated, and the ND solution will not be obtained. To avoid this situation in MOPSO+, we first prepare a list of all solutions whose neighbours are to be explored. Subsequently, we treat the members of this list one by one until the entire list is covered. The new solutions are then used to update the ND set.
For the WDS problems we are considering here, the decision variables are the pipe diameters. Each decision variable can take on values 1, 2, , (where is the number of possible diameters, see Table 1), with 1 corresponding to the smallest diameter, 2 to the next larger diameter, etc. For the BLA network, for example, there are 23 decision variables, each taking a value between 1 and 14. We consider two solutions to be neighbours if they differ by in one of the decision variables. In the BLA case, therefore, a general solution not lying on a decision space boundary would have = neighbours. If there are solutions in the current ND set, we need to perform function evaluations, each involving a call the the EPANET program for computation of cost and resilience. The results before and after local search are shown in Fig. 3 for the BLA network at an early stage of the MOEA. The improvement arising from local search can be clearly observed.
The above procedure is obviously expensive. One way to reduce the number of function evaluations is to keep track of points which have already been evaluated. For these points, it is not necessary to repeat function evaluation. However, comparison of a given point with a list of points (in the decision space) is also expensive if it is compared with each point in the list one by one. To make the comparison more efficient, we employ a tree structure to represent the ND set. An example is shown in Fig. 4. In this case, there are three decision variables , , . Nodes at level are labeled by the values of the solutions. The path traced in checking if the solution =, =, =, is present in the ND set is shown with a thick grey line.
4 MOPSO+ algorithm
In this paper, we demonstrate the usefulness of a hybrid algorithm using the MOPSO algorithm along with local search for the WDS design problem. Several variations of the MOPSO algorithm have been reported (see [28],[23] for a review). Here, we use the MOPSO algorithm proposed by Coello [5], with some modifications.
The velocity update formula in MOPSO is similar to that used in singleobjective PSO, and is given by
(3) 
where denotes the PSO iteration number, and are the position and velocity of the particle, respectively, and and are random numbers between 0 and 1. The algorithm parameters are (inertia weight), (cognitive learning factor), and (social learning factor). In this work, we have used =, =, =. It should be pointed out that, although the decision variables in the WSD benchmark problems considered here take on discrete values, they are treated as real numbers in the velocity and position update equations. In computing the particle fitness, each decision variable is converted to the nearest integer (which gives the pipe diameter index for the concerned pipe).
The MOPSO algorithm differs from the singleobjective PSO algorithm in the computation of (the personal best position of the particle so far) and (the position of the leader). In [5], is updated in every iteration by comparing its current position with the previous value of . If the current position dominates, it replaces . If it is nondominated with respect to , then one of them is selected randomly as the next .
In MOPSO [5], assignment of is made using the ND set stored in an external archive (repository). The objective space is divided into hypercubes, and each ND solution, depending on its position in the objective space, is assigned one of these hypercubes. For each particle, in each PSO iteration, a leader is selected from the archive giving preference to ND solutions which occupy lesscrowded hypercubes. Roulettewheel selection is used to first select a hypercube, and one of the ND solutions in that hypercube is picked randomly as the leader. This procedure helps to ensure that the ND solutions are well distributed in the objective space.
In addition, a mutation operator is used in [5] to enhance exploration of the search space in the beginning of the search. The mutation rate is made zero as the algorithm converges.
With this background, we now describe modifications made in the proposed MOPSO+ algorithm.

Leader selection: The leader selection process in MOPSO [5] is illustrated in Fig. 5 (a). In each PSO iteration, each particle is assigned one of the ND particles in the archive, preferring less crowded hypercubes. In MOPSO+, we continue to use the roulettewheel selection procedure of the MOPSO algorithm. However, to intensify exploration of the lesscrowded regions of the archive, we assign the same leader to all particles, as shown in Fig. 5 (b) and keep the same leader for iterations.

Mutation: In PSO, when the velocity and position update steps fail to generate new ND solutions, mutation can be useful [23]. In the context of the WDS benchmark problems, we have observed that there is an initial phase of MOPSO in which the ND set is improved relatively rapidly. However, beyond a certain point, the rate of generation of new solutions drops significantly. For this reason, different mutation schemes have been implemented in MOPSO+ (see Fig. 6).
In the “constant” option, the mutation probability remains constant (a low value such as 0.01). In the “pulse” option, the probability is made nonzero only for iterations in the early stages and zero otherwise. In the “periodic” option, the probability is made nonzero for iterations in every iterations, thus periodically encouraging enhanced exploration. The mutation process itself is common in the three cases and involves changing one of the decision variables of the particle randomly.

Local search: The local search operation has been described in Sec. 3. We will refer to that procedure as a “unit local search” (ULS) step. In MOPSO+, local search is implemented as follows.

As seen in Sec. 3, a ULS step can lead to some improvement in the ND set. If it is applied again on the new ND set, further improvement is possible [16]. For this purpose, MOPSO+ allows the ULS to be repeated times at a given PSO iteration. If, after some ULS steps, it is found that no further generation of new ND solutions is taking place, the LS step is discontinued.

As explained in Sec. 3, local search is expensive, and it is not practical to perform it in every PSO iteration. In MOPSO+, therefore, LS is performed periodically instead of every iteration. Furthermore, it was observed in the context of the benchmark WDS problems that LS is more effective in the early stages. Based on this observation, a twostage LS strategy is implemented. From iteration to , LS is performed every iterations, and after , it is performed every iterations.

5 Results for medium problems
In this section, we present results obtained with the MOPSO+ algorithm for the “medium” problems (MP) listed in Table 1. We first consider the HAN network to understand the effect of the various algorithm parameters. We then choose the set of parameters which gives the best performance for the HAN network and use it for the other three benchmark problems in the MP category, viz., the NYT, BLA, and GOY networks. We will refer to the ND sets reported in [5] as “UExeter PFs” and the ND sets given by MOPSO+ as the “MOPSO+ PFs.”
MOPSO+  UExeter (PFA)  MOPSO+ (PFB)  

options  
LS/Leader/M  10,000  575  534  1  41  90 M  750  748  215  2  74.6 M  533 
LS/Leader/M  10,000  575  537  4  38  90 M  726  714  181  12  71.3 M  533 
LS/Leader/M  10,000  575  552  88  23  90 M  629  581  117  48  59.6 M  464 
LS/Leader/M  10,000  575  536  13  39  90 M  721  703  180  18  69.1 M  523 
LS/Leader/M  10,000  575  535  13  40  90 M  748  729  207  19  71.0 M  522 
LS/Leader/M  10,000  575  558  187  17  90 M  561  452  81  109  40.0 M  371 
LS/Leader/M  20,000  575  553  92  22  90 M  660  577  116  83  80.0 M  461 
Table 3 shows the results (i.e., , , etc. as defined in Sec. 2) for the HAN network with = (number of particles). For the first six rows of the table, the number of PSO iterations = while for the last row, it is 20,000.
The MOPSO+ options mentioned in the first column have the following meaning.

LS: Local search is performed with a period of 100 between = and 5,000, and with a period of 1,000 thereafter. The parameter (see Sec. 3) is set to 50.

Leader: The new leader assignment scheme (Fig. 5 (b)) is used with =.

M: No mutation is performed.

M: Constant mutation (Fig. 6 (a)) is performed with =.

M: Pulse mutation (Fig. 6 (b)) is performed with =, =, starting with =.

M: Periodic mutation (Fig. 6 (c)) is performed with =, =, =, starting with =.

LS: Local search is not used.
In each case, 20 independent runs of MOPSO+ were performed, and the total number of function evaluations (using the EPANET program) over all independent runs was recorded. We can make the following observations from Table 3.

For each set of MOPSO+ options, a substantial number of new ND solutions (given by ) are found.

Comparing the first two rows, we find that the new leader assignment scheme used in MOPSO+ gives a larger . It also results in a smaller , which is desirable, as explained in Sec. 3.

Among the different mutation schemes (see rows 1, 3, 4, 5), the periodic scheme gives the best PF.

Local search plays a very important role (see rows 1 and 6) in MOPSO+. Without LS, we see that a large number of UExeter solutions (=) are missed out by MOPSO+. Also, there is a substantial reduction (from 215 to 81) in the number of new solutions found by MOPSO+ when LS is not used.

Although LS has improved the ND set, it has taken a larger number of function evaluations (=M with LS as against =M without LS, see rows 1 and 6). To make a fair comparison, we have run MOPSO+ without LS with a larger =M by doubling , the number of PSO iterations. In this situation (row 7), and are worse than row 1 which implies that improvement due to LS is not simply because of increased .
It is instructive to look at the LS process for one specific run for the first row of Table 3. At a given PSO iteration at an early stage (=), the number of accepted solutions (), the number of rejected solutions (), and the number of solutions in the updated ND set () are plotted versus , the LS iteration number, in Figs. 7 (a)(c). At the first LS iteration (=), the number of solutions in the ND set is about 550, which means that there are (where 34 is the number of decision variables), i.e., 37,400 neighbours if all of the 550 points are interior. Some of these points would lie on the boundaries of the decision space, and they would have less than neighbours each. Nevertheless, it is clear that the number of solutions to be evaluated at the first LS iteration is rather large. Of these, only 138 solutions (Fig. 7 (a)) get accepted in the ND set. About 25,000 (Fig. 7 (b)) are found to be dominated by other solutions in the ND set and are therefore rejected.
During the first LS iteration (i.e., the ULS step described earlier), a list of all solutions being evaluated is prepared using the tree format shown in Fig. 4. In the subsequent LS iterations, the neighbours of the updated ND solutions are compared against this list, and only those not found in the list are evaluated. Therefore, a much smaller number of solutions need to be tried in the subsequent LS iterations. For example, for =, only about 5,000 solutions need to be evaluated. This process continues until either iterations are completed or when an LS iteration fails to generate any new ND solutions. For the example shown in Fig. 7, the latter situation is reached (at =). With each LS iteration, the ND set gets improved, and in this specific example, the number of solutions in the ND set also increases (Fig. 7 (c)) with each LS iteration.
Fig. 8 (a) shows the variation of the size of the ND set (), and Fig. 8 (b) shows the total number of LS iterations tried () as the MOPSO+ algorithm proceeds. Note that, as mentioned earlier, the period of local search is low (100) up to = and high (1,000) thereafter, and that is the reason for the larger density of data points up to = in Fig. 8 (b). For the first LS step (at =), ==. Beyond this point, generation of ND solutions due to LS slows down, and after about 2,000 PSO iterations, we see that mostly remains zero (which means only one LS iteration was tried which failed to produce new ND points). This points to a side benefit of using LS: it provides a necessary (but not sufficient) condition for convergence of an MOEA. If LS yields additional ND points, we can say that the MOEA has not converged.
MOPSO+ options  UExeter (PFA)  MOPSO+ (PFB)  

200  20,000  30  627  599  20  28  90 M  648  639  60  9  150.7 M  579 
300  20,000  30  627  597  28  30  90 M  670  630  61  40  211.3 M  569 
100  60,000  10  627  596  25  31  90 M  655  635  64  20  76.7 M  571 
100  60,000  17  627  595  4  32  90 M  661  656  65  5  130.3 M  591 
Network  MOPSO+ options  UExeter (PFA)  MOPSO+ (PFB)  

NYT  100  60,000  17  627  595  4  32  90 M  661  656  65  5  130.3 M  591 
BLA  200  10,000  10  901  849  0  52  90 M  1045  1045  196  0  44.1 M  849 
HAN  200  10,000  20  575  534  1  41  90 M  750  748  215  2  74.6 M  533 
GOY  100  20,000  10  489  444  3  45  90 M  571  570  129  1  37.9 M  441 
We now use the MOPSO+ algorithm parameters which have been found to work well for the HAN problem (option LS/Leader/M in Table 3) for optimisation of the other medium problems described in [5], viz., the NYT, BLA, and GOY networks. The results for the NYT network are shown in Table 4. We observe that the values are comparable in the four cases listed in the table. However, for =, =, and =, is substantially lower than the other three cases.
Comparing rows 1 and 2 of the table, we conclude that increasing the number of particles is not necessarily beneficial, which was also pointed out in [30] for other MOEAs. A systematic study of the effect of and on the performance of MOPSO+ has not been carried out in this work. Instead, we present the results obtained after limited experimentation with these parameters.
Table 5 shows the results for the four MP benchmark problems described in [5]. In each case, a substantial number of new ND solutions () have been found by MOPSO+. Furthermore, almost all ND solutions in the UExeter PFs are covered by MOPSO+ (see the column in the table). The improvement in the PFs have been obtained with a smaller for three networks and with a larger for the NYT network. These results overall point to the usefulness of the MOPSO+ approach for the WDS design problem.
Figs. 9 (a)(f) and 10 (a)(f) show expanded plots of the UExeter and MOPSO+ PFs. It can be seen that some of the MOPSO+ solutions dominate some of the UExeter solutions.
sol. 1  sol. 2  sol. 3  sol. 4  sol. 5  sol. 6  sol. 7  

Resilience  0.57931  0.57943  0.57984  0.58060  0.58077  0.58113  0.58129 
Cost  0.12518  0.12522  0.12530  0.12535  0.12540  0.12547  0.12552 
254.0  254.0  254.0  254.0  254.0  254.0  254.0  
101.6  101.6  101.6  101.6  101.6  101.6  101.6  
203.2  203.2  203.2  203.2  203.2  203.2  203.2  
203.2  203.2  203.2  203.2  203.2  203.2  203.2  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
25.4  25.4  25.4  25.4  25.4  25.4  25.4  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
101.6  101.6  101.6  101.6  101.6  101.6  101.6  
101.6  76.2  101.6  76.2  76.2  76.2  76.2  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
76.2  101.6  101.6  76.2  76.2  101.6  101.6  
50.8  50.8  50.8  50.8  50.8  50.8  50.8  
50.8  50.8  50.8  50.8  50.8  50.8  50.8  
76.2  50.8  76.2  25.4  76.2  25.4  76.2  
25.4  25.4  25.4  25.4  25.4  25.4  25.4  
101.6  101.6  101.6  152.4  101.6  152.4  101.6  
25.4  25.4  25.4  25.4  25.4  25.4  25.4  
152.4  152.4  152.4  152.4  152.4  152.4  152.4  
101.6  152.4  101.6  101.6  152.4  101.6  152.4 
Even though many of the unique MOPSO+ solutions (i.e., solutions present in the MOPSO+ PFs but not in the UExeter PFs) are close in the objective space to the previously known UExeter solutions, they can differ substantially in the decision space. For example, some of the ND solutions obtained by MOPSO+ for the BLA network are shown in Table 6. Of these, solutions 2, 3, 4 are also present in the UExeter PF while solutions 1, 5, 6, 7 are new. Although the cost and resilience values for solutions 1 and 2 differ only in the fourth decimal place, we observe that these two solutions vary substantially in the decision space – four diameter values, viz., , , , , are different. This could have important implications in practice because of nontechnical issues such as land acquisition, and a wider choice of ND solutions may well open up previously unknown options for the DM. The metrics proposed in Sec. 2 therefore seem to be more attractive from a practical perspective than metrics such as average cost, lowest cost, standard deviation, or spread.
6 Results for an intermediate problem
As seen in Sec. 3, local search is an expensive step, and with a large number of decision variables (in our context, the number of pipes) involved in the networks of IP and LP categories in [30], the repeated ULS procedure described for the MP category (see Sec. 4) can become prohibitively expensive. In this section, we will consider one of the IPs, the PES network (see Table 1). Since there are 99 decision variables, local search around one solution involves = function evaluations. As we shall see, the ND set for this network is much larger than the MP networks, making even one unit local search step involving evaluation of all neighbours of all solutions in the ND set consume significant computation time.
In view of the above difficulty, the ND set is only partially explored during local search (see [15],[11]), and only one ULS step is performed. In particular, solutions from the current ND set are randomly selected, and all neighbours of the selected solutions are evaluated. Apart from this change, other algorithm details remain the same as before.
For the PES problem, 10 independent runs of MOPSO+ were performed with = and =. Periodic mutation was employed, with =, =, =, starting with =. The new scheme for finding the leader (see Sec. 4) is used, with =. Local search was performed with a period of 1,000 between = and 10,000, and with a period of 10,000 thereafter. was set to 200, which means only 200 solutions from the current ND set are randomly picked for any LS step.
Fig. 11 shows (the number of ND solutions) versus PSO iteration number for one independent run. A much larger (about 25,000) is obtained for this network. Note that, in contrast to Fig. 8 for the HAN network, has not saturated in this case, indicating that may increase further with additional PSO iterations.
UExeter (PFA)  MOPSO+ (PFB)  

1  782  307  292  475  150 M  21,536  12,161  12,146  9,375  269 M  15 
2  782  190  175  592  150 M  20,736  14,856  14,841  5,880  538 M  15 
3  782  167  152  615  150 M  22,949  20,081  20,066  2,868  808 M  15 
4  782  166  152  616  150 M  22,962  20,094  20,080  2,868  1,077 M  14 
5  782  39  22  743  150 M  23,267  23,168  23,151  99  1,348 M  17 
6  782  37  19  745  150 M  24,643  24,547  24,529  96  1,619 M  18 
7  782  36  19  746  150 M  25,607  25,510  25,493  97  1,889 M  17 
8  782  36  19  746  150 M  26,145  26,048  26,031  97  2,160 M  17 
9  782  35  18  747  150 M  26,488  26,399  26,382  89  2,431 M  17 
10  782  35  18  747  150 M  30,267  30,178  30,161  89  2,703 M  17 
Table 7 shows the effect of (the number of independent runs) on the performance of MOPSO+. With a single run, = new solutions are generated. However, = solutions from the UExeter PF are missed out by MOPSO+. The number of function evaluations for this run is 269 M out of which = M are due to the PSO part and 69 M are due to local search. This value of is already larger than the UExeter value of 150 M. As more independent runs are performed, increases proportionately. However, the number of additional ND solutions increases with , and the number of UExeter ND solutions () missed by MOPSO+ decreases, both factors being favourable. For ten runs, for MOPSO+ is about 18 times larger than UExeter. Considering the large number of new ND solutions (about 38 times larger than UExeter), the additional computational effort seems worthwhile.
Some of the ND solutions obtained by MOPSO+ are shown in Figs. 12 (a)(c), along with UExeter solutions. It can be seen that the density of solutions obtained by MOPSO+ is larger. Also, for some values of network resilience (see Fig. 12 (c)), the cost for some of the MOPSO+ solutions is lower by as much as 25 %.
7 Conclusions
In summary, a hybrid algorithm (MOPSO+) comprising a multiobjective particle swarm optimisation algorithm and local search is presented for water distribution system design. Four medium and one intermediate benchmark networks are considered, and for each of them, a significant number of new nondominated solutions have been obtained using the proposed algorithm, thus improving the bestknown Pareto fronts for these problems. For medium networks, the total number of function evaluations for MOPSO+ was of the same order as in [30]. For the intermediate network considered here, although the total was much larger (18 times) for MOPSO+, the large number of new solutions generated by MOPSO+ justifies the additional computational effort.
Based on the results presented here, MOPSO+ seems to be an attractive option for WDS design. Some important issues for future research in this area are listed below.

Performance of MOPSO+ for other networks in the IP and LP category needs to be investigated.

A unified set of algorithmic parameters for all problems of a particular complexity (MP, IP, LP) is desirable. More work is required to explore this possibility.

The LS mechanism introduced here can be applied to other multiobjective optimisation problems with discrete decision variables.
References
 [1] (2013) A mixed discretecontinuous variable multiobjective genetic algorithm for targeted implementation of nonpoint source pollution control practices. Water Resources Research 49 (12), pp. 8344–8356. Cited by: §1.
 [2] (2014) Multiobjective memetic algorithm applied to the optimisation of water distribution systems. Water resources management 28 (8), pp. 2229–2242. Cited by: §1, §3, 4.
 [3] (2016) Use of domain knowledge to increase the convergence rate of evolutionary algorithms for optimizing the cost and resilience of water distribution systems. Journal of Water Resources Planning and Management 142 (9), pp. 04016027. Cited by: §1, 4.
 [4] Centre for water systems. External Links: Link Cited by: §1, §2, Figure 10, Figure 9, Figure 12.
 [5] (2004) Handling multiple objectives with particle swarm optimization. IEEE Trans. Evol. Comput. 8 (3), pp. 256–279. Cited by: §1, §2, Figure 5, §4, §4, §4, §4, Table 5, §5, §5, §5, 1, 2, 3.
 [6] (2009) Multiobjective optimization using evolutionary algorithms. J. Wiley and Sons, Ltd., Chichester. Cited by: §2.
 [7] (2013) Combining two search paradigms for multiobjective optimization: twophase and pareto local search. In Hybrid Metaheuristics, pp. 97–117. Cited by: §3.
 [8] (2005) Tradeoff between total cost and reliability for anytown water distribution network. Journal of Water Resources Planning and Management 131 (3), pp. 161–171. Cited by: §1.
 [9] (2014) Queued pareto local search for multiobjective optimization. In International Conference on Parallel Problem Solving from Nature, pp. 589–599. Cited by: §3.
 [10] (2002) Hybrid evolutionary multiobjective optimization algorithms.. In HIS, pp. 163–172. Cited by: §3.
 [11] (2018) Manyobjective pareto local search. European Journal of Operational Research 271 (3), pp. 1001–1013. Cited by: §3, §3, §6.
 [12] (2008) Optimal design of water networks using a modified genetic algorithm with reduction in search space. Journal of Water Resources Planning and Management 134 (2), pp. 147–160. Cited by: §1.
 [13] (1996) Multiobjective optimization by means of the thermodynamical genetic algorithm. In International Conference on Parallel Problem Solving from Nature, pp. 504–512. Cited by: §3.
 [14] (2005) Memetic algorithms for multiobjective optimization: issues, methods and prospects. In Recent advances in memetic algorithms, pp. 313–352. Cited by: §3.
 [15] (2012) On dominancebased multiobjective local search: design, implementation and experimental analysis on scheduling and traveling salesman problems. Journal of Heuristics 18 (2), pp. 317–352. Cited by: §3, §3, §6.
 [16] (2016) Double archive pareto local search. In 2016 IEEE Symposium Series on Computational Intelligence (SSCI), pp. 1–7. Cited by: §3, §3, 1.
 [17] (2017) Chebyshevbased iterated local search for multiobjective optimization. In 2017 13th IEEE International Conference on Intelligent Computer Communication and Processing (ICCP), pp. 163–170. Cited by: §3.
 [18] (2010) Multiobjective particle swarm optimization applied to water distribution systems design: an approach with human interaction. Mathematical and Computer Modelling 52 (78), pp. 1219–1227. Cited by: §1, §2.
 [19] (2016) Nondominated sorting differential evolution algorithms for multiobjective optimization of water distribution systems. Journal of Water Resources Planning and Management 143 (4), pp. 04016082. Cited by: §1.
 [20] (2015) Efficiency of evolutionary algorithms in water network pipe sizing. Water resources management 29 (13), pp. 4817–4831. Cited by: §1.
 [21] (2014) A comparison of populationbased optimization techniques for water distribution system expansion and operation. Procedia Engineering 89, pp. 13–20. Cited by: §1.
 [22] (2012) Local search based hybrid particle swarm optimization algorithm for multiobjective optimization. Swarm and Evolutionary Computation 3, pp. 1–14. Cited by: §3, §3.
 [23] (2008) Multiobjective particles swarm optimization approaches. In Multiobjective optimization in computational intelligence: Theory and practice, pp. 20–42. Cited by: §4, 3.
 [24] (submitted to Sādhanā) Using external archive for improved performance in multiobjective optimization. Cited by: 1.
 [25] (2011) Multiobjective optimisation of water distribution systems design using metaheuristics. Ph.D. Thesis, Stellenbosch: University of Stellenbosch. Cited by: §1.
 [26] (2007) An efficient multiobjective optimization algorithm based on swarm intelligence for engineering design. Engineering Optimization 39 (1), pp. 49–68. Cited by: §1.
 [27] (2000) EPANET 2: Users Manual. Cited by: §1.
 [28] (2005) Improving PSObased multiobjective optimization using crowding, mutation and dominance. In International Conference on Evolutionary MultiCriterion Optimization, pp. 505–519. Cited by: §4.
 [29] (2009) Local search based evolutionary multiobjective optimization algorithm for constrained and unconstrained problems. In 2009 IEEE Congress on Evolutionary Computation, pp. 2919–2926. Cited by: §3.
 [30] (2014) Twoobjective design of benchmark problems of a water distribution system via MOEAs: towards the bestknown approximation of the true pareto front. Journal of Water Resources Planning and Management 141 (3), pp. 04014060. Cited by: Table 1, Table 2, §1, §1, §1, §1, §1, §2, §2, §5, §6, §7.
 [31] (2017) An adaptive convergencetrajectory controlled ant colony optimization algorithm with application to water distribution system design problems. IEEE Transactions on Evolutionary Computation 21 (5), pp. 773–791. Cited by: §1.