Water distribution systems have become immensely complex due to ever increasing water demand and sporadic availability of water at the sources. Generally, water management issues are handled through human intervention, which naturally leads to incompetent trial and error procedures. Moreover, the non-linear system dynamics and sequential pumping operations involved in the water distribution systems make the manual control tricky. Past studies suggest that system efficiencies can be tremendously enhanced by automatic control of such systems. However, the use of control algorithms is not yet in vogue among water engineers. Thus, to demonstrate the efficacy of auto-operated water distribution network, a water supply system dataset available in the literature has been considered as a case study in the present research. This study attempts to evaluate the performance of particle swarm optimization and Ziegler–Nichols tuned dynamic inversion control of pumps in maintaining the target water levels in a series of reservoirs, wherein the controller proficiency has been gauged for different initial conditions and non-constant demand inputs. The results indicate good performance and convergence characteristics of both the controllers and thus they can be used for real time operations.

## INTRODUCTION

The role of a water distribution system is to supply water from rivers, reservoirs and boreholes to the industrial and domestic users through a network of pipes via gravity flow or pumping. The efficient operation and control of such water supply systems is vital for meeting inconsistent water demand. Operation of water supply systems requires a range of decisions such as development, management, operational, etc. (Brdys & Ulanicki 1994). Human controlled operations are a common practice in water supply systems. For example, a pump configuration can be changed quite effectively by an experienced pump station operator on the basis of the observed water level in a reservoir and current water demands. However, continuously controlling the pipe flow with a changing demand pattern is an expensive, time consuming and tedious task. Also, it may not be feasible if a system is very complex, such as in the case of sequential pump operations and a series of storage reservoirs. A computer model is then certainly desirable for simulating the system to first understand its dynamics and then optimize it by performing complicated iterations to aid optimal operations.

Benallouch *et al.* (2014) applied model predictive control in a drinking water supply system for regulating the storage in each water tank. Their results depicted good convergence characteristics of the control. Fiorelli *et al.* (2012) applied predictive control for handling water storage in a four reservoir water supply system of Luxembourg. They compared the results of proportional integral derivate (PID) and global predictive control and concluded that global predictive control leads to better management than PID and mechanical control. Miyaoka & Funabashi (1984) developed a pressure and flow regulation technique for a nonlinear water distribution system. Kumar & Kumar (2009) studied linear (PID) and nonlinear (dynamic inversion (DI)) controls, and found the nonlinear technique better for attaining target flows and levels. Rai *et al.* (2013) developed a constraint tuning approach for optimal pump operations in a water distribution system. Bezerra *et al.* (2012) used a fuzzy approach for controlling pressure in a water supply system. Eker *et al.* (2003) used the H∞ control technique for designing a level controller in a water system. PID controllers are one of the most frequently used controllers in the industry (Campisano & Modica 2002). They are easy to adopt and provide satisfactory performance for a wide range of processes. The design of such controllers requires the specification of three parameters: proportional gain, integral gain and derivative gain. So far, great effort has been made to develop methods to reduce the time spent on optimizing the choice of controller parameters (Kuo 1987). The linear control theory provides many tools for tuning a PID controller (Åström & Hägglund 1995). However, when the system consists of many interconnected non-linear subsystems such as pipes, tuning is a challenging task and can lead to non-optimal solutions. In such cases, an efficient tuning method is required.

Particle swarm optimization (PSO) is the latest evolutionary computation technique and has found applications in many areas as it is very simple to execute because of its exceptional computational efficiency (Abido 2002; Gaing 2004; Lin *et al.* 2008). PSO is considered as a simple concept that is easy to apply because it requires adjusting of very few parameters, and delivers good computing speed and precision when compared with other methods such as machine learning, neural network and genetic computation (Castillo *et al.* 2008). Its algorithm connects the social psychology principles in socio-cognition human agents and evolutionary computations. The algorithm is initiated on the basis of the behavior of organisms such as fish schools and bird flocks. In the beginning, the algorithm arbitrarily generates an initial population, which comprises of a number of candidate solutions. Each candidate solution is called a particle or an individual and such a particle is moved by a velocity updated in accordance with its own experience and other particles’ experiences. PSO is a very versatile technique for improving the global and local investigation capabilities as compared to the other optimization techniques. It possesses useful support and shares information reasonably between particles of the population (Abido 2002). The PSO algorithm technique has been used to solve the system identification for a class of nonlinear rational filters (Lin *et al.* 2008). The algorithm properly predicts the rational filter parameters. The PSO technique has also been employed for multi-machine power system stabilizers to solve its optimal design problem (Abido 2002).

