Self-adaptive Co-evolution differential evolution algorithm for optimal design of water distribution networks

By Kathleen Austin,2014-09-06 22:00
8 views 0
Self-adaptive Co-evolution differential evolution algorithm for optimal design of water distribution networks


    Self-adaptive Co-evolution differential evolution algorithm

    for optimal design of water distribution networks

    5 FAN Qinqin, CHEN Kuilin, YAN Xuefeng

    (Key Laboratory of Advanced Control and Optimization for Chemical Processes of Ministry of Education, East China University of Science and Technology, Shanghai, 200237) Abstract: In this

    paper a new DE algorithm, i.e. self-adaptive co-evolution differential evolution (SACDE), is

    proposed to improve the optimization performance. SACDE employs the idea of

    10 co-evolution, which means that population and control parameters are coded together. While individuals are evolving, control parameters are changing in coordinately. The automatically updated parameters can adapt to different optimization problem and thus increase the robustness of the algorithm. Simulation results show that the proposed algorithm is better than classic DE algorithm. In comparison with other self-adaptive DE algorithm, the proposed algorithm is better in some conditions.

    15 Finally, the proposed algorithm is employed to optimize a discrete complex constrained optimization problemwater distribution network problem. SACDE achieves better results in benchmark water distribution network problems than previous methods. It demonstrates that SACDE not only can handle discrete constrained optimization problem but also shows its robustness due to the exponentially increasing complexity of different water distribution networks.

    20 Key words: Differential evolution algorithm;parameter co-evolution; water distribution networks

    0 Introduction

    Optimization problems raise many concerns in the reality, because industry problems that have been modeled are converted to find the optimal value of object function. The research of

    25 optimization problem weights much in the field of process control and process system engineering. Many engineering optimization problems contain multiple optimal solutions. Typically, the global optimal solutions are desired, rather than local optimal solution. In early years, deterministic optimization algorithms were mostly used such gradient descent method. However, deterministic optimization algorithms have many requirements for optimization problem and cannot adapt to

    30 different kinds of optimization problems.

    In recent years, evolutionary algorithms (EA) are popular for find the optimal solutions for

    non-linear multimodal problems. Inspired by the basic idea of sharing information among individuals of animal groups, E. Eberhart etc. designed the Particle Swarm Optimization (PSO)

    [1]algorithm . Based on the theory of biology evolution, Holland etc. proposed the Genetic

    [2] 35 Algorithm (GA), which imitates the natural selection and heredity in the field of biology. Inspired by the information sharing among ants in searching food, A. Colorni etc. proposed the

    [3]Ant Colony Optimization (ACO) algorithm , which demonstrates the information positive

    [4] feedback phenomenon. N. Metropolis proposed Simulated Annealing (SA) algorithm, which

    simulates the annealing process of materials.

    40 Differential evolution (DE) is a kind of reliable and versatile evolutionary optimization

    [5]algorithm proposed by Storn and Price . DE, unlike other EAs, generates offspring by perturbing

    the solutions with a scaled difference of two randomly selected population vectors, instead of

     recombining the solutions under conditions imposed by a probabilistic scheme. In addition, DE

    Foundations: National Natural Science Foundation of China (No.21176073); Doctoral Fund of Ministry of Education of China (No.20090074110005); Program for New Century Excellent Talents in University (No.NCET-09-0346); Shu Guang project (No.09SG29); the Fundamental Research Funds for the Central Universities.

    Brief author introduction:FAN Qinqin(1987), Male, PhD

    Correspondance author: YAN Xuefeng(1972), Male, Professor. The research interests include complex chemical process modeling, optimizing and controlling, and intelligent information processing. E-mail:

    - 1 -


    employs a one-to-one spawning logic which allows replacement of an individual only if the

    45 offspring outperforms its corresponding parent. DE is proved to be the most powerful optimization algorithm in 1996 IEEE evolution optimization contest. Although DE is a little slower than the other two deterministic optimization algorithms in terms of speed, those tow deterministic optimization algorithms have a very narrow application domain. The structure of DE is very

    [5]simple . DE only has two control parametersscale factor F and crossover possibility CR.

    0 These two parameters determine the performance of DE to a large extent. Due to its simple 5

    structure and fine convergence, DE is widely used in solving engineering optimization problem.

    However, DE shares some similar problems that other intelligent algorithm that also has, for instance, premature and unfavorable convergence. It is difficult to choose differential evolution algorithm’s control parameters, which plays a critical role in DE’s optimization performance. To

    55 improve the performance of DE, the algorithms controlling parameters scale factor F and

    ility CR, are modified, by evolving coordinately with the algorithm evolution crossover possib


