The main operational objectives for water utilities include the reduction of water leaks, and management of drinking-water quality. Leakage reduction is achieved by controlling average zone pressure (AZP) within water distribution networks (WDNs), while satisfying minimum service requirements [UBRR00, WAPS15] . Pressure control schemes are implemented through pressure reducing valves (PRVs), which reduce pressure at their downstream node. The problem of optimal placement and operation of PRVs in WDNs has been formulated in previous literature, and solved using both evolutionary algorithms [ARC06, NZ09] and mathematical optimisation methods [EM12, PAS19] .
Monitoring and control of disinfectant residuals in drinking water distribution networks is critical to maintain the water quality and eliminate the risks of contamination with pathogens such as bacteria and viruses in distribution [PUW + 02, ASGK14, SLA20] . This is particularly critical during the current COVID-19 pandemic as leaking sewage from sewer networks could allow potentially harmful contaminants into drinking water networks [IRFS17, QWM + 20] . In order to deactivate any pathogens that might exist in distribution networks, disinfectant is typically added at water sources (e.g. water treatment plants), with chlorine being a commonly used water disinfectant. Because chlorine is reactive, it is depleted over time as it travels across the pipe networks, causing a reduction in the ability to prevent microbial contamination. Water utilities aim to maintain a target chlorine concentration, which is sufficient to safeguard public health, while avoiding excessive chlorination, resulting in taste and odor problems, as well as the growth of disinfection bi-product. In addition, the objective is to maintain optimal and constant chlorine concentrations, as variations in chlorine concentration are perceived as water quality problems by customers. Chlorine booster stations are used to deal with this challenge [BTU + 98, PU04, LPP + 07, OO14] . Using booster chlorination, disinfectant is re-applied at selected locations within the network, leading to a more uniform spatial and temporal distribution of chlorine residuals. Previous literature has modeled the operation of booster stations assuming known flow velocities across network pipes -see as examples [BTU + 98, PUW + 02, PU04, LPP + 07]. However, this can lead to sub-optimal design and operation of WDNs. In fact, in order to mitigate disinfectant decay reactions, network operators should aim to reduce travel time from water sources to demand nodes. This can lead to sub-optimal pressure management schemes, where minimum pressure constraints are not satisfied, as observed in [KL10] . Therefore, we consider the joint optimisation of hydraulic pressure and flows, together with chlorine residual concentrations in WDNs.
We investigate the problem of minimising average zone pressure, while penalising deviations from chlorine target concentrations, and satisfying regulatory constraints on pressure and chlorine concentration levels. [Ost05] and [KL10] implemented genetic algorithms to solve problems of optimal operation of WDNs, where optimisation unknowns include network flows and chlorine concentrations, while locations of PRVs and chlorine booster stations are fixed. However, pressure reducing valves and booster stations should be optimally placed for a more effective pressure control and management of chlorine residual concentrations.
In this manuscript, we propose a new mathematical framework for optimal placement and operation of pressure reducing valves and chlorine booster stations. The considered objective is the minimisation of average zone pressure, while penalising deviations from target chlorine concentrations at demand nodes. The transport of chlorine through each pipe is modeled by a one dimensional first-order advection PDE [RB96] , where flow velocity corresponds to the one-dimensional velocity field, and a linear function is used to represent chlorine decay [HWF + 02]. We implement an implicit upwind scheme to discretise the considered PDE. Optimisation constraints include quadratic equations modelling head loss due to pipe friction [EM15, PAS17] , and bilinear terms due to the presence of unknown flow velocities within the discretised advection PDE. In addition, binary variables are used to model the direction of flow across pipes, and the placement of valves and booster stations. The resulting optimisation problem is a non-convex mixed integer nonlinear program.
In comparison to previous literature [Ost05, KL10] , which relied on genetic algorithms, we investigate the application of mathematical optimisation methods to compute good quality feasible solutions for the considered problem, with guaranteed bounds on their level of sub-optimality. We propose a solution method based on polyhedral relaxations of the non-convex terms. A tailored optimisation based bound tightening algorithm is implemented to tighten the relaxations. Furthermore, we develop a heuristic method to compute feasible solutions for the original problem, after solving the relaxed linear problems. In comparison, we show that off-the-shelf global optimisation solvers failed to generate feasible solutions for the considered problem. The performance of the proposed methods is investigated using multiple problem instances for different WDN case studies.
A WDN with n n demand nodes, n 0 source nodes (e.g. water sources, water treatment plants), and n p links is modelled as a directed graph with n n + n 0 vertices and n p edges. Define P := {1, . . . , n p } and N := {1, . . . , n n }. Consider network operation within a time interval [0, T ], where T > 0. We formulate the problem of optimal placement and operation of pressure reducing valves and chlorine booster stations, with the objective of minimising average zone pressure, while penalising deviations from target chlorine concentrations at demand nodes. In particular, we assume flows across links as unknowns within the formulation of the water quality modelling equations. Therefore, it is not possible to use the same formulation for optimal booster placement considered in previous literature [PU04, LPP + 07], where flow velocities are assumed as known.
In this manuscript, we consider the following PDE to govern the chlorine concentration r l (x, t) (measured in mg/l) along link l:
where L l and D l are length and diameter of link l (measured in meters), respectively, and q l (t) is the unknown flow across link l at time t (measured in m 3 /s). Moreover, α l > 0 is the first order decay coefficient associated with pipe l [HWF + 02], measured in s −1 . Equation (1) models the advective transport of a constituent with first order decay. The problem of water quality modelling in WDNs has been extensively investigated in previous literature. As a result, adequate simulation tools have been implemented -see as examples [RB96, MHL14] . Our aim is not to derive a simulation model for water quality analysis, rather to investigate approximate models suitable for implementation within mathematical optimisation problem formulations. Water quality simulations in WDNs require solving equation (1) for each pipe. Lagrangian solution methods have been shown to be accurate and stable, and they are implemented in various modelling software packages [Ros00, MHL14] . These schemes track water parcels as they move through the network, resulting in a set of linear equations whose unknowns are nodal concentrations. In order to write these relations in analytical form, Lagrangian approaches require flow velocities to be known in advanced. In comparison, the discretised set of equations resulting from Eulerian methods can be written considering flow velocities as unknowns.
We implement an Eulerian implicit upwind scheme, where backward differences are used to approximate both temporal and spatial derivatives. The implemented scheme is first-order accurate and unconditionally stable [LG87, ICC99] . First, we consider a time discretisation t k = k∆ t , k ∈ {0, . . . , n t }, with ∆t = T nt . Then, for each link l ∈ P, we introduce a space discretisation x j,l = j∆x l , j ∈ {0, . . . , J l }, with ∆x l = L l J l . We write r j,l (t k ) := r l (x j,l , t k ) for all j ∈ {0, . . . , J l } and k ∈ {0, . . . , n t }. For all j ∈ {1, . . . , J l }, l ∈ P, and k ∈ T = {1, . . . , n t } :
with given initial concentrations c 0 ∈ R nn+n 0 at network nodes, measured in mg/l.
In this section, we describe the problem formulation for the optimal placement and operation of pressure reducing valves and chlorine booster stations. First, we introduce decision variables and constraints related to network hydraulic properties. Let k ∈ T . Source nodes are assumed to have known hydraulic heads h 0 (t k ) ∈ R n 0 , measured in meters. Let d(t k ) ∈ R nn be the vector of known demands, measured in m 3 /s. We denote by q(t k ) ∈ R np the vector of unknown flow rates (measured in m 3 /s), while h(t) ∈ R nn is the vector of unknown hydraulic heads (measured in meters). Moreover, let A 12 ∈ R np×nn be the following link-node incidence matrix:
In addition, let A 10 ∈ R np×n 0 be a link-known head node incidence matrix. We introduce the vector θ(t k ) ∈ R np to model the unknown frictional head losses across network links. Frictional head losses are often represented by either the Hazen-Williams (H-W) or the Darcy-Weisbach (D-W) equations [DLWB15] . Since both formulae involve non-smooth terms, quadratic approximations have been proposed and used in previous literature [EM15, PAS17] . Let a ∈ R np and b ∈ R np be vector of coefficients of these approximations. We model frictional head losses as
Pressure reducing valves reduce pressure at their downstream node, introducing additional head losses, which are represented by η(t k ) ∈ R np . Vector of binary variables v ∈ {0, 1} 2np models the placement of control valves. We have:
v l = 1 a valve is placed on link l in the positive flow direction 0 otherwise,
and v np+l = 1 a valve is placed on link l in the negative flow direction 0 otherwise.
These binary variables are subject to the following physical and economical constraints:
The following constraints formulate energy and mass conservation laws, and the placement of pressure reducing valves on network links:
Where, for all l ∈ P and k ∈ T , we set
Next, we describe variables and constraints associated with water quality. Let c(t k ) ∈ R nn+n 0 be a vector of nodal chlorine concentrations, for k ∈ T , and c max ∈ R nn+n 0 be the vector of maximum allowed chlorine concentration at network nodes. We also introduce auxiliary variables w j,l (t k ) for all j = 0, . . . , J l and l ∈ P, such that:
and rewrite (2) as:
for j = 1, . . . , J l , l ∈ P, k ∈ T , where γ l = 10 −3 (4∆t)/(πD 2 l ∆x l ). Furthermore, r 0,l (t k ) is assumed to be equal to the concentration of the upstream node, depending on the flow direction:
Given a node i ∈ {1, . . . , n n + n 0 }, let I in i and I out i be the index sets corresponding to links with assigned direction entering and leaving the node, respectively. Our problem formulation considers as free decision variables concentrations at source nodes c i (t k ), i ∈ {n n + 1, . . . , n n + n 0 }. Moreover, let v b ∈ {0, 1} nn be a vector of binary decision variables, modelling the placement of a chlorine booster station at network nodes, i.e. v b i = 1 if a chlorine booster station is placed at node i, v b i = 0, otherwise. The number of chlorine boosters to be installed is set by the constraint
Chlorine concentration at unknown head node i ∈ N is governed by the following mixing equations:
where slack variable y i (t k ) ≥ 0 is introduce to model the additional constituent mass injected by a booster, and y max i (t k ) are sufficiently large positive constants, for all i ∈ N , and k ∈ T . In addition, vectors of auxiliary variables ρ(t k ) ∈ R np are subject to the following linear constraints:
Let j ∈ {0, . . . , J l }, l ∈ P, k ∈ T , and define r max
. Finally, we include the following lower and upper bounds:
In order to write the overall problem formulation, we define the following vectors of hydraulic variables q : (1 + J l ) and r(t k ) = [r 1,1 (t k ), . . . , r J 1 ,1 (t k ), . . . , r Jn p ,np (t k )] T . We set r = [r(t 1 ) T , . . . , r(t nt ) T ] T . Analogously, let w ∈ R ntJ be a vector whose components are w j,l (t k ). Finally, let y := [y(t 1 ) T , . . . , y(t nt ) T ] T .
The objective of this study is to minimise average zone pressure in water distribution networks, while penalising deviation from target chlorine concentrations. Average Zone Pressure (AZP) is defined as the following weighted sum of nodal pressures [WAPS15] :
where ξ ∈ R nn is the vector of known nodal elevations (measured in meters). Weights in the definition of f AZP (·) are defined as follows:
Let c * ∈ R nn be a vector of target chlorine concentration at network nodes (measured in mg/l). Moreover, setd
We define the Average Target Deviation (ATD) as
The absolute value within f ATD (·) can be reformulated as a linear function by introducing auxiliary variables s ∈ R ntnn , where s = [s(t 1 ) T , . . . , s(t nt ) T ] T satisfy the following linear constraints:
(27b)
Function f ATD (·) can be rewritten in terms of these new variables as
We formulate the following non-convex Mixed Integer Nonlinear Program (MINLP) for optimal placement of pressure reducing valves and chlorine booster stations: 
Problem (29) has n t (8n p + 4n n + 2J + n 0 ) continuous variables, n t n p + 2n p + n n binary variables, and n t (2n p +J) non-convex quadratic terms.
The considered MINLP (29) combines difficulties in handling non-convex constraints with the presence of integer decision variables. In addition, the formulation of Problem (29) includes a discretised PDE for each network link, resulting in a large number of continuous variables and non-convex constraints, even for small-medium size WDNs -see Section 4. As shown in the experiments reported in Section 4, off-the-shelf solvers BARON and SCIP have failed to compute feasible solutions for various instances of Problem (29), and the resulting lower bounds to the optimal value are not tight. In this section, we propose a tailored algorithm for computing feasible solutions to Problem (29), and certified lower bounds to the optimal value. Note that the non-convex terms in (12) and (16) are the only non-linear terms within the formulation of Problem (29). Therefore, it is convenient to consider linear relaxations of these non-convex terms. First, we implement polyhedral relaxations of constraints (12). Let l ∈ P and k ∈ T , and define
. Consider the following constraints defined using linear over and under estimators of function φ l (·):
for all l ∈ P and k ∈ T . An example of these relaxations is shown in Figure 1 . The bilinear terms (16) are relaxed via the Reformulation Linearisation Technique (RLT) [SA99] . These relaxations are given by:
for all j ∈ {0, . . . , J l }, l ∈ P, and k ∈ T . Finally, observe that Problem (29) includes a large number of binary variables, which results in impractical computational effort for medium-large water networks. Therefore, we consider the following continuous polyhedral relaxation of Problem (29), where we also relax the binary constraints:
subject to (31), (32) 
The tightness of relaxations (31) and (32). depend on the bounds on the flow variables q l (t k ), l ∈ P, k ∈ T . Hence, we consider an optimisation based bound tightening (OBBT) scheme, to reduce the domain of the flow variables. A direct application of OBBT to Problem (29) would rely on the solution of a sequence of relaxed optimisation Problems (33), where the objective function is substituted by µq l (t k ), with µ ∈ {−1, 1}, l ∈ P, and k ∈ T . However, this would require a significant computational effort even for small-medium size water networks. We investigate an alternative approach, aimed at solving smaller linear relaxations of Problem (29). We expect flow variables to be primarily influenced by hydraulic variables and constraints in Problem (33). Hence, we consider a relaxation of Problem (33) where water quality related variables and constraints (17) 
where µ ∈ {−1, 1}. A further relaxation is obtained considering only a single time index in constraints (9), (11), (13), (14), and (31), which are all separable with respect to time indices. The proposed OBBT scheme relies on the solution of 2n t n p linear programs of the form:
subject to (9), (11), (13), (14), (31) only for
where µ ∈ {−1, 1}, l ∈ P, and k ∈ T . We propose the following method to compute feasible solutions for Problem (29), together with bounds on their level of sub-optimality:
Step 1. Implement Algorithm 1 to tighten the bounds on flow variables and compute lower bounds to the optimal value of Problem (29). The algorithm repeatedly performs OBBT until the computed lower bounds for Problem (29) np+l is considered for each link l ∈ P. We define a vectorv ∈ R 2np asv l := 1 if l ∈ J nv 0 otherwise.
For all i ∈ {1, . . . , iter}, solve the following nonlinear program using a local optimisation solver, e.g. IPOPT [WB06] :
Let (q best ,ĥ best ,η best ,θ best ,v best ) be the solution corresponding to the lowest AZP value among those computed. If none of the considered NLPs is feasible, then do not go to Step 3 and return LB best and UB best := +∞, with no feasible solutions.
Step 3. Defineû best := 10 3 |q best | andẑ best ∈ R ntnp such that
Solve the following MILP to optimally locate chlorine booster stations, minimising ATD: minimise z,c,r,w,y,v b ,s k∈T i∈Nd 
The solution of Problem (39) corresponds to optimal booster station locations and operational settings, which minimise ATD for the given pipe velocities.
Let (ĉ best ,r best ,ŵ best ,ŷ best , (v b ) best ,ŝ best ) be the computed solution. Define
The outlined method terminates with a lower bound LB best and an upper bound UB best to the optimal value of Problem (29). A worst-case estimate on the degree of sub-optimality of the computed solution is given by:
We evaluate the developed methods on different benchmark water distribution network models, with varying size and level of connectivity. All LPs and MILPs are solved using GUROBI (v9.0) [Gur20] , while the nonlinear programs are solved using IPOPT (v3.12.9) [WB06] . We consider a published benchmark network, referred to as 2loopsNet [OS93] . In addition, we formulate and solve the problem of optimal valve and chlorine booster placement using pescara and modena, originally presented by [BDL + 12]. In order to obtain more realistic problem instances, we have introduced temporal and spatial variability of demand profiles and Hazen-Williams roughness coefficients, respectively. We have also added first-order chlorine decay coefficients to network pipes. All case study models consider 24 hours of network operation, with a time step of one hour (i.e. n t = 24). Network hydraulic models and bounds on hydraulic heads and pipe flows are provided at [PSO20] . The layout of 2loopsNet is presented in Figure 2a . The network has n n = 6 demand nodes, n p = 10 links, and n 0 = 3 water inlets. Case study pescara include n n = 68 nodes and n p = 99 links, and n 0 = 3 water inlets -see Figure 2b . Finally, for modena we have n n = 268, n p = 317, and n 0 = 4 -see Figure 2c . Note that the considered case study networks result in large non-convex MINLPs, with a significant number of binary variables and non-convex termssee Table 1 . using the software for hydraulic and water quality analysis EPANET [RB96] , with fixed chlorine concentrations at inlets equal to 0.5 mg/l. We define c 0 in equation (2) as the nodal concentrations at 24:00 hours computed by the EPANET simulation. In the formulation of Problem (29), maximum allowed concentration at demand nodes is set 2 mg/l, while chlorine concentrations at network inlets are not allowed to be greater than 0.5 mg/l. The target concentration at demand nodes is 1 mg/l. Finally, in the discretisation of the PDE (1), we set a temporal time step ∆t = 3600s (1 hour) and ∆x l = L l 2 , where L l is the length of link l, for all l = 1, . . . , n p . We set m = 5 in (31) as preliminary experiments have shown that this setting results in sufficiently tight polyhedral relaxations, for all case studies -see also Figure 1 .
We have implemented the off-the-shelf solvers BARON (v1.89) [KS17] and SCIP (v5.0.1) [GEG + 17] and compare their performance with the method proposed in Section 3. We consider Problem (29) with n v = 2 and n b = 2, as these are among the problem instances where the proposed method results in larger relative optimality gaps. For the off-the-shelf solvers, we set a maximum time limit of one day (86400 s). As reported in Table 2 , BARON and SCIP have failed to compute a feasible solution for all three case studies. Moreover, SCIP failed to compute a lower bound in the case of modena. Figure 3 compares the lower bounds computed at each iteration of Algorithm1, with those obtained by off-the-shelf solvers after a day of computations. We observe that the lower bounds computed by Algorithm 1 are tighter than those computed by off-the-shelf solvers, for both pescara and modena. Figure 3 also shows that the implemented bound tightening algorithm improves the computed lower bounds.
We formulate Problem (29) for n v ∈ {1, 2, 3} and n b ∈ {0, . . . , 3} in 2loopsNet, and n v ∈ {1, . . . , 5} and n b ∈ {0, . . . , 5} in pescara and modena. Hence, for 2loopsNet we consider 20 different problem instances, while 30 different problems are formulated for each case study network pescara and modena. The solution method described in Section 3 is implemented to compute feasible solutions for the considered problem instances, together with bounds on their levels of sub-optimality. The full results are reported in Tables 1-3 Figure 4 . In all instances, the implemented algorithm computed bounds on the level of sub-optimality smaller than 7% for 2loopsNet, while the worst-case sub-optimality gaps obtained for pescara and modena are larger. As expected, the computational time increases with the size of the considered case study. This is mainly due to the bound tightening scheme, where two linear programs are solved for all network links and times steps, at each iteration. More computationally efficient implementations could reduce the time spent in bound tightening, for example by opportunely executing some computations in parallel, or exploiting network structure to propagate bounds on particular links without solving an optimisation problem.
In Figure 5 , we report the computed AZP values for the feasible solution obtained by the proposed heuristic. For the same number of installed valves n v , the computed AZP values for n b = 0, . . . , 5 are the same. As it should be expected, feasible solutions computed for increasing number of valves correspond to decreasing values of AZP.
Analogously, the average target deviation for nodal concentrations is reduced as additional chlorine booster stations are installed -see Figure 6 . Without any chlorine booster station, there is limited ability to control chlorine concentrations, using only the injected concentrations at water sources, which, in our formulation, can not be greater than 0.5 (mg/l). As we install more booster stations, the system is able to maintain nodal concentrations closer to the target.
Observe that the best possible value of ATD is 0. However, several hours are required for nodal concentrations to reach the optimised level, following the installation and operation of chlorine booster stations. If the travel time between a newly installed booster station and a specific node is T age hours, we expect nodal concentrations to reflect the action of the booster station after T age hours. In addition, network topology and spatial distribution of decay coefficients can affect the ability to control chlorine concentrations at selected locations. Therefore, we do not expect nodal concentrations to be exactly equal to the target at all time steps. However, these are very close to the target concentration after they reach their optimised level. This is shown in Figure 7 , which reports the chlorine residual concentrations at two nodes in pescara and modena, experiencing the largest average target deviation in most problem instances. 
We have proposed a new mixed integer nonlinear programming formulation for the problem of optimal placement and operation of pressure reducing valves and chlorine booster stations in water distribution networks. The numerical experiments reported in this manuscript show that off-theshelf global optimisation solvers fail to compute feasible solutions for the considered problem, and the computed lower bounds to the optimal value are not tight. We have implemented polyhedral relaxations and a bound tightening scheme resulting in improved lower bounds compared to off-theshelf solvers, which rely on spatial branch and bound algorithms. Furthermore, we have proposed a tailored heuristic method to compute feasible solutions for the considered problem. The developed algorithms have been evaluated by solving multiple problem instances for three case study networks. Tightened relaxations and heuristic method are shown to outperform off-the-shelf solvers for the considered case studies. In addition, the developed heuristic has enabled the computation of good quality feasible solutions for the vast majority of the considered problem instances.
The proposed method enables the joint optimisation of pressure and disinfectant dosage in water distribution networks. This allows water utilities to implement integrated and efficient schemes for pressure and water quality management, in order to minimise leakage and protect public health. Future work should extend the problem formulation to include the operation of pumps and water tanks within the same optimisation framework. Moreover, the proposed polyhedral relaxations could be tightened, for example implementing semidefinite or second-order cone relaxations of the non-convex quadratic constraints. 