Optimization of pipe networks including pumps by simulated annealing

OPTIMIZATION OF PIPE NETWORKS INCLUDING PUMPS BY SIMULATED ANNEALING

A.L.H.Costa, J.L.de Medeiros and F.L.P.Pessoa

Department of Chemical Engineering, Federal University of Rio de Janeiro, UFRJ,

Escola de Química, Bloco E, Centro de Tecnologia, Cidade Universitária ,

Phone: [phone], CEP 21949-900, Rio de Janeiro - RJ, Brazil

(Received:October 14, 1999 ; Accepted: May 18, 2000 )

Abstract - The objective of this work is to present an application of the simulated annealing method for the optimal design of pipe networks including pumps. Although its importance, the optimization of pumped networks did not receive great attention in the literature. The proposed search scheme explores the discrete space of the decision variables: pipe diameters and pump sizes. The behavior of the pumps is describe through the characteristic curve, generating more realistic solutions. In order to demonstrate the versatility of the approach, three different examples of optimal design problems are solved.

Keywords: pipe networks; simulated annealing; optimization;

INTRODUCTION

Pipe networks are important structures in urban and industrial activities. These networks are composed by a set of elements, e.g. pipes, pumps or compressors, valves, etc., interconnected in order to transport a fluid from supply sites to demand locations. The water distribution in a city is an important example of a service provided by pipe networks. The design of a network is a complex engineering task. Due to a high number of alternatives, the traditional procedure, a trial-and-error based scheme, is usually not efficient. Thus, optimization algorithms can be important tools for this kind of project.

The least cost design (LCD) of distribution networks consists in the selection of the network elements associated to minimum costs such that each consumer is supplied according to hydraulic head constraints. Many works approached the LCD as a pipe sizing problem, describing the pipelines in different ways: (i) pipe diameters as continuous variables (Lansey and Mays, 1989); (ii) pipe diameters as discrete parameters and pipelines can present sections with different diameters (split-pipe), where the optimization variables are the lengths of these sections (see Bhave and Sonak, 1992; and, Sherali and Smith, 1997; for several references about this topic); and (iii) pipe diameters as discrete variables and split-pipes are not allowed (Cheng and Mah, 1978; Dandy et al., 1996; Savic and Walters, 1997; Cunha and Sousa, 1999). In this work, it is addressed the latter alternative where the network design is formulated as a combinatorial problem. Dynamic programming (Cheng and Mah, 1978), genetic algorithms (Dandy et al., 1996; Savic and Walters, 1997) and simulated annealing (Cunha and Sousa, 1999) are optimization techniques that have already been used to solve the LCD. The optimal design of a pressure relief header network is a similar combinatorial problem where the simulated annealing method has also been used (Dolan et al., 1989; Cardoso et al., 1994).

This paper presents a simulated annealing (SA) approach for the LCD of pipe networks including pumps. Despite its importance, optimization of networks with pumps did not receive great research attention (Murphy et al., 1994). The literature on network optimization is concentrated in gravity fed systems, without pumps. The SA search explores the discrete space of pipe diameters and pump sizes toward the design configuration with minimum costs. The versatility of the optimization scheme is illustrated with the resolution of three different examples of LCD.

NETWORK MODEL

Consider a network composed by S elements (Spi pipes: k = 1,...,Spi and Spu pumps: k = Spi + 1,..., S) and N nodes (interconnection points: t = 1,..., N). The volumetric flow rate in an element k is q_k (elements are associated to an arbitrary flow direction: if the flow is along the orientation then q_k > 0, otherwise, q_k < 0). The head at a node t is given by H_t and its elevation is z_t. Heads at the terminal nodes of an element k are also represented by H_k1 and H_k2 (according to the assigned orientation: 1® 2). The external flow rate at a node t is w_t.

There are several possible formulations to describe the network behavior (Mah and Shacham, 1978). The adopted model is composed by equations of mass conservation at network nodes and energy conservation equations along the elements:

(1)

(2)

where is the set of elements that are pointing at node t and is the set of elements that are pointing in the opposite direction. The network equations above are valid for an incompressible fluid, however, they can be adapted to a simplified description for compressible fluid flow (Hansen et al., 1991).

In a pipe k, the function F_k(q_k) is the empirical relation that describes the fluid head loss due to friction. The Darcy-Weisbach equation evaluates the head loss as:

(3)