The DI controller synthesis methodology cancels the undesirable dynamics of the system and substitutes it by specified desirable dynamics. This task is achieved by cautious algebraic selection of the feedback function. Hence, the DI technique is also known as feedback linearization. Nonlinear DI controls provide better results than linear PID controller (Kumar & Kumar 2009). Hence, in the present study, a PSO tuned DI based PID controller is developed for a water distribution system to control pumps and maintain the target water levels in the reservoirs. Furthermore, its performance is compared with the well-known standard Ziegler–Nichols (Z–N) tuned DI based PID controller for different initial conditions and non-constant demand inputs.

## WATER DISTRIBUTION SYSTEM COMPONENTS

The major components of a water distribution network are pipes, valves, pumps, reservoirs, etc.

### Pipe

*Q*= flow in the pipeline,

*g*= acceleration due to gravity,

*D*= diameter of pipeline, = length of pipeline, = cross-sectional area of pipe,

*Δh*= head difference between two ends of a pipeline and = total loss in a pipe section which is given by Equation (2): where denotes friction loss and denotes local losses.

*f*= Darcy–Weisbach coefficient.

*HWC*is the Hazen–Williams roughness coefficient.

### Pump

*h*

_{p}developed by the variable-speed pumps is given by Equation (5): where

*,*and are the constants for a particular pump depending on its characteristics,

*N*is the speed of the pump, is the pump flow rate and

*n*the number of pumps running in parallel. Here,

_{p}is*n*= 1.

_{p}### Reservoir

*V*is the volume of a particular reservoir , is the variable head of particular reservoir, and is the area of particular reservoir.

## LINEAR AND NONLINEAR CONTROLLERS

### PID controller

*u*), is generally based on the error (

*e*) between the actual output and the desired output. The control variable (

*u*) is given by Equation (7): The terms on the right hand side of Equation (7) are the expressions for proportional, integral and derivative actions, respectively and is the time of intergration. The constants

*K*and

_{p}, K_{i}*K*are the gains (parameters) of the PID controller.

_{d}*K*affects the rise time, i.e. the time required for the system output to rise beyond 90% of the desired level for the first time.

_{p}*K*reduces the overshoot (i.e. how much the peak level is higher than the steady state) and settling time (i.e. the time taken for the system to converge to its steady state).

_{d}*K*eliminates the steady-state error, i.e. the difference between the steady-state output and the desired output.

_{i}### DI controller

*x*is the state variable

*. F*is a nonlinear state dynamic function and

*G*is a nonlinear control distribution function. It should be noted that

*F*and

*G*are both vector functions whereas

*h*is a scalar function. A scalar function is a function that gives a real number to a set of real variables.

*y*(

*t*) is the system output and

*x*

_{des}(

*t*) is the desired response.

*u*. Since the control variable does not appear in the error dynamics equation, a second order error dynamics is considered as shown in Equation (13). Second order error dynamics can be used to shape the system response (Piazzi & Visioli 2001):

## METHODOLOGY

The applicability of PSO tuned DI based PID control of pumps is demonstrated on a water supply system by Eker *et al.* (2003).

### Description of the water network

*et al.*(2003). The adopted water distribution network consists of three pumps and three reservoirs along the supply line. It can be observed from Figure 1 that the manual control of this water system is a difficult task due to a series of pumping stations and successive intermediate storage reservoirs.

*H*indicate the variable heads and

_{D1}, H_{D2}, H_{D3}*H*indicate the static heads in reservoirs, respectively, H

_{SH1}, H_{SH2}, H_{SH3}*is the static head in pipe 4 and*

_{SH4}*L*

_{pipe1}

*, L*

_{pipe2}

*, L*

_{pipe3}

*, L*

_{pipe4}denote the lengths of pipes. The water flow rates through pipes are denoted as

*Q*and

_{1}, Q_{2}, Q_{3}*Q*

_{out}for pipes 1, 2, 3 and 4, respectively.

*A*is the area of reservoir and

_{R}*D*is the diameter of pipe.

*A*are the pump constants.

_{0}, B_{0}, C_{0}Pipes | Reservoirs | Pumps |
---|---|---|

L_{pipe1} = 669.27 m | H = 113.4 m _{SH1} | A = 0.0001433 _{0} |

L_{pipe2} = 13805.04 m | H = 210.4 m _{SH2} | B = 0.005015 _{0} |

