On-line version ISSN 0378-4738
Water SA vol.35 n.4 Pretoria Jul. 2009
Heber Pimentel Gomes*; Saulo de Tarso Marques Bezerra; Paulo Sergio Oliveira de Carvalho; Moisés Menezes Salvino
Laboratory of Power and Hydraulics Efficiency in Water Supply, Federal University of Paraíba, Joao Pessoa, Paraíba, Brazil
This study is aimed at developing a pipe-sizing model for a water distribution system. The optimal solution minimises the system's total cost, which comprises the hydraulic network capital cost, plus the capitalised cost of pumping energy. The developed model, called Lenhsnet, may also be used for economical design when expanding existing hydraulic networks. The methodology developed includes an iterative dynamic calculation process as well as a hydraulic simulation model.
The performance of the method is tested against 4 benchmark examples in the literature. The results obtained show the feasibility of this model, presenting it as a viable alternative for water distribution systems. The method is easily used, once it is performed under EPANET2 software interface.
Keywords: pipe-sizing model, economic optimisation, energy
cost per unit length of the ist pipeline with diameter Di
hydraulic network cost
total distribution system cost
initial cost, related to initial diameter
cost of the diameter which is immediately larger than the diameter related to C1
internal diameter of the pipe
annual rate of increase in the unit cost of energy
unit cost of electrical energy
energy cost gradient
annual interest or discount rate
number of pipeline
number of annual pumping hours
cost gradient value
optimum value of the cost gradient
present worth factor
pipe flow rate
expected period of service for the network
co-efficient of Hazen-Williams equation
pressure gain in the most unfavourable node
pump-motor unit efficiency
Water distribution systems (WDS) design optimisation have been receiving special attention from engineers and researchers of water resources and other related areas, due to the high implementation and operational costs of such systems. Generally, the variables that determine WDS conception and expansion project optimisation are the hydraulic network pipeline diameters and pumping head.
WDS dimensioning is mathematically undetermined, thus allowing for innumerable solutions. Throughout history, several dimensioning methods have been proposed. The first of them, namely classical ones, were restricted to choosing network diameters that provided the hydraulic balance of system. However, the scientific community has been looking to minimise WDS cost for decades.
In the late 1960s, the consolidation of micro-computing in global research centres enabled the development of techniques focusing on water network dimensioning, therefore the first optimisation systems appeared based on mathematical linear programming models (LP) (Karmeli et al., 1968), as well as nonlinear programming (NLP) (Jacoby, 1968) and dynamic programming (DP) (Liang, 1971).
More recently, models based on genetic algorithms (GA) (Dandy et al., 1996), which are known as evolutionary algorithms (EA), have come into use. Some researchers have been using methods based on the organisation and/or evolution of other living species. Eusuff and Lansey (2003) proposed shuffled frog leaping algorithms (SFLA), a meta-heuristic algorithm based on the transformation of frogs and information exchange among the population. Maier et al. (2003) and Zecchin et al. (2006) used a new technique called ant colony optimisation (ACO) to WDS optimisation, based on the analogy of the foraging behaviour of a colony of searching ants, and their ability to determine the shortest route between their nest and a food source. Suribabu and Neelakantan (2006) and Montalvo et al. (2008) applied the particle swarm optimisation (PSO) algorithm. PSO is an EA which utilises swarm intelligence to achieve the goal of optimising a specified objective function. This algorithm uses the cognition of individuals and social behaviour in the optimisation process.
Numerous algorithms have been tested on distribution systems by researchers to obtain the most reliable solutions, using the shortest possible computational time (Biscos et al., 2003). EA methods have presented good results, but they require much more computer time. Cui and Kuczera (2003) highlight the problem of long computation times of some models and propose that such analyses be done using super-computers or by parallel computation. This makes their application by technicians difficult.
Abebe and Solomatine (1998) implemented global optimisation algorithms; 2 algorithms, adaptive cluster covering with local search (ACCOL) and GA, yielded promising solutions enabling a choice between accuracy and required computer time. The proposed optimisation set-up can handle any type of loading condition and neither makes any restriction on the type of hydraulic components in the network nor does it need analytical cost functions for the pipes. Liong and Atiquzzaman (2004) applied a powerful optimisation algorithm, shuffled complex evolution (SCE), in order to find solutions with low processing time. SCE deals with a set of population points and searches in all directions within the feasible space based on objective function. Gomes and Bezerra (2007) and Gomes et al. (2008) proposed an iterative method with a relatively short processing time for the optimisation of the total costs for the expansion and rehabilitation of WDS.
However, in spite of considerable developments as detailed in literature, these techniques have not been accepted in practice (Savic, 2002). Dimensioning of new networks and trial-and-error analyses for extensions are frequent. Jimenez et al. (2007) states that optimisation of water networks is not a rule in engineering yet and programs with a user-friendly interface have only just begun to appear.
This study aims at presenting an optimisation model called Lenhsnet, which is designed to obtain an optimal solution of WDS and provide a friendly interface for engineers. This model, which is connected to EPANET 2.00.12 (Rossman, 2008), provides network pipeline diameters and pumping head as a response to dimensioning so as to find the total minimum cost of the system (implementation cost plus energy cost).
The methodology includes an iteration process, based on an initial solution, in which the distribution network is designed, according to minimum accepted diameters in the project. Such an initial solution has the minimum cost of network implementation, once it is made up by the minimum diameters. However, such a solution is usually not a feasible option, since it provides excessive head losses in network pipelines, resulting in high pumping heads.
Based on the initial solution, the calculation process develops iteratively, in a way in which each consequent solution depends on the previous one. The following solutions will be obtained by increasing, in each iteration, the diameter of one of the pipelines, in a way which will keep the additional cost as low as possible. The iterative process finishes when the configurations of network diameters comply with the restrictions imposed by the project (maximum velocity in pipes and the minimum pressure at nodes).
The algorithm method is associated with a hydraulic simulator, which will provide, at each iteration, the hydraulic balance of the system and the values of the variables of the outflow status of the network water flow (flows, velocities, head losses and pressures). EPANET2 was chosen as simulator as it is widely accepted as the world standard in hydraulic and water quality modelling of WDS. Most of the existing WDS have been modelled within EPANET2 (Biscos et al., 2003).
Once the initial solution is established, the simulation of network outflow is done in order to obtain the pressure in all nodes. Once the most unfavourable node is detected, the optimisation process begins. In each iteration, several diameter configurations will be tested. The effective diameter change, in one iteration, will be defined according to the lowest additional network cost in relation to pressure relief given to the network. The pipeline whose change is confirmed will be the one which provides the lowest cost gradient. The cost gradient (Pg) related to a particular pipeline is given through Eq. (1). It represents the marginal cost of the additional pressure of the most unfavourable node, brought by the change of diameter of the network pipeline by its superior adjacent.
Pg is the cost gradient value ($/m)
C1 is initial cost, related to current diameter ($)
C2 is the cost of the pipe diameter which is immediately larger than the current pipe diameter ($)
Δp is the pressure gain in the most unfavourable node (m)
In each iteration there will be 'i' cost gradients, corresponding to 'i' configurations of diameter changes in the 'i' network pipelines; the optimum value of the cost gradient (Pg*) will be the lowest one among all those calculated. The pipeline corresponding to Pg* will be called potential section (s*). s* will have a new configuration, in which the diameter will be the one which is immediately higher (tested). This last configuration will be the start configuration for the following iteration. The iteration optimisation process follows the aforementioned methodology until the optimal solution objective is obtained.
Once the dimensioning solution is obtained, the next step is to check whether the velocities are within the acceptable maximum limit. If the velocity is greater than the maximum permitted, a new diameter is calculated based on the flow of the pipeline. After the diameter is defined, in case it is available, the initial solution for the model will be changed, and such part will be configured with the new diameter and the other ones with minimum diameters. Dimensioning will finish when the iterative process solution does not present any pipeline with a velocity greater than the maximum established one.
This method deals with 2 dimensioning options. In the first one, the network is supplied by a fixed piezometric level in the reservoir. In this case, the system total cost will correspond to the network implementation cost. In the second alternative, the water is directly propelled to the system or to an elevated reservoir, through pumping and the total cost of the system is calculated from the total cost of the pipeline network plus the capitalised energy cost of the pump station. In this last case, the level in the pumping head will be an extra decision variable in the optimisation process.
Fixed piezometric level
In the first alternative, the condition to stop iterations will occur when the pressure in all of the network nodes reaches, at least, the minimum acceptable value. The total cost of the distribution system (objective function), which includes the network implementation cost, may be expressed by Eq. (2):
C(D) is the cost per unit length of the ist pipeline with diameter Di
Li is the length of ist pipeline
C(Di) is the unitary price of the pipeline 'i' under diameter Di
D is the internal diameter of the pipe
n is the number of the pipeline
Variable pumping head
In the system dimensioning, in which the network is pressurised through a pumping station, this method takes the pumping energy cost during the project service life into account. The total cost of the distribution system (objective function) is the network implementation cost plus energy cost, Eq. (3):
C(D,Q,H) is the total cost of the distribution system
Q is the pipeline flow rate (m3/s)
H is the pumping head (m)cx
ηis the efficiency of the pump-motor unit
Np is the number of the hours of pumping per annum
Ec is the unit cost of electrical energy (unit cost/kWh)
PWF is the present worth factor, Eq. (4)
The PWF for the project extent (t years), which does the conversion of several annual costs to a present value is provided by Eq. (4), for annual interest or discount rate j and annual rate of increase in the unit cost of energy e:
The dimensioning optimisation system takes into account the energy cost through the unit called energy cost gradient (Eg). The Eg represents the updated cost of water pressurisation per elevation meter and is provided by Eq. (5):
Similarly to the previous procedure, the iterative process is executed, following the methodology which has been previously described. At the end of each process iteration the Pg* is compared to the calculated Eg. In case the Pg* is lower than the Eg, the investment cost to reduce energy losses in network pipelines - and consequently to increase the pressure in the most unfavourable node - will be lower than the energy cost to increase its load in the network. The iterative process will continue increasing the diameters of the portions until the Pg* value exceeds the Eg value. Once the latter is obtained, the pumping head is determined in a way which the minimum pressure of the system is equal to the required minimum pressure.
This methodology is synthesised in the flow chart presented in Fig. 1.
Examples of application
This method was applied in the optimised dimensioning of the Hanoi Network problem (Fujiwara and Khang, 1990), New York Tunnels Problem (Schaake and Lai, 1969), Bessa Network (Gomes and Formiga, 2001) and the R-9 Network (Leal, 1995).
Example 1: The Hanoi Network problem
The water distribution system in Hanoi (Vietnam) comprises 3 rings, 34 sections, 31 nodes and 1 fixed-level reservoir (Fig. 2, next page). This network was originally investigated by Fujiwara and Khang (1990), and was later used by several authors (Cunha and Sousa, 1999; Eusuff and Lansey, 2003; Liong and Atiquzzaman, 2004; Suribabu and Neelakantan, 2006; Zecchin et al., 2006; Van Dijk, 2008). The system data are presented in Table 1 and the cost data for pipes in Table 2. The piezometric level is 100 m and the minimum network node pressure is 30 m.
To design the system by Lenhsnet 2 contour systems were considered, concerning the maximum velocity acceptable in network pipelines. In the first, only the diameters which were available through the original reference were considered (Fujiwara and Khang, 1990), which provided high velocity, which were higher than the values which are acceptable in practice. In the second situation, new commercial diameters were added to the original series (1 231.2 to 1 435.4 mm), which increased the possibility search space from 634 to 834, and which reduced the maximum velocity to acceptable limits in practice (v < 3.5 m/s). The unitary costs values of new diameters added in Table 2 were determined through the linear tendency curve of the series of original diameters.
The equation used for load loss was Hazen-Williams formula, Eq. (6), using α value equal to 10.67 (default value EPANET2) and the Hazen-Williams coefficient (C) equal to 130.
The execution of this example through the proposed model brought the results of dimensioning of pipeline diameters and of optimised cost, which are presented in Table 3. The last 2 columns of this table contain the diameters obtained for the dimensioning done with and without maximum velocities restriction of network sections. The table also presents the results of the dimensioning done through other models mentioned in the literature.
Table 4 (next page) shows the pressures and costs obtained with Lenhsnet and other solutions mentioned in the literature.
By using only the diameters originally used by Fujiwara and Khang (1990), the solution presented by Cunha and Sousa (1999) achieved the lowest cost ($6.06 m.), see Table 4. However, pressures originated from these references were lower than 30 m in some network nodes, for α equal to 10.67 (default value EPANET2). The dimensioning solution by Lenhsnet provided an optimal system cost of $6.42 m.
Considering the contour condition, in which a maximum velocity in pipeline is admitted; dimensioning by Lenhsnet provided an optimised cost of $5.50 m. The processing time for this example was less than 8 s, through the use of a 1.60 GHz Intel® CoreTM Duo processor with 2 GB of RAM.
Example 2: The New York tunnels problem
The 2nd example is the New York tunnels problem (see Fig. 3), which is gravity-fed from a single reservoir and comprises 20 nodes connected via 21 pipes. The reservoir is at an elevation of 91.44 m and all the nodes are at zero elevation. The objective of the New York tunnels problem was to determine the most economically effective design for addition to the existing system of tunnels that constituted the primary water distribution system of the city of New York. Pipe diameters are considered as design variables. There are 15 available discrete diameters and one extra possible decision which is 'do nothing' option.
All 21 pipelines are considered suitable for duplication. A full enumeration of all possibilities would require: 1621 = 1.9342 x 1025 evaluations. The system has been used as a benchmark network since 1969 (Schaake and Lai, 1969) to compare various optimisation procedures, particularly by Dandy et al. (1996), Savic and Walters (1997), Maier et al. (2003), Matías (2003), Zechin et al. (2006), Gomes et al. (2008) and Montalvo et al. (2008).
The results thus obtained show that the proposed solution meets with all the hydraulic conditions required for the operation of the system. The present study identifies 6 pipes to be doubled in size (see Table 5). A direct comparison of the optimum solution obtained with the developed procedure and that obtained by other researchers is shown in Table 6. The current optimal solution for this is $38.64 m. and no pressure deficit although this can vary slightly depending on the modelling software and parameters used. The cost obtained with the proposed method was $41.24 m. The processing time for this example was less than 3 s, through the use of a 1.60 GHz Intel® CoreTM Duo Processor with 2 GB of RAM.
Example 3: The Bessa Network
This project consists of projecting the Bessa Network, which was originally presented by Gomes and Formiga (2001). This example was adopted because it used electric energy cost, in addition to pipeline investment costs.
It is aimed at dimensioning, as economically as possible, the network sections and the height of the elevated reservoir, considering the pipeline investment costs and capitalised energy cost of the pump, which pumps water to the reservoir. The pumping head is 30 m. Table 7 contains network data. The minimum pressure required at each node is 25 m and the maximum acceptable velocity is 3.0 m/s. Figure 4 shows the numbers of nodes and network sections.
The equation used for energy loss was Hazen-Williams formula, Eq. (6), which assumes α value as 10.67 (default value EPANET2). Table 8 provides the available diameters, materials, Hazen-Williams coefficients C and unitary prices for pipe investment.
The costs and pumping conditions are: number of pumping hours per annum, Np = 7 300, efficiency of the motor-pump unit, η = 0.75, expected period of service for the network, t = 20 years, the average pumping discharge (Q), equal to 420.43 ℓ/s, power tariff, Ec = R$ 0.20/kWh, annual rate of increase in the tariff of electrical power, e = 6%, and annual discount rate, j = 12%.
Based on the supplied data, the energy gradient value (Eg) is R$89 323.97 (R$ - Brazilian Real (R$ 2.00 = USD 1.00 and R$ 0.24 = ZAR 1.00)). The iterative calculation is processed by increasing the sections diameters and reducing the piezometric level, until the solution which provides the lowest system cost is obtained. The iterative process is ended after the 37th iteration, when the optimal cost gradient (Pg*) reaches R$108 879.58, which is higher than the energy gradient. The dimensioning variables values are obtained in the last iteration.
Once the stop condition is established, this method calculates the final pumping head, which was 15.79 m. When the pumping head is multiplied by the Eg, the current cost of electrical power of the system is determined at R$1 410 488.40. Pipeline network capital cost is R$3 260 811.50. Consequently, the water distribution system total cost is R$4 671 299.90.
Figure 5 shows cost evolution during the iterative process, providing several feasible alternatives. In case the managing agency of the system decides that the initial network cost should be lower than R$2.7 m., the solution to be chosen should be that indicated by the 30th iteration, in which pipeline cost is calculated to be R$2.6 m. The dimensioning results are presented in Tables 9 and 10.
When the results obtained by Lenhsnet are compared with those obtained by Gomes and Formiga (2001), by using the nonlinear programming, it is concluded that the proposed model obtained a more favourable result in economic terms, once it reduces the total network cost value by 14.64%.
Example 4: The R-9 Network
The 4th example (see Fig. 6) is the R-9 Network, a medium-size municipal network in Joao Pessoa city - Brazil, implemented by an urban water company in 1982, which consists of 53 nodes connected via 72 pipes and a reservoir. The input data for this problem are given in FORMIGA (2005) and are presented in Tables 10 and 11 (opposite page). This network has been used previously by some researchers (LEAL, 1995; LOPES, 2003; FORMIGA, 2005). The network data, in EPANET2 solver format, are available from http://www.lenhs.ct.ufpb.br/html/benchmarks.html. Table 12 provides the available diameters, Darcy-Weisbach coefficients and unitary prices for pipe investment.
The optimum solution of $199.39 m. cost units is obtained if the pipes as listed in Table 13 are used resulting in a minimum pressure in the system at node 6 being 17.10 m. When the results obtained by Lenhsnet are compared with those obtained by Formiga (2005) ($202.80 m.), by using non-dominated sorting genetic algorithm - NSGA-II, it is concluded that the proposed model reduces the total cost value by 1.7%.
The method presented in this study, based on a dynamic programming process, enables engineers to follow several dimensioning alternatives through the developed program. In this study, the optimisation algorithm was applied to 4 water distribution systems, whose results were compared to the ones obtained through other models, which used different optimisation techniques (NLP, ACO, GA, SFLA, PSO, etc.). This study showed that Lenhsnet presented total optimised costs, which are similar to the results obtained in the literature.
For various reasons, it is important to point out that it is not reasonable to expect that a WDS project be solved in a totally automated way. Optimisation should be seen as a decision support tool. It is within such a paradigm that Lenhsnet is presented as an attractive and practical alternative for WDS design. The necessary computer processing time is very low, and the interface program with the user is easy, once the method runs in the EPANET2 program.
The authors would like to thank the Brazilian Government for the financial support granted through Centrais Elétricas Brasileiras S.A. (ELETROBRÁS), through Financiadora de Estudos e Projetos (FINEP) and through Conselho Nacional de Desenvolvimento Científico (CNPq).
ABEBE AJ and SOLOMATINE DP (1998) Application of global optimization to the design of pipe networks. Proc. 3rd Int. Conf. on Hydroinformatics. 24-26 August. Copenhagen, Denmark. 989-996. [ Links ]
BISCOS C, MULHOLLAND M, LE LANN MV, BUCKLEY CA and BROUCKAERT CJ (2003) Optimal operation of water distribution networks by predictive control using MINLP. Water SA 29 (4) 393-404. [ Links ]
CUI L and KUCZERA G (2003) Optimizing urban water supply headworks using probabilistic search methods. J. Water Resour. Plng. Mgmt. ASCE 129 5 380-387. [ Links ]
CUNHA MDC and SOUSA J (1999) Water distribution network design optimization: Simulated annealing approach. J. Water Resour. Plng. Mgmt. ASCE 125 (4) 215-221. [ Links ]
DANDY GC, SIMPSON AR and MURPHY LJ (1996) An improved genetic algorithm for pipe network optimization. Water Resour. Res. 32 (2) 449-458. [ Links ]
EUSUFF MM and LANSEY KE (2003) Optimization of water distribution network design using the shuffled frog leaping algorithm. J. Water Resour. Plng. Mgmt. ASCE 129 (3) 210-225. [ Links ]
FORMIGA KTM (2005) Otimização multiobjetivo de projetos de redes de distribuição de água. D.Sc. Thesis. University of São Paulo, São Paulo, Brazil. [ Links ]
FUJIWARA O and KHANG DB (1990) A two-phase decomposition method for optimal design of looped water distribution networks. Water Resour. Res. 26 (4) 539-549. [ Links ]
GOMES H and BEZERRA STM (2007) Reabilitação de Sistemas de Distribuição de Água. In: Gomes H, Garcia R and Rey LI (eds.) Abastecimento de Água O Estado da Arte e Técnicas Avançadas. Editora Universitária da UFPB, João Pessoa, Brazil. 221-240. [ Links ]
GOMES H and FORMIGA KTM (2001) PNL2000: Método prático de dimensionamento econômico de água. Revista Brasileira de Recursos Hídricos. Porto Alegre: ABRH 6 (4) 91-108. [ Links ]
GOMES HP, BEZERRA STM and SRINIVASAN VS (2008) An iterative optimisation procedure for the rehabilitation of water supply pipe networks. Water SA 34 (2) 225-236. [ Links ]
JACOBY SLS (1968) Design of optimal hydraulic network. J. Hydraul. Eng. ASCE 94 (KY3) 641-661. [ Links ]
JIMÉNEZ MR, RODRÍGUEZ K, FUENTES AO and DE LUNA F (2007) Diseño óptimo de redes utilizando um algoritmo genético. In: Gomes H, Garcia R and Rey LI (eds) Abastecimento de Água - O Estado da Arte e Técnicas Avançadas. Editora Universitária da UFPB, João Pessoa, Brazil. 241-260. [ Links ]
KARMELI D, GADISH Y and MEYERS S (1968) Design of optimal water distribution networks. J. Pipeline Division. ASCE 94 (1) 1-10. [ Links ]
LEAL AF (1995) Estudo comparativo de métodos de otimização de redes malhadas pressurizadas. M.Sc. Dissertation. Federal University of Campina Grande, Campina Grande, Brazil. [ Links ]
LIANG T (1971) Design of conduit system by dynamic programming. J. Hydraul. Division. ASCE 97 (HY3) 383-393. [ Links ]
LIONG SY and ATIQUZZAMAN M (2004) Optimal design of water distribution network using shuffled complex evolution. J. Inst. Eng. (Singapore) 44 93-107. [ Links ]
LOPES AV (2003) Otimização do dimensionamento e análise de confiabilidade de redes de distribuição de água. M.Sc. Dissertation. University of Brasília, Brasília, Brazil. [ Links ]
MAIER HR, SIMPSON AR, ZECCHIN AC, FOONG WK, PHANG KY, SEAH HY and TAN CL (2003) Ant colony optimization for the design of water distribution systems. J. Water Resour. Plng. Mgmt. ASCE 129 (3) 200-209. [ Links ]
MATÍAS AS (2003) Diseño de redes de distribución de agua contemplando la fiabilidad, mediante algoritmos genéticos. D.Sc. Thesis. Universidad Politécnica de Valencia, Valencia, Spain. [ Links ]
MONTALVO I, IZQUIERDOA J, PEREZ R and TUNG MM (2008) Particle swarm optimization applied to the design of water supply systems. Comput. Math. Appl. 53 (3) 769-776. http://www.sciencedirect.com/science/article/B6TYJ-4S1JK2S-3/2/2395b1ec52911440f451e424498e7bde (Accessed on 3 September 2008). [ Links ]
ROSSMAN LA (2008) EPANET 2.00.12. U S Environment Protection Agency, Cincinnati, USA. [ Links ]
SAVIC DA (2002) Single-objective vs. multiobjective optimisation for integrated decision support. In: Rizzoli AE and Jakeman AJ (eds.) Integrated assessment and decision support. Proc. 1st Bienn. Meeting Int. Environ. Modell. Software Soc. 24-27 June, Lugano, Switzerland 1 7-12. [ Links ]
SAVIC DA and WALTERS GA (2002) Genetic algorithms for least-cost design of water distribution networks. J. Water Resour. Plng. Mgmt. ASCE 123 (2) 67-77. [ Links ]
SCHAAKE J and LAI D (1969) Linear programming and dynamic programming applications to water distribution network design, Research Report No. 116, Department of Civil Engineering, Massachusetts Institute of Technology, USA. [ Links ]
SURIBABU CR and NEELAKANTAN TR (2006) Design of water distribution networks using particle swarm optimization. Urban Water 3 (2) 111-120. [ Links ]
VAN DIJK M, VAN VUUREN SJ and VAN ZYL JE (2008) Optimising water distribution systems using a weighted penalty in a genetic algorithm. Water SA 34 (5) 537-548. [ Links ]
ZECCHIN AC, SIMPSON AR, MAIER HR, LEONARD M, ROBERTS AJ and BERRISFORD MJ (2006) Application of two ant colony optimisation algorithms to water distribution system optimisation. Math. Comp. Modell. 44 451-468. [ Links ]
Received 15 October 2008; accepted in revised form 3 April 2009.