where f is the friction factor, L is the pipe length, d is the pipe diameter, v is the flow velocity and g is the gravitational acceleration. In water distribution systems, the Hazen-Williams equation is widely employed (water properties are implicitly included). This equation can be represented by:

(4)

where C is the Hazen-Williams coefficient (related to pipe roughness) and w is a numerical conversion constant.

For a pump k, - F_k(q_k) is the head characteristic curve. This curve may be represented by a polynomial equation (e.g. Gostoli and Spadoni, 1985):

(5)

The flow direction in a pump is determined by the orientation of the equipment (conventionally, q_k > 0). Additionally, there is a maximum pump flow rate (). To ensure that these conditions are not violated, the model is complemented by the following equations (Costa et al., 1998):

for

(6)

for

(7)

where a and b are slack variables.

The network model is solved by the Newton-Raphson method using analytical derivatives. Since, as it will be presented later, the SA method generates a new design through a specific alteration of the current design, the simulation of a new alternative can use as initial guess the solution of the former design, thus reducing the computational needs.

OPTIMIZATION PROBLEM

The network must be designed such that capital and operating costs are minimized and all demand nodes are supplied with specified fluid flow rates according to predefined head bounds. The decision variables are the pipe diameters and pump sizes. Each pipe diameter d_k is selected from a set SD_k with ND_k discrete diameters, . The size of a pump k is represented by a vector formed by the coefficients of its characteristic curves: . The available commercial pumps are organized in ordered sets SP_k with NP_k pumps, .

Objective Function

The objective function is composed by pipe capital costs and pump capital and operating costs. These costs are considered together using the concept of net present value (Edgar and Himmelblau, 1988). Capital costs of a pipe k are given by:

(8)

where c(d_k) is the cost per unit length of a pipe with diameter d_k. Capital costs of a pump k are evaluated based on the respective head characteristic curve (Walski et al., 1987):

(9)

where q^rated and H^rated are the discharge and head of the pump at the rated conditions (point of best pump efficiency). Operating costs are calculated by:

(10)

where r is the fluid density, h is the pump efficiency, Nop is the number of operating hours per year and pc is the power cost ($/KWh).

Constraints

Head bounds on demand nodes are given by

(11)

where, in the optimization context, nodal heads are functions of pipe diameters and pump sizes:

(12)

Equation 12 corresponds to the response of the network simulation model (Equations 1-7).

Mathematical formulation

Using the objective function and constraints that have been presented, the LCD can be formulated by

(13)

subject to

(14)

(15)

(16)

where F_NPV is the factor that transforms the operating costs ($/y) to the present worth ($).

PIPE NETWORK OPTIMIZATION

Simulated annealing

Simulated annealing (SA) is a stochastic optimization method based on the physical process of annealing. In this thermal process, a material is heated to an elevated temperature and cooled slowly to achieve a minimum energy state. The computational reproduction of the annealing process originated the SA method (Kirkpatrick et al., 1983). Starting from an initial configuration of the decision variables, a neighbor configuration is selected randomly. If there is a reduction of the objective function, the new configuration becomes the current configuration, otherwise, the new configuration is accepted or not according to a certain probability. This process of movement-acceptation is repeated until a specified stopping criterion is achieved. During the search, the acceptance probability is gradually reduced. Since the search can accept "wrong way" movements, the method can migrate among different local optima toward the global optimum.

Proposed scheme

The SA method was adapted to be applied to the LCD. Each trial configuration is a vector of pipe diameters and pump sizes and corresponds to an alternative of network design. If a new design violates the head bound constraints, it is discarded immediately, otherwise, its costs are evaluated and the acceptance tests are checked. The following subsections present some features of the proposed optimization scheme.

Movements

The generation of a new configuration is done by a pipe or a pump movement. The selection of which movement is applied obeys a stochastic criterion. The probabilities of a pipe or a pump movement are Spi / S and Spu / S respectively.

Pipe movements can be conducted through three options: (1) a pipe k is selected randomly and its diameter is substituted by the immediately lower diameter in SD_k set; (2) the pipe diameter is substituted by the immediately upper diameter in SD_k set; and (3) the simultaneous application of the two former movements to different pipes. The choice of which pipe movement will be executed is random. Computational experiments indicated that the introduction of movement 3 can bring a significant improvement on the method performance.