L_{pipe3} = 20094.69 m | H = 283.4 m _{SH3} | C = 3.98 _{0} |

L_{pipe4} = 4689.04 m | H = 279.7 m _{SH4} | n = 1 _{p} |

D = 1.4 m | A = 475 m_{R}^{2} | |

f = 0.0183374 | |

Pipes | Reservoirs | Pumps |
---|---|---|

L_{pipe1} = 669.27 m | H = 113.4 m _{SH1} | A = 0.0001433 _{0} |

L_{pipe2} = 13805.04 m | H = 210.4 m _{SH2} | B = 0.005015 _{0} |

L_{pipe3} = 20094.69 m | H = 283.4 m _{SH3} | C = 3.98 _{0} |

L_{pipe4} = 4689.04 m | H = 279.7 m _{SH4} | n = 1 _{p} |

D = 1.4 m | A = 475 m_{R}^{2} | |

f = 0.0183374 | |

### System dynamics of the water system

*X(t)*] and pump speeds (

*N*and

_{1}, N_{2}*N*) are taken as control variables [

_{3}*U(t)*]. There are seven state variables (reservoir levels and flows) and three control variables (pump speeds). The outputs for the present problem are three reservoir levels. Hence, the state and control vectors are defined in Equations (14) and (15), respectively: The system dynamics can be written using the above-discussed Equations (1), (3), (5) and (6) in the following equations:

### Control synthesis procedure

*Q*

_{out}

***. Since the outflow rate is constant after reaching steady state, the derivative of

*Q*

_{out}is zero in equation, Equation (16), and hence the only unknown,

*H*is calculated using Equation (23). Control gains have been derived through PSO as well as standard Z–N method. Internal dynamics of the system have also been checked. The error,

_{D3}**e*, defined as the difference between the reference and the actual value of the reservoir level, is calculated by Equation (24): After substituting the state variables and their derivatives in Equation (13), the modified form is as follows: In Equation (27), only one control variable,

*N*, is the unknown and is quadratic in nature, and hence the equation is solved for

_{3}*N*. Once

_{3}*N*is calculated, Equation (26) can be solved for

_{3}*N*. Similarly Equation (25) is solved for

_{2}*N*.

_{1}### Tuning methods

To obtain appropriate values for the controller parameters, *K _{p}*,

*K*and

_{i}*K*, tuning methods such as PSO, and standard Z–N method (Ziegler & Nichols 1942) are used. These are discussed in the following section.

_{d}#### PSO tuning

The PSO tuning method has been identified as a suitable technique to calculate appropriate values of controlling parameters such as *K _{p}*,

*K*and

_{i}*K*. PSO starts with a random population and searches for optima by updating the population. Every particle is assigned a randomized velocity and is then flown through the problem space. Then at the best-fitness time, the data are shared between the correlated particles. The best positions of the particles are determined by the associated particles over the course of multiple generations and the velocity of every particle is updated accordingly through an iterative process. Kennedy & Eberhart (1997) explained that the particles move in a three-dimensional medium in large numbers and every single particle keeps track of its position and velocity vector along with time until it reaches the maximum level of fitness.

_{d}A potential solution is represented by every particle and its position is represented by position vector, *X _{i}*(

*t*). The particles go towards the best solution initiated by the neighboring particles and also by the particle itself when the particles are flown into the searching space.

Assume that D-dimensional is the searching space and *n* is the number of particles forming the colony. The *i*th particle represents a D-dimensional vector, *X _{i}*(

*t*) for

*i*= 1, 2, … ,

*n*. In searching the space, the position of the

*i*th particle is located at

*X*(

_{i}*t*)

*=*[

*x*(

_{i1}*t*)

*, x*(

_{i2}*t*),

*…*,

*x*(

_{iD}*t*)]. The position is the potential result of every particle. By inputting the position of a particle into the designated objective function, the particle's fitness is calculated. The

*i*th particle's flying velocity is also a D-dimensional vector, represented by

*V*(

_{i}*t*)

*=*[

*v*(

_{i1}*t*),

*v*(

_{i2}*t*), … ,

*v*(

_{iD}*t*)], where

*i*= 1, 2, … ,

*n*.

*P*(

_{i}*t*)

*=*[

*p*(

_{i1}*t*)

*, p*(

_{i2}*t*), … ,

*p*(

_{iD}*t*)] represents the best positions of every particle. Among all the particles, the best particle found up to time

*t*is represented by

*P*(

_{g}*t*)

