The aim of the present paper was to move water through a reservoir network in such a way as to meet consumer demands and level constraints, minimise the cost of electricity, and minimise the loss of chlorine. This was to be achieved by choosing the switching intervals of reservoir inlet pumps and valves, at the same time complying with the allowed minimum interval size of each device. Switching combinations that threatened to exceed constraints were rejected heuristically. Flows were balanced by linear programming (LP). The genetic algorithm gave confidence in the near-optimality of its solutions, through the well-defined Pareto fronts between the competing objectives. The method was applied to a 16-reservoir water distribution system in Durban, South Africa. Comparison with an equivalent ‘dead-band’ control showed a 30% improvement in a weighted objective.
The technique reported here has been applied on a model of the Durban Southern Aqueduct trunk-main system, comprising 16 reservoirs and 85 pipe connections. Such distribution networks have a degree of interconnectivity, allowing alternative sourcing and routing of water. Present control of Durban's system is essentially ‘dead-band’, with pumps or valves independently admitting water to reservoirs when a set low level is reached, and closing once the level has risen to a set upper level. This study was spurred on by a crisis incident which occurred in March 2010, in which a combination of water demands threatened supply to part of the network, despite manual intervention in the central control room. The primary operational requirement is to maintain the level in every reservoir between a lower constraint (emergency level for fire-fighting), and an upper constraint representing the reservoir's maximum capacity. Additionally, a secondary objective aims to minimise chlorine losses, to reduce downstream re-dosage. Finally, pumps in the network should ideally make use of the periods of reduced electricity tariff. These three objectives place competing demands on the operation. A Pareto-type analysis which elucidates the associated compromises is thus proposed.
Superimposed on the three objectives detailed above is the need to prevent switching of the pumps too often. If the start-up power surge occurs too often, pumps overheat and are damaged. Though ‘on’ periods could be short, the approach in this study has been to stipulate the same minimum size for both ‘on’ and ‘off’ periods, usually 6 hours for the larger pumps.
Previous attempts at the operational optimisation of the Durban potable water distribution network include Biscos et al. (2002, 2003) and more recently Purdon et al. (2010). These studies recognised the ‘binary’ nature of the operational optimisation. It is not so much a problem of continuously feeding water through the reservoirs to satisfy demand, but rather of deciding when to do each transfer. Pumps are either on or off, and valves are either open or shut. These on/off periods could be coordinated in some fashion to meet the objectives mentioned above.
The main difficulties experienced with the single-objective approaches adopted by these workers concerned the binary nature of the controls, and the minimum interval sizes. Solution times were long, and often only relaxed solutions were possible. Biscos et al. (2003) employed mixed integer non-linear programming for a constrained optimisation of reservoir level control, chlorine concentration control and pumping electricity cost minimisation. The large computational demand of this solution meant that systems including over 10 reservoirs ran slower than real-time – not a useful basis for predictive control. Aiming to include all 265 reservoirs, the recent work of Purdon et al. (2010) focused on finding a faster solution, and included mixed integer linear programming, linear programming (LP) and non-linear programming (NLP). Addition of terms such as (1 − δ)δ, 0 < δ< 1, in the NLP minimisation objective function gave approximately binary switch values δ, and compliance with the minimum switching intervals. However, there was in any case no guarantee of the quality of a solution, owing to the non-linearity of the system. With this in mind, the investigation moved to the ‘random’ methods of simulated annealing, and ultimately a genetic algorithm (GA).
López-Ibáňez (2009) researched the operational optimisation of two small test networks, the larger of which had seven reservoirs, seven pumps and eight pipe junctions. The objective was to minimise the pumping electricity cost, and to minimise maintenance, this latter objective being based on the number of switchings and the shortest ‘off’ time. The Strength Pareto Evolutionary Algorithm (SPEA-2) provided the Pareto front for these two objectives, but this was out-performed by an Ant Colony method.
Cui et al. (2011) used a WATHNET model on a year basis to determine the optimal control settings, such as reservoir ‘trigger’ levels, for the Canberra water distribution system. The fast Non-dominated Sorting Genetic Algorithm (NSGA-II) was used in a multi-objective search minimising incidences of restrictions, incidences of low level and cost.
Pianosi et al. (2011) tackled a related problem, to do with optimal water releases from a catchment lake in Vietnam. The aim was to use an artificial neural network to implement the optimal control law. So, the task of the NSGA-II calculation was to find the best combination of weights on the six neurons used, as represented by the Pareto fronts of three objectives: flooding, irrigation and hydropower. These weights were randomly varied to create the individuals in the population.
Dogaru & Lavric (2011) found an optimal scheduling of water users in a network where contaminated water could be passed on to other users with less stringent requirements. A multi-objective genetic algorithm was used with objectives of fresh water minimisation and capital/operating cost minimisation.
Kang (2014) aimed at a SCADA-based predictive control system using a single-objective genetic algorithm to determine pump on/off settings up to a receding horizon, so as to minimise the electricity cost. The author remarks on the necessity to ‘re-ground’ the starting state (reservoir levels) on each time-step, and adjust the forward horizon for the consumer demand and electricity tariff profiles. A small system including three pumps and three reservoirs required 22 minutes to solve one controller time-step, for a 48-hour horizon. The use of EPANET for the hydraulic model no doubt contributed to this excessive computation time.
Several applications of multi-objective genetic algorithms have concerned the physical design of water networks – e.g. Cheung et al. (2003) used SPEA to target minimum cost, minimum pressure and maximum hydraulic benefit in a piping rehabilitation exercise. Reehuis (2010) used NSGA-II and SMS-EMOA to target minimum cost, maximum nodal pressure and maximum reliability. Babaei et al. (2012) outlined a method using an ant colony algorithm together with the EPANET model, with the objectives of meeting the hydraulic requirements, preserving chlorine and maximising pump performance, by choice of pump types and configuration. Both Li et al. (2012) and Zheng & Zecchin (2014) used NSGA-II optimisations pitting network cost against hydraulic objectives.
Some researchers have considered optimisation incorporating both design decisions and time sequences of network operation. Abunada et al. (2014) included pipe diameters, reservoir elevations and reservoir sizes in a costing which considered operating reliability over 24-hour sequences. Giustolisi et al. (2014) allowed optional additional pumps in their optimisation, and included pipe, pump and reservoir costs in a multi-objective optimisation with energy costs. Within this context, they found that operation based on dead-band level control was more adaptable and reliable than pump-scheduling focused on energy, with little extra energy used. In a similar exercise, Jin & Wu (2014) optimised first pump heads, then associated pump schedules, in a two-phase optimisation. Kurek & Ostfeld (2013) considered three objectives in a multi-objective analysis using SPEA-2: energy cost, chlorine concentration and reservoir sizes. The 24-hour operation sequences changed pump rates over a continuous range at 1-hour intervals.
Genetic algorithms are very popular in the above multi-dimensional problems involving decision-type variables. Certainly, they are capable of providing ‘good’ solutions with far less computation time than a branch-and-bound method. Additionally, they are able to define the trade-offs between multiple objectives. Cîrciu & Leon (2010) compared the features of several algorithms. The simplest form had a single-weighted combination of the objective values, but of course this gives no information about the competition between objectives. The Vector Evaluated Genetic Algorithm makes a selection of the best individuals in each objective to retain in the population. Cîrciu & Leon (2010) pointed out that this has the disadvantage of effectively weighting the objectives according to the selection ratio. The NSGA-II of Deb et al. (2002), on the other hand, has several advantages, which probably account for its popularity in the above studies: It has a low computational demand, preserves diversity and is elitist – i.e. it has a mechanism for eliminating dominated solutions. The extension ‘II’ to NSGA refers to a method for reducing ‘crowding’ of the non-dominated solutions, to get a better definition of the Pareto front.
The objective of the present study will be to explore the potential of a genetic algorithm in determining efficient coordination of water transfers within the Durban Metro's Southern Aqueduct distribution system. This will be developed within the context of model predictive control, with optimisation attempted from the present state forward to a horizon which recedes on each time-step. Single sample optimisations within such a horizon are reported here, not the full closed-loop application requiring state feedback and re-optimisation to the horizon on each controller step. The interactions between the objectives of level control, chlorine loss minimisation and electricity cost minimisation will be elucidated using the NSGA method to sharpen the Pareto fronts limiting the mutually achievable minimisations.
The minimisation or maximisation of an objective function determined by continuous variables follows a path based on local gradients. But what can be done when the objective function depends on a discrete variable, such as the binary choice of a pump switch in this study? Branch-and-bound techniques offer a means of sounding out the resulting combinatorial tree structure in an efficient manner, but are nevertheless computationally intensive. Genetic algorithms have been developed to deal with such problems, usually yielding ‘good’ solutions, but with no guarantee of optimality. These mimic the adaptation and improvement of the human genetic sequence through selective breeding, survival of the fittest, and random mutations which occasionally reveal new potential. Thus, each individual in a population is represented by a series of discrete properties. As these properties are swapped with alternatives in each category, the ‘fitness’ of each individual is monitored in terms of an objective function, and only the best are allowed to proceed to the next stage.
Stream flow switch sequences
The problem to be addressed here is the determination of an optimal sequence of binary switchings (pump on/off or valve open/shut) of a number of streams in the distribution network. Each such switch determines a characteristic flow (capacity) in that pipe connection, which has been determined from measured ramps in the reservoir level data. The full flow may not be attained across pipe junctions, owing to capacity restrictions or shutting on either side, but this balancing issue is dealt with separately by LP. In the meantime, consider that one individual candidate in the optimisation is completely determined by the switching sequence of all switchable streams in the system (Figure 1). A basic interval is used (1 hour in this study), and only integer multiples of this step are considered between switchings. For all pumps and valves, these sequences constitute the parameters determining each ‘individual’ in the genetic pool, and they determine, in combination and in a non-obvious way, a degree of optimality according to a defined objective, e.g. minimisation of deviations from desired reservoir levels.
The aim of generating an optimal sequence of future pump/valve switchings relates to the model predictive control. The system will receive updates of present reservoir levels, and with knowledge of typical future consumer water usage patterns, it will compute the best sequence of water transfers up until a time horizon. In real-time, this computation would be repeated, say, once per hour, and the horizon would simultaneously move forward 1 hour. At each step, only the first newly calculated valve/pump setting would be implemented. On the next controller time-step, the search effort may be reduced by starting with the previous sequence, shifted forward one interval.
Minimum switching interval
Here, nh is the horizon size, as a number of steps of the basic interval, nmin is the minimum number of steps and ε is a random number drawn from a uniform distribution over the range 0–1. The exponent b creates a bias in the distribution – presently it has been set to 2 in order to give a higher proportion of smaller intervals, somewhat more flexible than large ones.
Node volume balances by LP
There are three possibilities for the linking pipes of the network, as follows. (a) A pipe that simply transmits anything fed to it, up to its maximum capacity flow, and in accordance with upstream and downstream restrictions. (b) A pipe that must supply exactly a specified flow at each hour. These pipes are typically used to simulate consumer demand profiles, which must be complied with. (c) A pipe like (a) that allows a maximum capacity flow, but zero flow when the valve or pump is off. In the considered network, these pipes are at the inlets to reservoirs.
Reservoirs in the network constitute nodes, which are able to disconnect the rigid pipework mass-balances. However, pipe junctions in the network are nodes which effectively have no volume and no volume-variation – so they require an exact mass-balance, as in Figure 3. There are four possible mass-balances according to the case of pump and valve. The flow fraction δ = 0 for pipes where valve/pump is off and δ = 1 for pipes with fixed flow. These balances are solved by LP.
At each basic time-step Δt (index t − 1 to t), the pipe capacity vector ct is set for the M streams as: ct = (c1t,c2t, …. cMt)T. In the cases of fixed flow profiles (δ = 1), these capacities are set to the actual flow required at any time t.
for i = all nodes which are junctions (i.e. not reservoirs)
It is notable that the fractions for fixed flows (δ = 1) and pumps/valves that are off (δ = 0) are fixed, but the remaining fractions have values between 0 and 1. When a fractional flow passes through a valve or pump, this does not mean that the valve/pump is in some intermediate state. It just means that it is ‘on’, but because of constraints elsewhere in the system, the full capacity is not utilised.
Integration for volume
Heuristic deflection from constraints
The random on/off allocation for individual candidate p (at all times t), the LP flow solution for time t of individual p and the volume integration over Δt from index t −1 to t for individual p are related as in the calculation sequence of Figure 5. An advantage of the program is the early reduction of the number of bad individuals in the population. Reservoir levels that are moving towards a constraint, and have already crossed an approach threshold (e.g. 10% of scale), are likely to exceed the constraint, which would be represented by a heavy penalty in the objective function. Rather, the trend is aborted by curtailing the current mode (on or off) of the pump/valve feeding this reservoir. This is only done if the resultant aborted interval complies with the minimum switching interval. The change only occurs after the present integration time, so this stream change is correctly handled by all reservoirs, since the loop is arranged to integrate all nodes in parallel. The generation of the comparable dead-band control is achieved by choosing only the largest intervals in the random generation of the on/off strings. The heuristic deflection then causes the system to operate on dead-band control. For the meantime, however, a random collection of intervals over the allowed range is generated, and elimination by the described heuristic deflection detracts only in a minor way from the randomness required.
In this study, the objectives of interest included not only ‘minimum deviation from desired reservoir level’, but also ‘minimisation of pumping cost’ by use of lower-tariff electricity periods, and ‘minimisation of chlorine loss’, which occurs with first-order kinetics in reservoirs.
All three of these are affected in different ways by the chosen pump and valve switching sequence. Genetic algorithms also provide the means of defining the interdependence of such multiple objectives, in a single execution. For example, consider the level control vs. chlorine loss objectives. In general, minimisation of chlorine loss will require smaller reservoir inventories that will push levels below set-point. So, for any one level control objective value, there will be a limited best possible chlorine loss objective value. This point on the two-dimensional objective value plot is called a non-dominated solution. But there will be points with a worse chlorine loss objective value for the same level control objective value, e.g. where the level control value arises more from maintaining levels in the larger reservoirs than the smaller ones. Such points are called dominated solutions. If minimisation is sought for both objective values, these dominated points will lie in an area where either objective value could still be reduced. Conversely, the non-dominated solutions will lie on a two-dimensional curve for two objectives, three-dimensional surface for three objectives, etc. This is called the Pareto front. The Pareto front is not usually obvious in terms of the parameters defining each individual. However, there are useful techniques allowing gradual movement of dominated points towards the Pareto front, to create a clearer representation of it.
Chlorine balances and integration
Here, k is the first-order decay rate constant for the chlorine, and Ct is a vector of reservoir chlorine inventories at time t. The main difficulty is to interpret the nodal chlorine concentrations Xit as stream concentrations xjt for the next time-step. First, xjt = Xit for stream j leaving reservoir i. Next, the Xit for pipe junctions (Vi = 0) are evaluated by dividing chlorine inflow by volume inflow. Finally, the concentration in a stream j leaving a junction i can be set as xjt = Xit. To allow this to propagate through a series of pipe junctions, several iterations are required for each time-step. In the present study, five iterations were adequate, and a value k = 0.05 h−1 was taken as an average (van der Walt 2002). The feed stream to the network was set at xFEED = 2 ppm chlorine, and the reservoir Xi0s were all initialised at a fixed pattern determined by a typical running of the system to ‘steady-state’.
Objective functions for minimisation
All objective values are scaled up by a factor of 100 for plotting purposes.
In Equation (14), Vit is the volume of water in reservoir i at time index t, so that is the final volume at the horizon index t = nh. VSPi is the current set-point for the volume in reservoir i. Vmaxi and Vmini are the high- and low-volume limits for reservoir i, respectively. Ri = Vmaxi − Vmini is the range of reservoir i. In this study, values were taken as follows: nh = 24 (Δt = 1 h steps), wi =1 (all i), α = 3, β = 100. The large value of β aims to discourage violation of constraints.
A refinement of this criterion can be done by minimising xFEED, yet maintaining desired chlorine concentrations at all of the consumer outflows of the network. This refinement was not considered in this study.
Here, fjt = cjtδjt, the flow-rate of stream j during the interval t. A relative electrical power tariff structure Tt in Figure 6 was obtained for large power consumers in Durban, from Ethekwini Municipality (2014). For this study, the optimisation period of 24 hours was synchronous with the period midnight-to-midnight on the tariff. The function fe effectively averages the total flow over the solution period, with a weighting according to Tt. Then, it is adjusted to get a volume turnover frequency per day, and multiplied by τDAY, the time-length of a day. Thus, fe is the weighted total flow expressed as the number of total volume turnovers per day. An assumption has been made that all transfers require electrical power in proportion to this rate.
The above description sets the scene for a random creation of NPOP individuals in a population, each of which can be ranked in terms of the three different objective functions. Now, the NSGA method is used to identify the trade-offs between these objectives, and ultimately choose a compromise optimal switching sequence.
Figure 7 illustrates the non-dominated sorting algorithm of Deb et al. (2002) for a two-objective minimisation problem. The original first non-dominated front is stripped off to reveal a second non-dominated front, and so on until all of the individuals in the population are accounted for. The order that these points occurred in the fronts then determines their overall fitness ranking in the population, i.e. the ‘Sorted’ list. Algorithmically, this is done as follows. For each individual I in the population, a list is made of all of the individuals (B,E) which it dominates (i.e. for which it is better in all objective values). Additionally, the total number of individuals which dominate I is recorded (1). Individuals for which the latter count is zero then form the first non-dominated front. This front is then stripped away, effectively by reducing the number count of each remaining individual, according to the number of dominating individuals removed. The dominating number count of these remaining individuals is then examined for zeros to form the next front, and so on.
A selection of the better individuals towards the top of the sorted list are retained in the population, and also offered the chance to propagate random variations of themselves. As a result, points on a multi-objective value plot gradually move onto a hypersurface defining the best simultaneously attainable objective values, i.e. the Pareto front. The methods used to alter the population are common in GA methods, and include crossover and mutation (Figure 8). In this study, parents are randomly paired in the top fraction (up to 50%) of the ranked individuals. The two children produced by crossover replace the equivalent number of individuals in the bottom fraction – so the parents are retained. For mutation, a specified fraction of the children (e.g. 20%) is subjected to a random change of a specified fraction (e.g. 30%) of their intervals. The crossover and mutation operations are applied to all streams j in an individual, with the crossover to the children randomly changing direction from stream-to-stream (i.e. Child 1 will get the 1,3,5, … intervals of Parent 2 instead of Parent 1 for some streams). On completion of crossover and mutation, a new population of individuals as in Figure 1 exists, and these are now ready for evaluation of their flows, volumes, concentrations and objective values as before.
The method has been applied to a section of the Durban Southern Aqueduct trunk-main water supply system, comprising 16 reservoirs and 85 pipe connections (Figure 9). A 24-hour period was considered (midnight to midnight) with a time-step of 1 hour. Pipe flows at full capacity (pump/valve ‘on’) were established from volume changes in reservoirs, and consumer demand flows out of the network were defined at 1-hour intervals by averaging measured flow records. Consider that the immediate stream feeding each reservoir is a switchable stream with a pump or valve, e.g. stream 15 feeding the reservoir at node 9 is switchable. All of these streams are restricted to a minimum switching interval of 6 h, except stream 37 (3 h) to the reservoir at node 26, and stream 71 (1 h) to the small reservoir at node 22.
Optimisations were conducted with all three objectives (volume, chlorine and electricity) as in Figure 10, and then with selected pairs of the three objectives as in Figure 11. The highlighted dot is the point with the best value of an objective sum: V + 20C + E in each solution. In Figure 10, 41 separate solutions, each with 50 points, have been superimposed. Thus, 41 ‘best value’ dots are shown. Conversely, each plot in Figure 11 is a single, separate solution of 200 points, optimising only the two-axis variables. The thick non-dominated zone (solid dots) in Figure 10 arises from the fact that one is effectively viewing a projection of the three-dimensional surface of the Pareto front.
Why are the three objectives competitive? The volume–electricity conflict is to do with making the best use of the low-tariff pumping period. The volume–chlorine conflict arises from minimising the chlorine inventory by reducing the amount of water held up in the system – i.e. below set-point. The chlorine–electricity conflict again arises because electricity is being minimised by increasing volume hold-up in the low-tariff period. The usefulness of the multi-objective analysis in highlighting objective conflicts is illustrated in Figure 11. In the volume–chlorine solution, the inclusion of electricity in the V + 20C + E criterion for the ‘best’ point (highlighted dot) doubles the V penalty with little decrease in the C penalty.
To illustrate a typical solution, Figures 12 and 13 show the results for ‘best’ point on the volume–chlorine plot, Figure 11, according to the single V + 20C + E criterion (Run 172). In the right-hand panel of Figure 12 are the known 24-hour consumer draw profiles for each reservoir. The left-hand panel shows the determined ‘optimal’ stream switching sequences. In the central panel, with the resultant reservoir level trends, only the minimum to maximum level range is shown for each reservoir, so excursions outside of this range are constraint violations. The nominal set-point is shown as a straight line between these two constraints. The tendency to drop levels, despite the endpoint weighting, is obvious. This is caused by the simultaneous attempts to reduce electrical power and preserve chlorine.
The associated chlorine concentration trends in the reservoirs are shown in Figure 13. The objective function for chlorine sought to minimise the total decay rate, by minimising the chlorine inventory. This economic strategy would require penalty constraints to be effective. Unfortunately, no constraints were set for the chlorine concentration in this study, with the result that most concentrations declined. A further refinement of this strategy would be to target set-point chlorine concentrations at the reservoirs, with minimisation of the total chlorine fed to the system.
For comparison, the 41 best weighted objective values for the three-objective runs (highlighted dots on Figure 10) were 30% less with a mean value of 490, and a standard deviation of 23 (5% of the mean).
The methods used here have given some confidence in the solutions for a 16-reservoir system with 85 pipe links. The computational requirements for the genetic algorithm solutions are modest (e.g. a 50-point solution for 24 control steps, using 50 NSGA iterations, took 5 minutes on a notebook PC). Indications are that much larger systems could be handled for real-time predictive control. At each predictive controller time-step, only the first-control step to the horizon would actually be. The entire calculation is repeated for each step, with the horizon simultaneously moving forward at the same rate. Strictly, the method relies only on having state feedback on each step. This is fine for the reservoir levels, but in most situations, there are few chlorine concentration measurements, so an observer model for chlorine would need to run simultaneously in real-time.
The obtained regular Pareto fronts give credence to the notion that optimal solutions have indeed been found, both for the three-way and two-way objective contests (Figures 10 and 11). This is quite remarkable considering that the system had around 1070 switching interval permutations (16 streams × 5 slots av. × 8 interval sizes av.). The GA technique is clearly well suited to this problem, which resembles DNA sequencing itself! The NSGA approach efficiently resolved the multi-objective Pareto fronts.
Interesting aspects of this study included: (a) the method of ensuring minimum switching intervals, by random selection from an appropriate range; (b) the balancing of junction flows by a LP solution at each time-step; and (c) an early heuristic rejection of individuals from the population which was based on predicted trajectories towards constraints.
The motivation for this work was the need to control an increasingly complex water distribution system, by better means than the existing dead-band switching and manual intervention. A weighted sum of the three simultaneous objective values, based on volume control, chlorine loss and electricity cost, produced a nominal ‘best point’ on the Pareto front of each objective of a large number of solutions. These point positions were clustered near the ideal point of the Pareto front (Figure 10), and produced an average objective value 30% lower than dead-band control.
The authors are grateful for the support they have received under the France-South Africa Scientific Cooperation Agreement, and from the Durban Metro.