Pump movements are similar to pipe movements: a pump k is selected randomly and its size is substituted (1) by the immediately lower size or (2) by the immediately upper size in SP_k set. However, since the reduction of a pump size causes a considerable head decrease along the entire network, pump movement 1 frequently leads to a design with head deficit, i.e., an unfeasible design that is immediately abandoned. Thus, the transition of the search toward lower pump sizes becomes difficult. In order to avoid this obstacle, movement 1 is applied, according to a probability of 50 %, simultaneously with the increase of the diameter of a group of pipes (the number of pipes that must be changed is a parameter that depends on the network size). To ensure the feasibility of the new configuration, the diameter enlargement is applied to the pipes that most affect the heads along the network. These pipes are identified using the significance index (Arulraj and Rao, 1995):

(17)

Cooling Schedule

In the SA search, the acceptance probability of a new worse configuration is given by exp( -DOF/ T) (Metropolis et al., 1953) where DOF is the difference between the objective functions of the new and the current configurations and T is the temperature, the parameter that controls the process. As the search advances, the temperature is reduced (simulating a cooling process) and the method becomes more restrictive.

The initial value of the temperature is defined according to an empirical proposal based on an estimate of the fraction of the movements accepted in the beginning of the search. Assuming ND_k = ND, "K :

(18)

where X = 0.95 is the initial expected acceptance probability. The temperature reduction obeys the equation proposed by Aarst and Van Laarhoven (1985):

(19)

where s is the standard deviation of the objective function values of all new configurations analyzed in the former temperature level and d is the cooling parameter. Temperature is reduced when the number of movements in a temperature level reaches a predefined limit or when a movement implies a reduction of the objective function value, according to the concept of nonequilibrium proposed by Cardoso et al. (1994) (in the first execution of the inner loop, the temperature is only reduced when the number of movements reaches the predefined limit). The adopted stopping criterion interrupts the search if no movements are accepted in an entire temperature level.

Algorithm

Let d be the vector of the network pipe diameters and let be the matrix formed by the pump size vectors. With these definitions, the proposed algorithm is presented below:

Step 1: Initialization

Define algorithm control parameters: initial temperature (T), maximum number of movements per temperature level () and the cooling parameter (d). Select the initial configuration, which must be a feasible point: , and . Initialize the best solution found: d^sol = d^old, and OF^sol= OF^old. Initialize the counter: N_mov= 0.

Step 2: Movement

If (random number Î[0,1] ) £ (Spi / S) then apply a pipe movement otherwise apply a pump movement. Generate the new configuration: d^new and . Update the movement counter: N_mov = N_mov + 1.

Step 3: New Configuration Acceptation

Solve the network model for the new configuration. If the constraints are violated then discard the new configuration and go to Step 4, otherwise, evaluate the objective function OF^new = OF(d^new, ) and DOF = OF^new - OF^old.

If DOF £ 0 then (a) the movement is accepted: d^old= d^new, and OF^old= OF^new; (b) update the best solution found: d^sol = d^old, and OF^sol= OF^old and (c) go to Step 5 (except for the first execution of the inner loop).

If DOF > 0 then the acceptance criterion is checked: if (random number Î[0,1] ) £ exp( -DOF / T) then the movement is accepted: d^old= d^new, and OF^old= OF^new; otherwise the movement is rejected.

Step 4: Inner Loop End

If then go to Step 5, otherwise go to Step 2.

Step 5: Temperature Cooling

Update the temperature according to Equation (19) (if s = 0then use the previous value of s) and reinitialize the counter: N_now = 0.

Step 6: Outer Loop End

If in the last temperature level, none of the movements was accepted then stop, the optimal design is d^sol and ; otherwise, go to Step 2.

NUMERICAL RESULTS

In order to illustrate the application of the optimization method, three different problems of water distribution networks were studied - Network 1: without pumps; Network 2: with one pump; and Network 3: with one pump and one extra reservoir. Figure 1 is the general structure that contains all analyzed networks.

All pipes are 2500 m long and the parameters of the Hazen-Williams equation are w = 10.5088, C = 130, a = 1.85 and b = 4.87. The head lower bound is H^min = 30 m for all demand nodes. Demand flow rates and node elevations are presented in Table 1 (1m³/ h = 2.7778m³/ s ).

The available discrete diameters for all pipes and the respective costs per meter are shown in Table 2. The head characteristic curves of the available pumps are represented by a polynomial of degree 2:

- F(q) = aoq2 + a1q + a2 (20)

- F(q) = aoq2 + a1q + a2

(20)