*=*[

*p*(

_{g1}*t*)

*, p*(

_{g2}*t*), … ,

*p*(

_{gD}*t*)].

*v*is called the velocity of particle

_{iD}*i*;

*x*represents the position of particle

_{iD}*i*with objective value fitness;

*p*is the best historical position of particle

_{iD}*i*itself;

*p*is the global best position; rand

_{gD}_{1}(), rand

_{2}() are uniformly distributed random numbers between interval 0 and 1.

*c*and

_{1}*c*are learning rates and have been considered to be equal to 2.

_{2}*w*denotes the inertia weight and its value is dependent on

*w*

_{max}

*, w*

_{min}and number of iterations. For this study,

*w*

_{max}has been taken as 0.9 and

*w*

_{min}as 0.4. Five hundred iterations were carried out in the present study. The inertia weight is a user-defined parameter that controls, along with

*c*and

_{1}*c*, the previous values of particle velocities on the current one.

_{2}On the basis of the particle's best solution, *P _{i}*(

*t*) and particle's neighbors best solution,

*P*(

_{g}*t*), as calculated by means of Equations (28) and (29), the velocity vector of every particle is adjusted. The objective function guides the PSO algorithm to converge towards the global optimal solution. Hence, before executing the PSO algorithm, the objective function should be defined properly.

In the present study, the integral of the square error (ISE) is used to define the objective function as:

The following steps summarize the procedure adopted in PSO:

Step 1: Initially, the population is generated randomly at certain interval that contains ‘

*n*’ particles (population size).Step 2: The algorithm stops when the prescribed number of generations is achieved.

Step 3:

*p*_{best}and*g*_{best}are recorded by evaluating the corresponding objective function, ISE of every particle.Step 4: Velocity updating of Equation (28) and the position updating of Equation (29) is performed for each particle.

Step 5: Each particle's position is checked using the following formula and the search interval is between 0 and 1:

Step 6: Again, it goes back to Step 2.

#### Z–N tuning

The Z–N frequency response tuning technique is a heuristic method of tuning a PID controller. It is done by setting the integral and derivative gains to zero value. The proportional gain is then augmented until it reaches the ultimate gain, at which the system output has steady and regular oscillations. Ultimate gain and oscillation period are then used to set the controller parameters.

A description of all the symbols used in this paper can be seen in Appendix I.

## RESULTS AND DISCUSSION

The actual time integration of the state equations has been carried out using fourth order Runge–Kutta method with a fixed time step of 0.5 sec. Simulations are run for 15 minutes duration. The tunable parameters of a controller such as proportional gain (*K _{p}*)

*,*integral gain (

*K*) and derivative gain (

_{i}*K*), which are very sensitive, are obtained from PSO and Z–N tuning methods. Following the Z–N method, the values for the ultimate gain (

_{d}*K*) and period of oscillation (

_{u}*P*) are found to be 0.0015 and 288, respectively. Values of

_{u}*K*,

_{p}*K*and

_{i}*K*obtained from the PSO method are shown in Table 2.

_{d}K _{p1} | K _{p2} | K _{p3} | K _{i1} | K _{i2} | K _{i3} | K _{d1} | K _{d2} | K _{d3} |
---|---|---|---|---|---|---|---|---|

0.0656 | 0.0033 | 0.0032 | 0.0017 | 3.3208e-006 | 1.0000e-007 | 0.6963 | 0.0929 | 0.09 |

K _{p1} | K _{p2} | K _{p3} | K _{i1} | K _{i2} | K _{i3} | K _{d1} | K _{d2} | K _{d3} |
---|---|---|---|---|---|---|---|---|

0.0656 | 0.0033 | 0.0032 | 0.0017 | 3.3208e-006 | 1.0000e-007 | 0.6963 | 0.0929 | 0.09 |

Eker *et al.* (2003) performed nonlinear simulations and took observations on the example system considered in this study. From the data obtained from this study, the mean water flow rate, *Q*_{out} was found to be approximately 2.83 m^{3}/s, and the heads in the reservoirs 1, 2 and 3 were measured as 4.2, 2.15 and 3.2 m, respectively. Hence, to check the developed control model, the system was simulated by considering an initial condition based on the Eker *et al.* (2003) study:

1st initial condition: *X* (*0*) = (2.83, 3.20, 2.83, 2.15, 2.83, 4.20, 2.83)

Further, to validate the model, an arbitrary second initial condition was chosen:

2nd initial condition: *X* (*0*) = (2.20, 3.50, 2.20, 2.70, 2.20, 3.80, 2.20)