1 Basic Differential Evolution Algorithm

1.1 Operation of basic DE algorithm

    60 In this section the basic operation of differential evolution is introduced. Necessary notations and terminologies are also described in order to facilitate the explanation of proposed SACDE algorithm later.

    Generally speaking, differential evolution algorithm follows the basic procedure of an

    [6] evolutionary algorithm . 0 x (i 1, 2, ..., NP) 65 Initiation: The initial population is randomly generated according to a i ( L) (U ) (j=1,2,…D), where D is the dimension of the problem andx x xuniform distributionj j j ( L ) (U ) NP is the population size. is the low boundary of the domain;is the up boundary of thex x j j

    domain. G 1 Mutation: At each generation G, a mutation vector v is created by following ways: i G G G G 1 x F ( x x) v 70 i r1 r 2 r 3

     Where r1, r2 and r3 are distinct integers randomly selected from the set {1,2,NP}, which G G G x x requires that NP>4. is the difference vector to mutate the parent vector x . F is the r 2 r 3 r1

    scale factor which usually ranges on the intervals [0, 2]. In class DE, F is fixed.

    Crossover: In order to increase the diversity, crossover operation is introduced. The trial

     75 vector is :

    G 1 G 1 G 1 G 1 , ..., u) (uu , u i 2i 1i Di ; G 1 if (randb( j) ? CR) or j rnbr(i) v ;jiG 1 u? ji Gx if (randb( j) ( CR) or j ? rnbr(i) ?ji ?

    Where randb(j) is a uniform random number on the interval [a, b] and independently generated for each j, rnbr(i) is an integer randomly chosen from 1 to D and newly generated for

     80 each i, and the crossover possibility in classic DE is a fixed parameter. CR ([0,1]

    Selection: The selection operation selects the better one from the parent vector from parent

    - 2 -


    G G 1 vector according to their fitness value. x and trial vector u i i

    It should be noted that some components of the trial vector may violate the predefined

    boundary constraints. Possible schemes include resetting schemes and penalty schemes. A simple

    5 solution is to randomly generate a trial vector form the domain to replace the origin one. 8

    1.2 Other variations of DE algorithm

    The basic operation of DE algorithm is stated above. In the practical application of DE, some

    variations of DE algorithm are developed to enhance the performance of DE. The notation

    DE/x/y/z is used to identify different variations of DE, where x defines the type of mutation vector, 90 for example random or ‘best’, y means the number of differential vectors, and z indicates the

    way of crossover. The basic DE algorithm can be described as DE/rand/1/bin.

    [7] Other famous variations of DE are stated below:

     (1) DE/best/1 v x F ( x x) i best r 2 r 3

    95 (2) DE/best/2 v x F ( x x) F ( x x) i best r 2 r 3 r 4 r 5 (3) DE/rand/2

     v x F ( x x) F ( x x) i r1 r 2 r 3 r 4 r 5

     (4) DE/current-to-best/1

    100 v x F ( x x) F ( x x) i i best i r 2 r 3

     2 Self-adaptive co-evolution differential evolution algorithm

     2.1 Basic ideas and flow of the algorithm The study shows that the scale factor F and crossover possibility CR are two important parameters for the performance of DE. Therefore, many researchers put their emphasis on the

    study of self-adaptive selection of F and CR. Due to the stochastic searching property of DE, total 105

     random or deterministic strategies cannot obtain the optimal parameters for the DE algorithm every time. Therefore, self-adaptive co-evolution differential evolution algorithm (SACDE) is

     proposed to update the parameters self-adaptively by evaluating the weight of the parameters. In the proposed algorithm, the control parameter scale factor and crossover factor are

     designed as symbiotic individuals of its original individuals. Each original individual has its own110 0 0 is the original individuals, which contains the solutions, symbiotic individuals. S S is 1 2 while the symbiotic individuals which contains the control parameterscale factor and crossover

     possibility. By individuals’ differential searching and accessing the weight of control parameter, the proposed algorithm obtains the weighted value of control parameter. Control parameters

    adaptive adjustment is realized by regenerating control parameter according to normal distribution. 115

    Thus, more proper control parameter for algorithm is obtained to improve the optimization


    Tab. 1 Original population and symbiotic population GG S S 12

     GG G F CRx 11 1

    - 3 -


     G G G F CRx 22 2

    ??? ???

     G G F CR x NPG NP NP Flow of algorithm: 0 0 120 (1) Initialization: Initialize origin population matrix S and parameter matrix S , 2 1

    determine number of individuals of population NP and number of iteration G. The boundary of m ( L) (U ) .the origin population is set as:x xxj j j GGG , are (2) Evolution of origin population S : Three mutually distinct individuals , , xxx1 r r r 1 2 3 G 1 v is generated by the pseudo-randomly selected from the population. A provisional offspring i

    following mutation: 125

    G1 G GG G v x F ? ( x x)r i r r i 2 3 1 GG G F is a scale factor which controls the length of exploration vectorWhere . Thex x i rr 2 3 GG effective value for is F ( (0,1) in most cases. F i i G1 G1 G 1 , then is selected randomly Boundary management: if v low or v( high v ij i ij i ij

    from population. 130 G 1 (3) Crossover operation: Each gene of individual v is exchanged with corresponding it G G 1 xgene of . The final offspring u is generated as following equation. it it G 1 G 1 G 1 G 1 , ..., u) (uu , u Di i 2i 1i G G x, rand(0,1)>CR,G 1 it i u t 1,2, ,n.? it G 1 otherwise, v it ?

    135 Where rand(0,1) is a random number between 0 and 1; is the index of gene under t G CR( (0,1) examination; is the crossover probability. i G 1 (4) Selection of the origin population: According the greed principle, u’s fitness value is i G x compared with ’s to decide which one will be selected into the next generation. i G 1 G 1 G u,f (u) ? f ( x)i i i G+1 ; Selecting operation:x ?i Gotherwise x, i

    140 (5) Control parameter evolution:

     In the standard DE and many other improved DE algorithms, the controlling parameter scale factor and crossover possibility are constant, who do not change throughout the whole evolution process. In the proposed algorithm, the scale factor and crossover possibility are evolving coordinately with the solution in each generation.

    145 Control parameters weight calculation: the weight of scale factor is calculated by the ration

    of the difference between the value of objective function upon one feasible solution in the

    population to and the maximum objective function value to the sum of such differences. So is the

    weight of crossover possibility.

    - 4 -


    G+1 ) f f ( xmaxi G ;Weight G NPFi G+1 ) f f ( x;maxi i 1 G+1 f ( x) f maxi G ;150 Weight G NPCRi G+1 f ( x) f ;maxi i 1

    Control parameters’ weighted average value calculation:

     NP;;weight Fi i i 1 G G G F ;(Weight ? F ) NPG G G CR) ?CR(Weight ;i weight CRii 1 Control parameters evolution: the next generation parameters are generated from the normal

    distribution. 155 G1 G F N F( , );i weight ; G1 G CR N CR(, ); iweight

     1.2 (G / G) Where ; N is the normal distribution function. m G 1 G 1 F 1F 1(i i Setting of control parameters’ boundary: ? G 1 G 1 F 0F 0?? i i G 1 G 1 CR 1CR( 1 ?; i i ;160? G 1 G 1 CR 0CR 0??; i i

    (6) Repeat (2) and (6) until the number of generation is bigger than the maximum allowable

    iteration number G. m 2.2 Simulation test of benchmark functions To test the performance of SACDE based on control parameter normal distribution, 20

    benchmark functions (shown in Tab.2), which are commonly used in literatures both at home and 165

     abroad, are employed. The results are compared with a famous self-adaptive differential evolution [8] [8]algorithm JADE and a traditional differential evolution algorithm DE/rand/1. The population and maximum iteration number can be seen in Tab.3. F=0.5 and CR=0.9 [5](recommended by Storn and Price ) are set for traditional DE. Each function is run

    independently for 50 times by SACDE, JADE and DE/rand/1. 170

    Tab.2 Benchmark functions

    Name Benchmark function Domain D D2Sphere [100,100] f ( x) ; x1 i i D DDf( x) | x| | x|Schwefel 2.21 [10,10] )2 i i i i 2D i DSchwefel 1.2 [100,100] f( x) (x )3 j i j

    Df( x) max{| x|,1 ? i ? D} Schwefel 2.21 [100,100] 4 i

    - 5 -


     D1 D2 2 2Rosenbrock f( x) [100( x x) (x1) ][30, 30] 5 i 1 i i i D 2DStep f ( x) ; (| x 0.5 |)[100,100] 6 i i D D4 Noisy Quartic [1.28,1.28] f ( x) ; ix rand[0,1)7 i i D Df( x) xsin | x| D;418.98288727243369Schwefel 2.26 [500, 500] 08 i i i D D2 Rastrigin [5.12, 5.12] f ( x) ; [x10 cos(2 x ) 10]9 i i i D 1 2 ( x) 20 exp(0.2 x)f 10 i i D D Ackley [32, 32] D 1 exp( cos(2 x )) 20 e ii D D D 1 x2 i D x cos( ) 1 f ( x) Griewank [600, 600] 11 )i i i 400 i D 1 2 2 2 ??f( x) {10 sin ( y) ( y1) 1 10 sin ( y) 13 i i i 1 ? ; D i 1 D 2 ( y1) } u( x,10,100, 4) D i i 1 x 1 Penalized Di ;[50, 50] y 1 i function1 4 m k ( x a)x ( a i i ? u( x , a, k , m) 0 a ? x ? a? i i ?m k ( x a)x a? i

     n-1 2 2 2 ; ?? f ( x) 0.1{sin 3 x ; ( x -1)1 sin (3 x ) ; 14 1 i i 1 ? i 1 n 2 ?? x -11 sin (2 x )} ; u( x , 5,100, 4) ; n n i ? Penalized Di 1 [50, 50] Function2 m k ( x a)x ( a i i ? u( x , a, k , m) 0 a ? x ? a? ;i i ?m k ( x a)x a? ;i 5.1 5 1 4 2 [5,10]?[0,15] Branin f ( x) ( x x x 6) 10(1 ) cos x 10 15 2 1 11 2 4 8 2 2 (19 14x 13xf ( x) (1 ( x x 1)15 1 2 1 1 2 14x 6x x 3x)) 22 1 2 2 Goldstein-price [2, 2] 2 2 * (30 (2 x 3x )(18 32 x 12 x 48x1 2 1 1 2

    2 36x x 27 x))1 2 2 4 3 23f( x) cexp[a( x p) ]Hartman3 [0,1] 16 i ij j ij i 1 j 1 4 6 26f( x) cexp[a( x p) ]Hartman6 [0,1] 17 i ij j ij i 1 j 1

    - 6 -


     5 4T 1[0,10] Shekel5 f( x) [(x a)(x a) c]18 i i i i 1 7 4T 1Shekel7 [0,10] f( x) [(x a)(x a) c]19 i i i i 1 10 4T 1Shekel10 [0,10] f( x) [(x a)(x a) c]20 i i i i 1

    Benchmark functions Branin, Shekel5, Shekel10, Hartman3, Hartman6 are low dimension


     Tab.3 Parameters of test functions Literatures SACDE Global Dimension optimum NP ? G NP ? G m m Spere 100 ?1500 50 ?1000 30 0

     Schwefel 2.22 100 ? 2000 50 ?1000 30 0

     Schwefel 1.2 100 ? 5000 50 ?1000 30 0

     Schwefel 2.21 100 ? 5000 50 ?1000 30 0

     Rosenbrock 100 ? 20000 100 ? 20000 30 0

     Step 100 ?1500 50 ?1000 30 0

     Noisy quartic 100 ? 3000 50 ? 6000 30 0

     Schwefel 2.26 100 ? 9000 50 ?1000 30 0

     Rastrigin 100 ? 5000 50 ?1000 30 0

     Ackley 100 ? 2000 50 ?1000 30 0

     Griewank 100 ? 3000 50 ?1000 30 0

     Penalized function1 100 ?1500 50 ? 3000 30 0

     Penalized Function2 100 ?1500 50 ? 3000 30 0

     Branin 30 ? 200 30 ? 200 0.397887 2

     Goldstein-price 30 ? 200 30 ? 200 3.00000 2

     Hartman3 30 ? 200 30 ? 200 3 3.86278

     Hartman6 30 ? 200 30 ? 200 3.32237 6

     Shekel5 30 ? 200 30 ? 200 10.1532 4

     Shekel7 30 ? 200 30 ? 200 10.4029 4

     Shekel10 30 ? 200 30 ? 200 10.5364 4


    Tab.4 Comparison of optimization results JADE without SACDE JADE with archive DE/rand/1 archive

    - 7 -


    Spere 0(0) 1.8E-60(8.4E-60) 1.3E-54(9.2E-54) 9.8E-14(8.4E-14)

    Schwefel 2.22 0(0) 1.8EE-25(8.8E-25) 3.9E-22(2.7E-21) 1.6E-9(1.1E-9)

    Schwefel 1.2 0(0) 5.7E-61(2.7E-60) 6.0E-87(1.9E-86) 6.6E-11(8.8E-11)

    Schwefel 2.21 0(0) 8.2E-24(4E-23) 4.3E-66(1.2E-65) 0.42(1.1)

    Rosenbrock 0(0) 8E-2(5.6E-1) 0.32(0.1) 8.0E-2(5.6E-1)

    Step 0(0) 0(0) 0(0) 0(0)

    Noisy quartic 1.8E-3(8.0E-4) 6.4E-4(2.5E-4) 6.8E-4(2.5E-4) 4.7E-3(1.2E-3)

    Schwefel 2.26 0(0) 0(0) 7.1(28) 57(76)

    Rastrigin 0(0) 0(0) 0(0) 71(21)

    Ackley 8.8E-16(0) 4.4E-15(0) 4.4E-15(0) 9.7E-11(5.0E-11)

    Griewank 0(0) 0(0) 2.0E-4(1.4E-3) 0(0)

    Penalized1 1.6E-32(5.5E-48) 1.6E-32(5.5E-48) 1.6E-32(5.5E-48) 1.1E-14(1.0E-10)

    Penalized2 1.3E-32(1.2E-47) 1.4E-32(1.1E-47) 1.4E-32(1.1E-47) 7.5E-14(4.8E-14)

    Branin 0.397887(3.5E-11) 0.39788(0) 0.39788(0) 0.39788(0) Goldstein-price 3.00000(1.6E-15) 3.00000(1.1E-15) 3.00000(5.4E-16) 3.00000(6.4E-16)

    Hartman3 -3.86278(2.8E-15) -3.86278(0) -3.86247(1.3E-16) -3.86278(0)

    Hartman6 -3.31488(0.026) -3.31044(3.6E-2) -3.30806(3.9E-2) -3.24608(5.7E-2)

    Shekel5 -9.6028(1.909) -10.1532(9.4E-14) -9.5469(1.6) -10.1532(4.2E-16)

    Shekel7 -10.2502(0.0800) -10.4029(9.4E-16) -10.4029(2.4E-12) -10.4029(2.5E-16)

    Shekel10 -10.4291(0.7581) -10.5364(8.1E-12) -10.5364(6.1E-14) -10.5364(7.1E-16)

     From Tab.4, SACDE is more precise in 13 30-dimension benchmark function, except Rosenbrock and Noisy Quartic. Besides, SACDE has an advantage in times of evaluation of functions. When it comes to 7 low-dimension benchmark functions, SACDE is obviously worse in

    180 Shekel5 than results reported in literatures. Meanwhile, SACDE obtains at least as good as results

     in literatures for other 6 low-dimension functions. As a conclusion, the proposed algorithm is improved a lot when compared with JADE and DE/rand/1 for high-dimension problems. SACDE

     is more steady and robust based on standard deviation. 3 Water distribution networks optimization

185 3.1 Problem introduction Water distribution network is a complex system comprising of many elements, including reservoirs, pipes and nodes. A proper water distribution network design should meet the requirements of all nodes water consumption and keep the pressure head of every node above given value. Today’s societies require the maximum benefit with minimum cost. In order to

    achieve this goal, some optimization methods should be employed. This paper mainly focuses on 190

    - 8 -


     selecting from commercial available pipes of different diameters for water distribution networks and aiming at minimizing the cost on condition that all water consumption and pressure head

     requirements are satisfied. Due to the complexity of the water distribution network problem, the number of available combinations can be very huge, which make this problem very difficult to [9]195 . solve

     Because of the nonlinear relationship between flow and head pressure, the difficulty of ear programming and designing a reliable water distribution network is further increased. Lin nonlinear programming were commonly used previously, however, those methods cannot obtain global optimal effectively. With the development of computer technology, many intelligent

    200 heuristic searching algorithms are widely used to solve water distribution network problems ,for

     instance, evolutionary algorithm (EA). But these algorithms still have some flaws, for example, excessive computational complexity and violation of constraints of water distribution networks.

     In this chapter, SACDE is used to solve water distribution network problems. EPANET is employed to build the model of water distribution network, because EPANET can handle static

    and dynamic data loading and have interfaces with MATLAB, which make it easy for SACDE to 205

     solve the problem in MATLAB. 3.2 Problem formulation [10] The least cost design of a water distribution network can be stated as follows : Minimize: Cost of the water network design

    Subject to ; 210

     1. Continuity equation 2. Conservation of energy equation

     3. Minimize pressure requirements 4. Other Constraints (maximum pressure; flow velocity,; reliability)

    The objective function (cost of the water network design) is mathematically assumed to be a 215

    cost function of pipe diameters and length,


    C f (D, L) i i i 1

    f (D, L) Where Dand length Lis the cost of pipe i with diameter , and N is the i i ii

     number of pipes in the network. The function is to be minimized under the following constraints.

    220 3.2.1 Continuity equation

    For each node, a continuity constraint should be satisfied,

    Q Q Q in out e

    where is the flow rate into the node, QQQis the flow rate out of the node, and is inoute

     the external inflow rate or demand at the node.

    225 3.2.2 Conservation of energy equation

    For each loop in the network, the energy constraint can be written as follows:

    h E 0 f P

    h Eis the energy added to the water by a pump. where is the head loss and fP

    In this study, the h can be computed by Hazen-Williams Equation: f

    - 9 -


     L ;230 h ~;Qf CD;

    ;is the Hazen-Williams roughness coefficient, Q is the flow rate , is the pipe Where C D ;

    ;diameter, and is the pipe length.L ;

    ;~ is a numerical conversion constant; is a coefficient equal to 1.85; and is a ;coefficient equal to 4.87; The higher the constant ~ , the great the head loss, because higher ~;;

     235value requires larger diameters , thus more expensive water network design.

    3.2.3 Minimum pressure requirement

    For each node in the network, the minimum pressure requirement is given in the following


    min H ? H ; j 1,......, M j j min 240 is the minimum required pressure head Where H is the pressure head at node j , H jj

    at node j , and is the number of nodes in the network. M

     3.2.4 Handling constraints The penalty function is often introduced to effectively guide solution vector from an infeasible solution area that only slightly violates the constraints to a feasible solution area. The

    penalty cost for an infeasible solution is calculated by the distance away from feasible solution 245

    area. In the water network design problem, a penalty function is introduced to prevent the

    algorithm from searching the infeasible solution area, where pipes with small diameters that do not

    satisfy the minimum pressure head requirement at each node are located. The penalty function is

    [11] in the form : min 2 min ) a{max(0, )}250H H b sgn{max(0, H H )}f (Hpenalty j j j j j 3.3 Application Mechanism SACDE will be used to optimize several benchmark water distribution network models. A [12] brief description of the steps of the SACDE for network optimization is given below and illustrated in Fig.1.

    1. Generation of population of points by SACDE. Each point represents a combination of 255

     pipe diameters of the pipe network. 2. Computation of the network cost for each of the solutions after converting the randomly generated pipe sizes to market sizes. 3. Read the network data from the EPANET input file.

    4. Adjust the pipe diameter in the input file. 260

     5. Perform hydraulic analysis of each network with EPANET network solver. 6. Read the nodal pressure from the output file of EPANET. Check the pressure at nodes

     required to meet certain nodal pressures. Compute penalty cost if the nodal head at any node is less than the required minimum.

    7. Calculate the objective function which is the sum of the network cost and the penalty cost. 265

    The value is returned to SACDE for comparison.

    8. Repeat step2~7 until the termination condition is satisfied.

    The flow chart of optimization of water distribution network by SACDE in MATLAB is

    shown in Figure 4-2.

    - 10 -

Report this document

For any questions or suggestions please email