where q is in m³/s and F is in m. The coefficients of the curves are shown in Table 3 where pump size 1 corresponds to a "no pump" option. It is adopted only one curve that relates pump efficiency and flow rate for all pumps:

h(q) = 695.4q2 + 418.3q + 2.857 (21)

h(q) = 695.4q2 + 418.3q + 2.857

(21)

where q is in m³/s and h is in %.

The constant for the pump capital costs (Equation 9) is K^pump = 700743 where q is in m³/s and H is in m. The power cost is 0.12 $/KWh. The interest rate is 12 % and the project life is 20 years. It is considered that the pump operates during the entire year. All economic data were extracted from Walski et al. (1987). The head characteristic curve of pump size 10 and the efficiency curve were also obtained through the interpolation of the data presented in Walski et al. (1987).

The SA control parameters employed are d = 0.01 and number of elements. The lower pump movement is associated to the enlargement of one pipe diameter. The initial configuration adopted is formed by the largest diameters and pump sizes applied to all elements. Due to the stochastic nature of the SA method, each network example was solved in 10 different runs. The algorithm performance is shown through the objective function of the best solution found, the average objective function in the solution set and the average number of simulations executed.

Network 1

The first network analyzed does not have pumps, the driving force of the flow is the reservoir elevation (node 1). Nodes 10 and 11; and, elements 12 and 13 are not present, thus the network is formed by N = 9 nodes and S = 11 pipes. Since there are 10 available discrete diameters for each pipeline, the search space contains 10¹¹ elements. The best solution found through the application of the proposed algorithm presents a cost of $2,610,500. The average cost in the 10 runs is $2,611,500 and the average number of simulations employed is 13454. It is important to note that the cost of the initial configuration is $9,517,750, i.e., the method reached an optimal design with a cost reduction of 72 %. The pipe diameters and the nodal heads of the optimal design are presented in Table 4 and Table 5.

Ten additional runs were conducted using only pipe movements 1 and 2 to confirm the importance of the introduction of pipe movement 3. The comparison of the results is displayed in Table 6. It can be observed that pipe movements 1, 2 and 3 can reach better results with equivalent computational effort when compared to pipe movements 1 and 2.

Network 2:

This network includes one pump (element 13) connected to the unique source node available (node 11). Node 10 and pipe 12 are not present. The network is formed by N = 10 nodes and S = 12 elements. The search space contains 10¹² discrete points and the initial configuration has a cost of $13,861,716. The cost of the best solution found is $5,505,050 (60 % of cost reduction). The average cost in the 10 runs is $5,532,900 and the average number of simulations is 13342. Despite the increase in the problem complexity, due to the inclusion of the pump, the optimum was found with no increase of the computational effort. The optimal pump size is the alternative 4 and the optimal pipe diameters are shown in Table 7. Nodal heads are presented in Table 8.

The progression of the SA search can be verified in Figure 2 where the values of the objective function of the accepted configurations are plotted.

Network 3:

In this network, there are two source nodes (10 and 11), one of them connected to a pump (element 13 is connected to node 11). In this case, the pump operating point is initially unknown, i.e., the optimization cannot handle the pump as a fixed head increase independent of the pump flow, the pump behavior must be represented by the head characteristic curve. The network has N = 11 nodes and S = 13 elements. The best solution found has a cost of $5,040,013, the average cost in the run set is $5,063,199 and the average number of simulations is 11817. The optimal design selected the pump size 2 and its pipe diameters are shown in Table 9. Nodal heads are shown in Table 10.

An interesting extension of the LCD is the inclusion of the study of the pump optimal location. Establishing several pumps along the network and including "no pump" options in the available pump sets, the proposed algorithm could find the better location to the pump.

CONCLUSIONS

In the literature, the problem of optimal network design is frequently formulated as a search for optimal diameters, whereas the analysis of the pumping stations is relegated to a less important position. The present paper describes a SA approach for the solution of the least cost design problem of pipe networks with pumps. The adopted network model presents a realistic representation of the pump behavior, including the head characteristic curve. Starting from an initial configuration of pipe diameters and pump sizes, the method explores the search region through random movements, generating a sequence of designs. New designs with higher costs are accepted or not according to a stochastic criterion. Network simulations check if the hydraulic constraints are violated. In order to illustrate the application of the proposed optimization scheme, three network design problems were solved. Further studies must be done to optimize the values of the SA control parameters and to investigate the application of the proposed algorithm to classical test cases in order to compare its performance with previous works.

ACKNOWLEDGEMENTS