*Q*

_{out}) is 2.4 m

^{3}/s and the target reservoir levels (

*H*,

_{D1}*H*,

_{D2}*H*) are 4.0, 2.5 and 3.91 m, respectively. Simulations have been made for 15 minutes duration. Depending on this target outflow rate, the third target reservoir level is calculated. Figures 2 and 3, subplots (a), (b) and (c) show how the reservoir levels (

_{D3}*H*and

_{D1}, H_{D2}*H*) vary with time and how they are reaching targets. The control variables, i.e. pump speeds (

_{D3}*N*,

_{1}*N*and

_{2}*N*), are presented in Figures 2 and 3 subplots (d), (e) and (f).

_{3}*Q*

_{out}) (from 2.83 to 2.7 m

^{3}/s, from 2.7 to 2.8 m

^{3}/s and from 2.8 to 2.6 m

^{3}/s) have been imposed on the supply system, for every 30 min duration. Simulations have been made for 1.5 hours duration. The target level of the third reservoir has been calculated from these outflow changes. The remaining two target reservoir levels have been kept constant intentionally (i.e. first reservoir level at 4.0 m and second reservoir level at 2.5 m). In all these results, a stable process was obtained, showing the ability of the control scheme to reach final set points. As can be expected, the control action is more intense during the first steps. Later, the pump speeds show small changes to reach final targets. Figure 4(a)–(c) show the first and second reservoir levels and outflow (state variables), respectively, and Figure 4(d)–(f) show the pump speeds (control variables) variation with time. It can be said that the designed controllers are working properly for all cases. However, in the case of the Z–N method, oscillation is more and it needs more time than the PSO method to reach the steady state.

The results show (Figures 2–4) that in all the cases, Z–N and PSO tuned DI controllers are capable of suitably rejecting the disturbances caused by changes in the state of the system and bring the system back to the desired condition.

The DI controller is a nonlinear controller that give almost zero steady-state error. As can be seen from Figure 4, control action is intense in the initial stages, but later the pump speeds display minor changes in reaching the final targets. It is observed from Figures 2–4 that as compared to the Z–N DI controller, the PSO DI controller is able to achieve the targets quickly without creating undue transients. It is also seen that the PSO DI controller performs better in all the considered cases. Thus, it can be said that the PSO DI is quite robust and can be used for extended period simulation as can be seen from Figure 4.

The results obtained from the developed PSO tuned DI controller have also been compared with Eker *et al.* (2003). In Eker *et al.* (2003), the maximum water level in Reservoir 1 was to be kept below an overflow limit of 4.3 m and hence a level controller was required for a stable operation. The speed of Pump 1 was changed for controlling the water level in Reservoir 1, whereas the other two pumps were assumed to operate at a nominal speed of 985 rpm. A square wave output flow disturbance of ±0.142 m^{3}/s was applied after Reservoir 1 with a frequency of 0.000005 Hz every 27.7 hours up to 138.5 hours.

*et al.*(2003). It can be seen that the output flow disturbance results in negligible change in the head of Reservoir 1 and the developed controller is also able to successfully maintain the desired 4 m water level in Reservoir 1. Eker

*et al.*(2003) reported substantial changes in the water level of Reservoir 2. They observed that the head reached 3.8 m in the increased cycle of flow disturbance and 0.25 m in the decreased cycle. A similar pattern can be seen in the existing results too. In our case, the head reached 3.678 m in the increased cycle of flow disturbance and 0.486 m in the decreased cycle. Also, Eker

*et al.*(2003) found a 0.5 m water head difference between the peaks of increased and decreased cycle for Reservoir 3, while this difference was 0.528 m in the present study. Thus, it can be seen that the results obtained by using PSO tuned DI control are comparable with that of Eker

*et al.*(2003).

## CONCLUSIONS

A PSO and Z–N tuned DI controller have been synthesized for a typical water supply system consisting of a sequence of pumping stations and intermediate storage reservoirs. The automatic control design procedure proposed to achieve the target levels in the reservoir produced satisfactory results. The controller is working effectively for both the problems: (1) constant outflow rate and constant reservoir levels and (2) non-constant demand input, i.e. outflow and constant reservoir levels. The results indicate that the performance of the controllers is good for any given initial conditions. However, when the PSO tuned DI controller is used in place of the Z–N tuned controller, targets can be achieved quickly without creating undue transients. PSO has very good convergence characteristics too. Hence, it could be a very good controller for real time operations.