The authors are grateful to the Research Support Foundation of the State of Rio de Janeiro (FAPERJ) and the Brazilian National Council for Scientific and Technological Development (CNPq) for their financial support.

REFERENCES

• Aarst, E. H. L. and Van Laarhoven, P. J. M., Statistical Cooling: a General Approach to Combinatorial Optimization Problems, Philips J. Res., 40, p. 193 (1985).
• Arulraj, G. P. and Rao, H. S., Concept of Significance Index for Maintenance and Design of Pipe Networks, Journal of Hydraulic Engineering, 121, No. 11, p. 833 (1995).
• Bhave, P. R. and Sonak, V. V., A Critical Study of the Linear Programming Gradient Method for Optimal Design of Water Supply Networks, Water Resources Research, 28, No. 6, p. 1577 (1992).
• Cardoso, M. F., Salcedo, R. L. and Azevedo, S. F., Nonequilibrium Simulated Annealing: A Faster Approach to Combinatorial Minimization, Ind. Eng. Chem. Res., 33, p. 1908 (1994).
• Cheng, W. B. and Mah, R. S. H., Interactive Synthesis of Pipeline Networks Using Pigraph, Computers and Chemical Engineering, 2, p. 133 (1978).
• Costa, A. L. H., de Medeiros, J. L. and Pessoa, F. L. P., Steady-State Modeling and Simulation of Pipeline Networks for Compressible Fluids, Brazilian Journal of Chemical Engineering, 15, No. 4, p. 344, (1998).
• Cunha, M. C. and Sousa, J., Water Distribution Design Optimization: Simulated Annealing Approach, Journal of Water Resources Planning and Management, 125, No. 4, p. 215 (1999).
• Dandy, G. C., Simpson, A. R. and Murphy, L. J., An Improved Genetic Algorithm for Pipe Network Optimization, Water Resources Research, 32, No. 2, p. 449 (1996).
• Dolan, W. B., Cummings, P. T. and LeVan, M. D., Process Optimization Via Simulated Annealing: Application to Network Design, AIChE Journal, 35, No. 5, p. 725 (1989).
• Edgar, T. F. and Himmelblau, D. M., Optimization of Chemical Processes, McGraw-Hill, New York (1988).
• Gostoli, C. and Spadoni, G., Linearization of the Head-capacity Curve in the Analysis of Pipe Networks Including Pumps, Computers and Chemical Engineering, 9, No. 2, p. 89 (1985).
• Hansen, C. T., Madsen, K. and Nielsen, H. B., Optimization of Pipe Networks, Mathematical Programming, 52, p. 45 (1991).
• Kirkpatrick, S., Gelatt, C. D. and Vechi, M. P., Optimization by Simulated Annealing, Science, 220, p. 671 (1983).
• Lansey, K. E. and Mays, L. W., Optimization Model for Water Distribution System Design, Journal of Hydraulic Engineering, 115, No. 10, p. 1401 (1989).
• Mah, R. S. H. and Shacham, M., Pipeline Network Design and Synthesis. In: Drew, T. B. et al. (editors), Advances in Chemical Engineering, 10, Academic Press (1978).
• Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A. and Teller, E., Equation of State Calculations by Fast Computing Machines, J. Chem. Phys., 21, p. 1087 (1953).
• Murphy, L. J., Dandy, G. C. and Simpson, A. R., Optimum Design and Operation of Pumped Water Distribution Systems, Proc., Conf. on Hydraulics in Civil Engineering, Institution of Engineers, Brisbane (1994).
• Savic, D. A. and Walters, G. A., Genetic Algorithms for Least-Cost Design of Water Distribution Networks, Journal of Water Resources Planning and Management, 123, No. 2, p. 67 (1997).
• Sherali, H. D. and Smith, E. P., A Global Optimization Approach to a Water Distribution Network Design Problem, Journal of Global Optimization, 11, No. 2, p. 107 (1997).
• Walski, T. M., Brill, E. D., Gessler, J., Goulter, I. C., Jeppson, R. M., Lansey, K., Lee, H. L., Liebman, J. C., Mays, L., Morgan, D. R., Ormsbee, L., Battle of The Network Models: Epilogue, Journal of Water Resources Planning and Management, 113, No. 2, p. 191 (1987).

Publication Dates

• Publication in this collection
  16 Mar 2001
• Date of issue
  Dec 2000

History

• Received
  14 Oct 1999
• Accepted
  18 May 2000