Optimal design-for-control of self-cleaning water distribution networks using a convex multi-start algorithm

The provision of self-cleaning velocities has been shown to reduce the risk of discolouration in water distribution networks (WDNs). Despite these findings, control implementations continue to be focused primarily on pressure and leakage management. This paper considers the control of diurnal flow velocities to maximize the self-cleaning capacity (SCC) of WDNs. We formulate a new optimal design-for-control problem where locations and operational settings of pressure control and automatic flushing valves are jointly optimized. The problem formulation includes a nonconvex objective function, nonconvex hydraulic conservation law constraints, and binary variables for modelling valve placement, resulting in a nonconvex mixed integer nonlinear programming (MINLP) optimization problem. Considering the challenges with solving nonconvex MINLP problems, we propose a heuristic algorithm which combines convex relaxations (with domain reduction), a randomization technique, and a multi-start strategy to compute feasible solutions. We evaluate the proposed algorithm on case study networks with varying size and degrees of complexity, including a large-scale operational network in the UK. The convex multi-start algorithm is shown to be a more robust solution method compared to an off-the-shelf genetic algorithm, finding good-quality feasible solutions to all design-for-control numerical experiments. Moreover, we demonstrate the implemented multi-start strategy to be a fast and scalable method for computing feasible solutions to the nonlinear SCC control problem. The proposed method extends the control capabilities and benefits of dynamically adaptive networks to improve water quality in WDNs.


Introduction
The management of water quality in water distribution networks (WDNs) presents a complex operational challenge. As a direct consequence of ageing and deteriorating infrastructure, the mitigation of discolouration incidents is becoming one of the key operational challenges within water quality management programmes. In addition to discolouration being the largest source of customer complaints (Vreeburg and Boxall, 2007, Husband and Boxall, 2011, Armand et al., 2017, there is a growing body of evidence suggesting its occurrence harbours increased microbial activity (Liu et al., 2013, 2014, van der Wielen and Lut, 2016. These conditions can accelerate biofilm growth and drastically reduce the efficacy of disinfectant residuals in protecting against waterborne illness and contaminants intrusion. Moreover, WDNs in the UK are highly sectorized and operated with fixed topology for purposes of leakage management. This type of network configuration, referred to as district metered areas (DMAs), has been demonstrated to exacerbate water quality deterioration and increase the risk of discolouration incidents (Machell andBoxall, 2014, Armand et al., 2018). With progressively stringent water quality regulations, water companies are seeking effective and cost-efficient operational control strategies to reduce the risk of discolouration.
Discolouration is primarily a consequence of resuspended material accumulated within WDNs (Vreeburg and Boxall, 2007). It can materialize from the cumulative impact of the following processes (Boxall and Dewis, 2005): (i) the ingress and/or development of particulate matter; (ii) the accumulation of particulates at the pipe invert and/or formation of cohesive layers at the pipe wall; and (iii) a hydraulic disturbance (i.e. trigger event), which mobilizes loose particulates and generates sufficient shear stress to overcome cohesive forces at the pipe wall. Such hydraulic disturbances can be generated from different phenomena, including pressure transients during unsteady hydraulic conditions (Aisopou et al., 2012). Apart from their origin, the physical pathways of discolouration are intrinsically connected to network hydraulics. In a recent study focusing on the impact of network sectorization on water quality, Armand et al. (2018) proposed a set of surrogate hydraulic variables for discolouration risk assessment. Central to their findings was the role of diurnal flow velocities on particle transport and fate. This connection between discolouration and hydrodynamic conditions has been supported by numerous experimental and theoretical studies; see van Summeren and Blokker (2017) and Armand et al. (2018) for reviews on the topic. These studies have mainly focused on the development of predictive tools for modelling particle transport and accumulation processes. Most notably, Boxall et al. (2001) developed the Prediction of Discolouration in Distribution Systems (PODDS) model, an empirically-based numerical tool which aims to characterize cohesive layer strength at the pipe wall. The PODDS model was later updated to account for material regeneration in Furnass et al. (2014), where both erosion and regeneration processes require calibration using continuous flow and turbidity data. Because such tools require extensive field testing and are generally limited to pipe-level assessments, their use in practice has not yet been widespread. Recognizing this limitation, van Summeren and Blokker (2017) presented a theoretical particle transport model, combining the effects of gravitational settling, hydraulic shear stresses, and bed-load transport. To complement this, several laboratory-based experimental studies have emerged to better understand the complex interactions between particle properties and pipe hydraulics (e.g. Sharpe et al., 2019;Braga et al., 2020).
In addition to predictive modelling, research has also focused on reducing the severity and frequency of discolouration incidents through network design, maintenance, and control. Water companies in the Netherlands have been conducting experimental research on the design and implementation of controls for self-cleaning networks. The self-cleaning capacity (SCC) of a WDN is defined as the ability for pipes to experience peak daily flow velocities above a threshold required to routinely re-suspend particles and thus prevent accumulation (Vreeburg et al., 2009). Previous experimental programmes have suggested resuspension velocities on the order of 0.2 m/s to 0.25 m/s in distribution pipes (Ryan et al., 2008, Blokker et al., 2010. This has been corroborated with a recent field study monitoring turbidity under various flow rates, where an increase in turbidity levels were observed at flow velocities greater than 0.2 m/s (Prest et al., 2021). Water companies in the Netherlands have demonstrated successful self-cleaning implementations by redesigning looped, oversized networks to branched layouts with smaller diameter pipes (Vreeburg et al., 2009). A recent study has also investigated the trade-off between self-cleaning velocities and fire flow capacity in North American WDNs (Gibson et al., 2019). However, since the redesign of WDN infrastructure becomes cost-prohibitive at scale, there have been recent forays in the reconfiguration of existing network topology to promote self-cleaning networks (Blokker et al., 2012, Abraham et al., 2016, 2018. Combining UK and Dutch experience, Abraham et al. (2016Abraham et al. ( , 2018 formulated an optimization problem for increasing SCC by redistributing flow through changes in network topology. More specifically, the optimization problem aimed to maximize the number of pipes with flow velocities above a self-cleaning threshold through two separate strategies: (i) optimal closure of isolation valves and (ii) optimal operational settings of existing pressure control valves (Abraham et al., 2016(Abraham et al., , 2018. Abraham et al. (2018) solved the problem of optimizing valve closures using a linear graph analysis tool. Following Schaub et al. (2014), a line-outage distribution factor (LODF) matrix was computed to estimate the flow redistribution resulting from an outage (closure) of an (or multiple) edge-to-edge relation(s). For the optimal control problem, Abraham et al. (2016) computed a local solution by approximating the nonsmooth objective function as a continuous nonlinear function, followed by application of a tailored sequential convex programming algorithm. While the benefits of the LODF solution method for optimal valve closures were demonstrated numerically using an operational network in the Netherlands, results from the control problem were limited to a small-scale theoretical network. Moreover, decision variables were restricted to the control of existing unidirectional pressure reducing valves (PRVs). Building on the SCC optimization problem posed in Abraham et al. (2016Abraham et al. ( , 2018, this manuscript considers both control and design-for-control problem formulations. The latter involves the simultaneous optimization of valve placement and operational settings for both existing and new control valves. In addition to unidirectional PRVs, this work also considers bidirectional dynamic boundary valves (DBVs) and automatic flushing valves (AFVs) as dynamic hydraulic controls. These hydraulic controls were developed to facilitate the novel operational framework of dynamically adaptive networks (Wright et al., 2014, Ulusoy et al., 2022. The resulting optimization problem is formulated as a nonconvex mixed integer nonlinear program (MINLP).
Both mathematical optimization and heuristic methods have been used to solve design and control problems in WDNs (see literature review in Mala-Jetmarova et al., 2017). For mathematical optimization methods, scalability is recognized as a current limitation in solving MINLP problems to global optimality (Koch et al., 2012, Sahinidis, 2019; that is, the implementation of global solvers become impractical for large problem cases. Consequently, heuristic approaches are often employed to compute satisfactory feasible solutions. A common heuristic method used for WDN optimization problems is the genetic algorithm (GA). While GAs have been successfully applied to design problems, the computational effort required to find solutions sufficiently close to the global optimum grows rapidly with problem size (Maier et al., 2014). In this manuscript, we develop a heuristic algorithm based on convex optimization and a multi-start scheme to compute feasible solutions to the considered MINLP problem. To handle integer variables, we first formulate a convex subproblem through polyhedral relaxations of nonconvex terms and the continuous relaxation of binary variables. We subsequently employ a randomization heuristic to sample N candidate valve configurations from the set of fractional values generated from the convex subproblem. We then fix binary variables for each sampled valve configuration and compute (locally) optimal operational settings from a nonlinear programming (NLP) control problem. This follows the heuristic algorithm presented in Pecci et al. (2022), extending its application to the SCC design-for-control problem and to the nonsmooth Hazen-Williams friction model. Since the degree of nonlinearity of the SCC problem is higher than the problem investigated in Pecci et al. (2022), we include a multi-start strategy and a feasibility restoration problem for selecting starting points. This step aims to minimize the risk of poor local optima as well as ensure hydraulic feasibility of the NLP control problem. Finally, the best feasible solution is selected from the set of sampled valve configurations. The proposed heuristic algorithm further increases the benefits from the implementation of dynamically adaptive networks as it expands their control capabilities to enhance water quality in WDNs.
This manuscript is organized as follows. In Section 2, we formulate the design-for-control problem of maximizing the network SCC through dynamic hydraulic controls. We then present the proposed heuristic algorithm in Section 3. Finally, in Section 4, we demonstrate the performance of the developed heuristic algorithm using three case study networks with varying size and degrees of complexity. To facilitate a broader discussion on heuristic approaches for the design and control of WDNs, we compare the results with an off-the-shelf GA implementation, which is a common approach used in the literature.

Problem formulation
We investigate a design-for-control problem to maximize the length of network pipes experiencing flow velocities above a given self-cleaning capacity (SCC) threshold. This is achieved by installing new valves and/or controlling their operational settings. For this purpose, our problem formulation considers three valve types as pressure and connectivity control actuators. First, pressure reducing valves (PRVs), which are modelled having unidirectional flow. Second, bidirectional dynamic boundary valves (DBVs), for which flow is permitted in both directions across discrete model time steps. Here, DBVs represent the operation of remote-controlled isolation valves, which modulate flow and pressure between adjacent zones. Third, automatic flushing valves (AFVs), whose flushing rate is bounded by a set maximum value. Throughout this manuscript, we refer to either PRVs or DBVs as control valves, as both have the capability of controlling pressure and are modelled at network links. On the other hand, AFVs are simply referred to as flushing valves and are modelled at network nodes. We consider operational scenarios, for which PRV locations have been fixed to minimize average zone pressure (AZP), and thus decision variables include only their operational settings. In comparison, both locations and operational settings of DBVs and AFVs are considered as decision variables. The operational settings of valves are modelled as continuous variables, whereas their placement (location) are modelled through binary variables. All network links and nodes are considered as potential locations of DBVs and AFVs, respectively. As the current stage of this work focuses on the selfcleaning capacity of pipes at the DMA or distribution level, we do not consider storage tanks or pumping activity as forms of hydraulic control. Therefore, we assume discrete and hydraulically independent model time steps. Finally, it is noted that this work relies on the availability of a calibrated hydraulic model.

Hydraulic variables and constraints
The problem considers a water distribution network (WDN) with n p links, n n demand nodes and n 0 known head nodes (e.g. water sources, reservoirs). The network is modelled as a directed graph with n p edges (links) and n n + n 0 vertices (nodes). A demand-driven hydraulic analysis is used to simulate steady-state network hydraulics over n t discrete time steps. For each time step t ∈ {1, . . . , n t }, known hydraulic conditions are given by vectors of nodal demands d t ∈ R nn and source hydraulic heads h 0 t ∈ R n0 . Moreover, vectors η t ∈ R np and α t ∈ R nn are included to model local losses introduced by the action of control valves and operational demands at flushing valves, respectively. Unique vectors of hydraulic states q t ∈ R np and h t ∈ R nn are computed by solving the following steady-state energy (1a) and mass (1b) conservation equations governing pipe flow: where A 12 ∈ R np×nn and A 10 ∈ R np×n0 are the link-node incidence matrices for demand and known head nodes, respectively; and the vector φ(q t ) = [φ 1 (q 1,t ) . . . φ np (q np,t )] T models frictional head losses associated with flows q t . Omitting time index t, φ j (q j ) is defined in general form for flow conveyed across link j as: where the resistance coefficient r j and exponent n j , both independent of time t, take different values depending on the link type (e.g. pipe or valve) and on the frictional head loss model. For valve links, n j = 2 and r j = 8Kj gπ 2 D 4 j , with K j and D j representing the valve loss coefficient and diameter, respectively (Larock et al., 1999). In this work, we apply the Hazen-Williams (HW) model to characterize frictional head losses across pipe links. The HW model is an explicit and empirical relationship between pipe flow and frictional head loss, with n j = 1.852 and r j is defined as follows for all j ∈ {1, . . . , n p }: where C is the HW coefficient, a dimensionless number representing frictional characteristics; L is pipe length in meters; and D is pipe diameter in meters (Larock et al., 1999). Similarly, explicit approximations of the Darcy-Weisbach formula (e.g. Valiantzas, 2008) could be used to model frictional head losses. Finally, it is convenient to isolate the nonlinear term φ(q t ) in (1a). Here, we introduce a vector of auxiliary variables θ t ∈ R np , which separates the energy conservation constraint into its linear and nonlinear components, as follows: Valve placement and operation are modelled as follows. For each time step t ∈ {1, . . . , n t }, the continuous variable η t ∈ R np presented in (1a) models the local losses introduced by the action of control valves and the continuous variable α t ∈ R nn presented in (1b) models the flow emitted at flushing valves. Moreover, binary variables z ∈ {0, 1} np are included to model PRV and DBV placement, and v + t ∈ {0, 1} np and v − t ∈ {0, 1} np to assign their control capabilities in the positive or negative flow direction, respectively, across each time step t. Thus, for all links j ∈ {1, . . . , n p } and time steps t ∈ {1, . . . , n t }, binary variables z j , v + j,t , and v − j,t are set as Analogously, the placement of AFVs at network nodes is modelled using binary variables y ∈ {0, 1} nn , defined as y i = 1 flushing valve placed at node i 0 no valve These binary variables are subject to the following physical and economical constraints, which limit pressure control capabilities in a single direction at each time step t (7a) and enforce a maximum number of control valves n v and flushing valves n f considered for installation (7b)-(7c): Since we assume that existing PRVs have fixed location and unidirectional flow, z j , v + j,t , and v − j,t are set a priori for all time steps t ∈ {1, . . . , n t } at the known set of PRV links N PRV ⊆ {1, . . . , n p }.
We introduce constant vectors to bound the continuous hydraulic variables and formulate big-M constraints for modelling the operation of control and flushing valves. For a given vector of maximum allowed velocities u max ∈ R np and vector of link cross-sectional areas A ∈ R np , let q L t = −Au max and q U t = Au max be the vectors of lower and upper bound flows across network links at time step t, respectively. Bounds on the auxiliary head loss vector θ t are set as θ L t := φ(q L t ) and θ U t := φ(q U t ). Moreover, let h min t and h max t ∈ R nn specify minimum and maximum heads at network nodes, respectively. The minimum head is set to a minimum regulatory pressure plus the node elevation and the maximum head is set to the largest available known source head. Bounds on η t for j ∈ {1, . . . , n p } are then defined as follows: We formulate big-M constraints to model valve placement and enforce energy conservation at control valve links, ensuring η t and q t act in the same direction. These constraints are written as follows: Finally, lower and upper bounds on hydraulic variables h t , q t , η t and θ t are set to define the feasible solution space, 2.2. Self-cleaning capacity objective function The objective of this study is to maximize the length of network pipes satisfying the self-cleaning capacity (SCC) threshold. The SCC objective function is defined as the following length-weighted sum over all pipes n p and hydraulic time steps n t (Abraham et al., 2018): where A is the link cross-sectional area; and κ j (·) is an indicator function which models the state of pipe velocities with reference to a minimum threshold. The indicator function is described for link j as follows with u min j representing the threshold flow velocity at link j, defined a priori (see presented literature in Section 1). Moreover, a weighting is included to normalize the length of link j to the entire network, w j = Lj np k=1 L k , where L ∈ R np is the vector of pipe lengths.
The SCC objective function f SCC is nonsmooth at ±u min , resulting in unbounded gradients. Therefore, in order to employ gradient-based optimization methods, the threshold function κ(·) is approximated with a continuous sum of sigmoids (or logistic) function, as proposed in Abraham et al. (2016). The sigmoidal function has positive and negative components, defined by ψ + j (u) : , respectively, where ρ is a parameter which sets the sigmoid function curvature. The following expression combines these sigmoid functions to approximate f SCC posed in (12): where velocity is defined as u j,t = qj,t Aj . In accordance with that reported in Abraham et al. (2016), we found ρ ≤ 100 provided a step-like objective function, whilst still being sufficiently smooth at the threshold boundaries for gradients to exist. An example of the SCC indicator function κ(·) and its continuous sum of sigmoids approximation f SCC are shown in Figure 1.

Mixed-integer nonlinear program
The SCC design-for-control problem aims to maximize (14), subject to hydraulic conservation laws and physical and economical valve constraints. The problem formulation includes continuous variables denoted by x := [q h η θ α] T and binary variables denoted by v := [v + v − ] T , y and z. Here, the objective function is replaced with its additive inverse to result in a minimization optimization problem. The resulting mixed integer nonlinear program (MINLP) is summarized by the following problem formulation. (14) subject to linear hydraulic conservation constraints (1b) and (4a) nonconvex HW head loss model constraints (4b) big-M constraints for control valve operation (9a) − (9f) big-M constraints for flushing valve operation (10) physical and economical valve constraints where Q is a rectangle representing upper and lower bounds for the continuous decision variables. The continuous decision variables are defined as: q := (q t ) t=1,...,nt , h := (h t ) t=1,...,nt , η := (η t ) t=1,...,nt , θ := (θ t ) t=1,...,nt and α := (α t ) t=1,...,nt . The binary decision variables varying with time step t are defined as: ..,nt . Problem (MINLP) has n t (3n p + 2n n ) continuous variables, 2n t n p + n p + nn binary variables and 2n t n p nonconvex terms. Observe that the problem grows rapidly with the size of the considered WDN (see Table 1), making it a difficult nonconvex MINLP problem to solve. To overcome these challenges, we develop a convex heuristic to compute feasible solutions to Problem (MINLP). The following section describes the solution algorithm and its implementation details.

Solution method
The proposed solution algorithm combines convex relaxations with a randomization heuristic and multistart solver to compute feasible solutions to Problem (MINLP). First, we formulate convex relaxations of Problem (MINLP), which yield a linear programming (LP) subproblem. Here, the nonconvex SCC objective function (12) and nonconvex energy conservation constraints (4b) are relaxed using polyhedral envelopes, and a continuous relaxation is applied to binary decision variables. We also include a domain reduction step where the resulting convex subproblem is tightened using both model decomposition and optimization-based bound tightening (OBBT) techniques. Then, a randomization heuristic uses the fractional valve placement values from the convex subproblem solution to form nonlinear programming (NLP) control problems. Local solutions to these NLP problems are computed using a strictly feasible sequential convex programming (SFSCP) solver. We also implement a multi-start strategy, which includes an optimization-based feasibility restoration problem to ensure hydraulic feasibility of the starting points. Figure 2 offers a detailed overview of the solution method, referred to as the convex multi-start (CMS) algorithm. The algorithm steps and overall implementation details are provided in the following subsections.

Convex relaxation
There are two sources of nonlinear nonconvexity that make Problem (MINLP) difficult to solve. These are the SCC objective function f SCC (14) and the HW head loss model φ(·) in equality constraint (4b). One approach to overcome nonconvexity is through polyhedral relaxations, which can be formulated as linear constraints and thus efficiently handled by state-of-the-art linear solvers. Since the resulting mixed integer Initialization Hydraulic model and SCC parameters

Output
Best feasible solution z' best ,y' best ,η best ,α best z' i , y' i , η

Multi-start solver Feasibility restoration and NLP control problems
Store solution and repeat for i = 1,...,N trials linear program (MILP) may still have a large number of binary variables, it is convenient to apply a continuous relaxation to binary variables. In addition to the computational advantages, a continuous relaxation increases the search space for the set of optimal binary variables. Here, we implement the aforementioned relaxation techniques to formulate a convex subproblem of Problem (MINLP), with its solution forming the basis of the heuristic algorithm.
We first relax the nonconvex objective function in (14) with a linear outer approximation. Let σ + t ∈ R np and σ − t ∈ R np be vectors of auxiliary variables introduced to model the positive ψ + and negative ψ − sigmoid functions, respectively. The objective function f SCC is then reformulated as the following set of inequality constraints: To ensure equivalence with (14), the SCC objective function becomes We construct concave envelopes for the positive ψ + and negative ψ − components of f SCC . This follows the methodology presented in Boyd (2014, 2016) for sigmoidal functions, resulting in piecewise linear relaxations. These relaxations are written as the following constraint: as well as the sigmoid function parameters. A detailed derivation of these relaxations is provided in Appendix A.1. Moreover, an example of the implemented relaxation is illustrated in Figure 3a, withψ denoting the set of linear relaxations.
We then implement polyhedral relaxations for the HW head loss model φ(·) in equality constraint (4b). This builds on the relaxation methods for monomials of odd degree introduced by Liberti and Pantelides (2003) by formulating polyhedral relaxations for the HW head loss model. Similar to (A.23), the formulated relaxations are written as the following linear constraint: where matrices R t and E t and vector r t are derived from flow bounds q L t and q U t and the HW model parameters. Further details are presented in Appendix A.2 and an example of the implemented relaxation is illustrated in Figure 3b, noting that the polyhedral relaxation is denoted byφ.
Lastly, we implement the continuous relaxation of binary variables v, y, and z. These are combined with continuous decision variables x, defined in Problem (MINLP), and the continuous auxiliary variables associated with the f SCC inequality constraints (A.22a) and (A.22b), denoted by σ := [σ + σ − ]. The resulting convex relaxation of Problem (MINLP) is represented by the following LP subproblem. minimize x, v, y, z, σ − f SCC reformulated in (16) subject to linear relaxations of HW head loss model constraints (18) linear relaxations of ψ + and ψ − sigmoid functions (A.23) linear hydraulic conservation constraints (1b) and (4a) big-M constraints for control valve operation (9a) − (9f) big-M constraints for flushing valve operation (10) physical and economical valve constraints (7a) − (7c) The optimal value to Subproblem (LP) yields a lower bound to the original nonconvex Problem (MINLP).
In particular, vectors y ∈ R nn and z ∈ R np of continuous variables for valve placement can be interpreted as the probability distribution from which random samples are drawn to solve the NLP control problem.
Moreover, the vector η ∈ R np from Subproblem (LP) is used as one of M starting points in the multi-start strategy. These steps are discussed in the subsequent sections.
Furthermore, we implement a domain reduction procedure to reduce the bound intervals and thus strengthen the convex relaxations formulated for Subproblem (LP). The main procedure implemented is an optimizationbased bound tightening (OBBT) algorithm (Belotti et al., 2009). The OBBT algorithm solves a series of LP optimization problems, with objective functions set to both maximize and minimize link flow. This is augmented with a forest-core decomposition scheme to reduce the number of flow variables whose bounds are tightened (Simpson et al., 2014). The forest of a WDN comprises the disjoint union of all outer branch (or tree) components (Deuerlein, 2008). The core represents the set of looped (or block) graphs, which contain the roots of all forest trees. Following Pecci et al. (2019), we only perform bound tightening on the set of core links. We include psuedocode for the OBBT algorithm in Appendix B.

Randomization heuristic
In searching for a good quality local solution to Problem (MINLP), we employ a randomization heuristic to sample candidate valve configurations. Here, the fractional values of y ∈ R nn and z ∈ R np yielded from the solution to Subproblem (LP) define a discrete probability distribution for the index sets {1, . . . , n n } and {1, . . . , n p }, respectively. We propose to randomly sample candidate pressure control n v and flushing n f valve configurations from the respective probability distributions over N sampling trials. This creates sampled vectors of binary variables y i=1,...,N ∈ {0, 1} nn and z i=1,...,N ∈ {0, 1} np , from which a local solution to the NLP control problem is obtained for each trial i ∈ {1, . . . , N } (see Section 3.3). To avoid redundant valve configurations, we store binary values from each trial in the set P ∈ R N ×(nv+n f ) , checking its intersection with {y i ∪ z i } before fixing valve placement values and proceeding to the NLP control problem. Psuedocode for the implemented randomization heuristic is detailed in Algorithm 1.

Multi-start solver
The final component of the CMS solution algorithm concerns the optimization of valve settings. This control problem is formulated to compute (locally) optimal operational settings for each candidate valve con-Algorithm 1 Randomization heuristic 1: Input: vectors of fractional values y ∈ [0, 1] nn and z ∈ [0, 1] np Problem (LP) 2: Output: sampled vectors of binary variables y i=1,...,N ∈ {0, 1} nn and z i=1,...,N ∈ {0, 1} np 3: Initialize P ← ∅ 4: for i = 1, . . . , N do 5: Sample y i from vector y with probability weights equal to fractional values 7: Sample z i from vector z with probability weights equal to fractional values 8: end while 9: Update visited valve locations: Vectors y i and z i Input to multi-start solver 11: end for figuration generated from the randomization heuristic. Because the objective function in Problem (MINLP) is highly nonlinear, we implement a strategy to use multiple starting points in order to avoid getting trapped in poor local optima. This includes a feasibility restoration subproblem to ensure hydraulically feasible starting points are passed to the NLP solver. We refer to the overall solution process as the multi-start solver, which has steps highlighted by the process diagram shown in  (14) subject to hydraulic conservation constraints (1a) and (1b) x ∈Q (NLP) whereQ represents the upper and lower variable bounds modified to account for directional constraints enforced by the fixed binary variables. Since the focus of this study is at the DMA-or distribution-level, which in practice limits the number of installed DBVs for consideration, we use an enumeration approach to test all DBV flow direction combinations for v + t and v − t . This approach was deemed appropriate on the basis of the relatively small computational (CPU) times recorded in Section 4. For example, the numerical experiment corresponding to the maximum number of installed DBVs on the largest network tested in this study (see Section 4.1 for details) resulted in a CPU time of less than five minutes. This CPU time is well within an assumed one hour limit for application in control strategies. We acknowledge, however, that this assumption warrants re-evaluation for system-wide studies which inevitably consider a larger set of control valves. Here, a total of T = 2 nDBV control problems are solved to find the best DBV direction combination at each time step t ∈ {1, . . . , n t } (see Figure 4). In this work, we use the strictly feasible sequential convex programming (SFSCP) method outlined in Wright et al. (2015) for solving a local optima to Problem (NLP). The SFSCP solver, also used to optimize SCC control valve settings in Abraham et al. (2016), is implemented because of its reliable convergence properties, namely its strict feasibility requirement for each iterate of the optimization process. The SFSCP solver linearizes the objective function and nonlinear hydraulic conservation constraints in Problem (NLP), from which a sequence of linear programs can be efficiently solved. Linear approximations are formulated using first-order Taylor expansions around the point x k := [q k h k η k α k ] T . Let f := f SCC and the hydraulic conservation equality constraints (1a) and (1b) be denoted by g(·). It follows that

Starting points
and are the linearized expressions used to formulate a convex subproblem for each iterate k. This subproblem computes step directions dη k and dα k at point x k for pressure control and flushing valve operational variables, respectively. Following Wright et al. (2015), a strictly feasible line search method is implemented to guarantee an improvement in objective function whilst ensuring hydraulic feasibility of the computed local optima. We define the hydraulically feasible region by the set F := x ∈ R np+nn+nv+n f x satisfies constraints (1a), (1b), and variable boundsQ . (21) The null space solver proposed in Abraham and Stoianov (2015) is used to solve hydraulic states q k+1 and h k+1 for the optimal controls obtained at iterate k. If x k ∈ F, then the optimization process moves to the next iterate k := k + 1. Otherwise, a backtracking procedure is implemented to reduce dη k and dα k until hydraulic feasibility is achieved. Termination of the SFSCP solver is triggered when relative improvements in the objective function are below a set tolerance tol . The SFSCP solver steps are detailed in psuedocode presented in Appendix C (Wright et al., 2015). Moreover, the subsequent section describes a multi-start strategy for selecting feasible starting points x 0 to be passed to the SFSCP solver. For completeness, we provide a brief discussion on the performance comparison with state-of-the-art nonlinear solver IPOPT. This is included as many large-scale WDN optimization studies have used IPOPT to compute locally optimal solutions for network controls (e.g. Pecci et al., 2019;Ulusoy et al., 2020;Nerantzis and Stoianov, 2022). Most of these studies, however, have applied a quadratic approximation (QA) to model the head loss function ψ(·) as the nonsmooth HW formula presents challenges for IPOPT, namely that its Hessian is unbounded at the origin. Preliminary efforts to introduce a QA model for the current SCC problem resulted in large inaccuracies compared to the HW model. This was expected as the current study aims to maximize pipe flow velocities, which increases the QA model error since a larger maximum flow value is required (Pecci et al., 2017;Equation 21). Therefore, we tested IPOPT using the HW model and relied on its feasibility restoration scheme to mitigate potential errors. Infeasibility errors were observed for the larger case study network when using the IPOPT solver, ultimately leading to the SFSCP solver as the more reliable and robust option for the current SCC control problem. The SFSCP and IPOPT solver comparison is detailed in Appendix D.

Multi-start strategy
As highlighted in Figure 4, a multi-start strategy is implemented to find good quality local optima to Problem (NLP). Here, we initialize starting points x 0 for multiple initial conditions j ∈ {1, . . . , M }, which are sequentially passed to the SFSCP solver. Starting points are formed by simulating hydraulic states for different initial control valve settings η 0 . For initial condition j = 1, we set η 0 j=1 equal to the control valve settings resulting from Subproblem (LP). Otherwise, η 0 j=2,...,M is randomly generated between the previously set lower and upper η bounds. Note that initial flushing operations are set as α 0 = {0} for all initial conditions j ∈ {1, . . . , M }. Since the initial control conditions do not guarantee a solution in the hydraulically feasible region, we also include a feasibility restoration step. That is, if starting point x 0 j / ∈ F for all j ∈ {1, . . . , M }, feasibility is restored through the solution to the following NLP optimization problem. minimize x η − η 0 2 2 subject to hydraulic conservation constraints (1a) and (1b) x ∈Q where the objective function is set to minimize the squared 2 -norm between the vector of initial settings η 0 and settings η which produce a hydraulically feasible starting point. Problem (FR) is solved using the nonlinear IPOPT solver. Apart from cases where the fixed DBV directions strictly cannot facilitate a hydraulically feasible solution, we did not encounter infeasibility errors when using IPOPT as the NLP solver for Problem (FR). The discussion in Appendix D suggests that the nonconvexity of the SCC objective function in Problem (NLP) is a key factor for the observed infeasibility issues.

Numerical experiments and discussion
We evaluate the SCC design-for-control problem using the proposed CMS solution algorithm. Since CMS is a heuristic, we also compare its performance with an off-the-shelf GA implementation, which is a common heuristic method for solving WDN design problems. This section is organized as follows. First, in Section 4.1 we describe the problem parameters, case study networks, and computational resources used to carry out the numerical experiments. Next, Section 4.2 discusses the results obtained when solving the control-only subproblem to Problem (MINLP), where valves locations are fixed. In this case, operational settings of existing PRV and DBV (if applicable) valve configurations are optimized to maximize selfcleaning conditions. We then investigate the solution to the design-for-control problem in Section 4.3, where the placement and operation of both pressure control and flushing valves are considered to improve SCC performance.

Computational setup
All numerical experiments were performed in MATLAB R2021b (64-bit) for Microsoft Windows 11, installed on a 2.50-GHz Intel(R) Core(TM) i9-11900H CPU with 8 cores and 32.0 GB of memory (RAM). The LP problems involved in the convex relaxation, OBBT algorithm, and SFSCP solver were solved using Gurobi (v9.5.0; Gurobi Optimization, 2022), accessed via its MATLAB interface. With the exception of large problem cases, which had solver parameters (e.g. NumericFocus, Crossover, and Presolve) modified to deal with numerical issues, Gurobi's default parameters were applied. Local solutions to the NLP feasibility restoration problems were solved using the state-of-the-art interior point solver IPOPT (v3.12.19;Wächter and Biegler, 2006). IPOPT was accessed in MATLAB via the OPTI toolbox interface (Currie et al., 2012) and was implemented passing the exact Jacobian and Hessian matrices as inputs to the solver. The single objective GA was implemented using MATLAB's off-the-shelf Global Optimization Toolbox with its default solver parameters (e.g. population size = 100; crossover fraction = 0.8; function tolerance = 1e−6). For each GA fitness evaluation, EPANET2.2 (Rossman et al., 2020) hydraulic solver was accessed via the EPANET-MATLAB Toolkit interface (Eliades et al., 2016). We set maximum time limits for the GA implementation associated with control and design-for-control problems of 6 and 12 hours, respectively. Furthermore, MATLAB's parallel computing toolbox was employed (with eight local workers) to speed up code in the OBBT algorithm and multi-start solver. Finally, network hydraulics were computed within CMS using the null space method proposed in Abraham and Stoianov (2015).
We evaluated the performance of the solution methods using three case study networks: Pescara, Modena, and BWFLnet. Problem data for each network, including the size of the resulting MINLP optimization problems, are summarized in Table 1. Moreover, network layouts are shown in Figure 5. Pescara and Modena are skeletonized versions of WDNs in medium-sized Italian cities, published by Bragalli et al. (2012) for benchmarking purposes. Since these are theoretical networks, unidirectional PRVs were placed at the outlet of each reservoir to represent existing valve configurations. BWFLnet is a large-scale operational network representing an area of the City of Bristol's WDN, located in southwest England. It currently operates with dynamically controlled PRV and bidirectional DBV valves, which modulate flow and pressure between adjacent zones to minimize AZP. BWFLnet's current control valve configuration was applied to model existing conditions -for more details on BWFLnet see Wright et al. (2014) and Waldron et al. (2020). Moreover, existing kept-shut boundary valves were opened in order to test the full capabilities of the proposed solution algorithm. Observe that the evaluated case study networks vary by size, network connectivity, and demand scenarios, enabling the proposed CMS algorithm to be tested on a range of network conditions. Pescara and Modena are relatively small-scale looped networks with theoretical (and quite large) demands. BWFLnet, on the other hand, represents a larger-scale operational network having a branched structure and demands representing water consumption in a typical urban area in the UK. For reference, we note that the size of the resulting MINLP optimization problems for BWFLnet are at the upper end of the range reported in the MINLPLib database (Vigerske, 2022); MINLPLib is a widely used library of benchmarking MINLP problem instances.  (14), we set the sigmoid function curvature parameter ρ to 50, unless stated otherwise. On the basis of the reviewed literature in Section 1, the self-cleaning velocity threshold u min j was set to 0.2 m/s for all j ∈ {1, . . . , n p }. Since we are concerned with the redistribution of peak flows, we define four discrete hydraulic time steps n t to represent peak demand scenarios within a typical diurnal demand profile. These were split between the morning and evening high demand periods and were unique to each case study. In relation to hydraulic feasibility constraints (11b), we applied a minimum regulatory pressure head of 15 m (UK regulations) at all nonzero demand nodes, for all case study networks. Moreover, the upper bound for flushing valve demands α U shown in constraint (11d) was set to 25 L/s. This value was selected with reference to the recommended fire flow capacity at fire hydrants in the UK (Local Government Association and Water UK, 2007). For the design-for-control problem, we evaluated the placement and operation of n v = 1, . . . , 3 bidirectional DBVs and n f = 0, . . . , 3 AFVs. Hence, for each network, a total of 12 valve configurations were considered. The number of randomization trials N in Algorithm 1 was set to the minimum of valve configuration combinations generated from the continuous relaxation in Subproblem (LP) or a maximum trial count; maximum trial counts were set to N = 50 for Pescara and Modena and N = 100 for BWFLnet. Finally, the number of initial conditions tested in the multi-start scheme (Figure 4) was set to M = 5.

Existing valve configuration results
We first investigated improvements to SCC performance achieved by optimizing the operational settings of existing pressure control valve configurations. This was considered a logical initial step as the implementation SCC-specific controls requires no additional capital investment and can thus be readily applied. As discussed previously, the theoretical case study networks, Pescara and Modena, were assumed to have PRVs assigned at each link leaving a source node (e.g. reservoir). On the other hand, BWFLnet had its existing operational valve configuration assigned. A subproblem of Problem (MINLP) was then formulated where binary variables z were fixed for the known control valve locations and directional variables v + and v − were set free for networks with existing DBV locations (e.g. BWFLnet). Valve control settings and DBV directions were then optimized over the four peak demand periods using the multi-start strategy described   Figure 6 through cumulative distribution plots of maximum flow velocities, weighted by pipe length. These plots offer a comparison between network SCC conditions with and without optimized controls. Considering f SCC is an approximation of the actual self-cleaning problem, we also include a sensitivity analysis on the sigmoid curvature parameter ρ to test changes in relative performance. Moreover, Table 2 summarizes the control problem results obtained with both the multi-start solver and the GA implementation described in Section 4.1. We set a maximum time limit of one hour for the off-the-shelf GA solver, which was assumed to reflect an upper bound for the practical application of control problems.
First, Figure 6 suggests that pressure modulation at existing control valves can significantly improve SCC conditions in the tested looped networks. This is highlighted by the numerical results for Modena, where the length of pipe experiencing self-cleaning velocities increased from approximately 60% to 95% for peak demand conditions. In contrast, results for BWFLnet, which has limited connectivity due to its  branched structure, showed only modest improvements. These results were expected as the reduced interconnectivity in BWFLnet leads to pipe flows being primarily a function of downstream demands. For this reason, Section 4.3 presents the results for a design-for-control problem where we consider the installation of additional DBVs and flushing valves, which aim to improve self-cleaning conditions in branched networks. Furthermore, we report modest changes in the flow velocity profiles for the ρ sensitivity analysis (see Figure 6). Since no pattern can be discerned between ρ values, we conclude ρ = 50 to be suitable for the design-for-control numerical experiments discussed in the subsequent section. As shown in Table 2, the multi-start solver computes solutions that are similar to those obtained by the GA, but at a fraction of the computational effort. While the GA consistently reaches the set timelimit of one hour, the longest computational time experienced by the multi-start solver is equal to 200 seconds. This is important as we are considering an optimal control problem, which may need to be solved in near real-time. The much longer CPU time required by the GA is likely due to the fact that it facilitates a search for the globally optimal solution. Thus, for highly nonlinear optimization problems, like that formulated for the current SCC problem, the GA retains a large solution space so as to avoid exclusion of the global optimum. This characteristic has been recognized as a challenge when using GAs for WDN control problems, especially when dealing with problems having a large number and range of free continuous variables (Maier et al., 2014, Ulusoy et al., 2022. Furthermore, we note that the GA implementation used in this work applies default parameters (see Section 4.1) and loads an external hydraulic solver (EPANET2.2) at each fitness evaluation. This off-the-shelf implementation is also likely to contribute to the relatively large CPU times. Extending the time limit to nearly 48 hours, the GA converges to a solution with an SCC value of approximately 24.8% for BWFLnet. This minimal performance gain suggests that the one-hour time limit is sufficient for the SCC control problem comparison.
In addition, Table 2 provides a comparison of average zone pressure (AZP) performance between control settings optimized for AZP and SCC, respectively. Note that the AZP objective function is formulated as a weighted sum of nodal pressures averaged over n t time steps (Wright et al., 2015;Equation 4). In accordance with findings in Abraham et al. (2016), we found a clear trade-off between the optimization of control settings for the SCC and AZP objectives. For example, a 16.2 m difference in AZP was recorded for BWFLnet. While this trade-off is not discussed further in the current study, we acknowledge the importance of having SCC and AZP objectives (among others) coexist in an overall operational network control strategy. Furthermore, the transition between such operational objectives may provoke large variations in hydraulic conditions (e.g. pressure transients) and should thus be managed accordingly. Future work is needed to investigate these considerations for the implementation of dynamically adaptive control strategies in operational WDNs.

Design-for-control results
This section considers the SCC design-for-control problem for optimal placement and control of new bidirectional DBVs and AFVs. We test a total of 12 valve configurations for each case study network, where Experiment 1 corresponds to n v = 1 and n f = 0 and Experiment 12 corresponds to n v = 3 and n f = 3. As stated in Section 4.1, each experiment uses hydraulic conditions from four peak demand conditions across a typical diurnal demand profile. Moreover, with the exception of existing DBVs in BWFLnet, the designfor-control numerical experiments include the existing PRV configurations shown in Figure 5, which were used for the optimal control subproblem in Section 4.2. We investigate the solution of Problem (MINLP) using three methods: (i) the proposed convex multi-start heuristic (CMS); (ii) CMS with domain reduction (as described in Section 3.1), referred to as CMS (DR); and (iii) an off-the-shelf GA implementation coupled with the EPANET2.2 hydraulic solver. Note that we set the following time limits for the GA on the basis of that recorded using CMS: 6 hours for Pescara and Modena; and 12 hours for BWFLnet. The SCC objective function results and corresponding CPU times from each solution method are shown in Figure 7 for the considered set of numerical experiments. Figures 7a, 7c, and 7e show that CMS finds feasible solutions for all numerical experiments. To highlight the improvements in SCC, we include a solid black line to represent initial hydraulic conditions and a dashed black line to report the SCC values obtained by optimally controlling existing valves, as done in Section 4.2. While SCC improvements are experienced for all case studies, we observe significant relative improvements compared with the control-only solution for BWFLnet. As previously discussed in Section 4.2, the SCC improvements for BWFLnet can be attributed to the deployment of AFVs, which facilitate controlled flushing demands at the optimally selected locations. Moreover, the additional DBVs act to further branch the network, thereby enhancing self-cleaning velocities. The increase in SCC is investigated further through network velocity plots shown in Figure 8. Here, we compare maximum pipe velocities resulting from initial hydraulics (i.e. no control) and, as an example, the design-for-control solution corresponding to n v = 2 and n f = 3. Figures 8b, 8d, and 8f clearly indicate a net increase in the overall extent of the network experiencing self-cleaning velocities from optimal valve controls. We observe optimal AFV placement to be near the network periphery in Pescara and BWFLnet. These results are intuitive as flushing demands facilitate self-cleaning velocities in areas that otherwise convey relatively low flow. In contrast, AFVs are concentrated near a source node in Modena. This is a result of existing PRVs modulating inlet pressure to shift the hydraulic balance point (or interface) between source nodes. In particular, the flow rates illustrated in Figure 8d suggest that source nodes near the top of the network supply a greater proportion of the demand. Additionally, we note that the optimal DBV locations, which increase self-cleaning velocities by redistributing flow paths, may result in sub-optimal conditions for controlling network pressure. Since we do not explicitly consider pressure management in the SCC problem formulation, the integration of selfcleaning and pressure management objectives should be a focus of future work. Lastly, Figure 8f highlights the limited influence of control features for the more branched BWFLnet network. Although the placement and operation of AFVs indeed improve self-cleaning hydraulic conditions, these are generally limited to direct upstream flow paths. In the case where self-cleaning velocities are sought for specific areas of the network (e.g. historical complaints), the SCC objective function can be modified to focus on a user-defined subset of pipes.
In comparing the two implementations of CMS, with and without domain reduction, we observe the former to yield better quality feasible solutions. This can be attributed to the tighter relaxations resulting from domain reduction. In fact, tighter relaxations may generate fractional valve placement values closer to the feasible domain, and therefore more likely to result in better quality design-for-control solutions. In addition, the solution of subproblem (LP) with tighter relaxations computes control valve settings η that are closer to the feasible set and the global optimum. As described in Section 3.3, these η values are used as a starting point for the multi-start solver. Hence, we expect these settings to be advantageous in finding a better solution. Nonetheless, the OBBT algorithm comes at an increase in computational effort. Figures 7b, 7d, and 7f compare CPU times across numerical experiments for each solution method. Although the increase in CPU times are modest for the smaller networks, OBBT does not scale well for the larger operational network. For example, the increase in CPU time reported for BWFLnet was, at minimum, approximately 50% when OBBT was applied. Even though Problem (MINLP) is considered a design problem, and is thus solved offline, OBBT can become impractical for large-scale and highly interconnected networks . Future work should investigate speed improvements to OBBT or the implementation of tighter relaxations for Problem (MINLP). Finally, we compare design-for-control solutions computed by CMS with an off-the-shelf GA implementation. As described in Section 4.1, the EPANET2.2 hydraulic solver was employed within the GA's fitness function to test hydraulic feasibility of the generated solutions. The solution was deemed feasible when the EPANET2.2 solver successfully converged and when the entire set of pressure heads across all time steps met the minimum pressure constraint. Otherwise, hydraulically infeasible solutions were discarded. Moreover, pressure set-points at control valves were bounded by the maximum head difference across the network (i.e. reservoir level less minimum node elevation). In Figure 7a, the GA solution is shown to produce comparable (if not better) SCC results for Pescara than CMS (DR). This is not surprising as the stochastic search procedure employed by the GA is well suited for highly nonlinear optimization problems, such as the current SCC objective function. Although the GA did not satisfy its termination criterion within the prescribed time limit (six hours), the relatively small number of free continuous variables in Pescara (e.g. pressure set-points and flushing rates; see Table 1) enabled the computation of good-quality feasible solutions. On the other hand, the GA was unable to find feasible solutions for numerous problem instances using Modena. This was particularly evident for experiments including the design and control of AFVs, which result in higher frictional losses materializing from flushing demands. Note that infeasibility is indicated by an SCC value of zero in Figure 7c. Additional GA experiments revealed such infeasibility issues to be linked to the relatively small range in hydraulically feasible pressures found in Modena. With the minimum regulatory pressure head relaxed from 15 m to 0 m (i.e. h min set to node elevations), the GA found feasible solutions comparable to CMS (DR). These results highlight the efficacy of the proposed convex heuristic for a wider range of network conditions. For BWFLnet, CMS (DR) finds similar quality design-for-control solutions to the GA. However, as observed in Modena, CMS (DR) is shown to consistently outperform the GA for experiments considering the implementation of AFVs. In light of the comparison between solution methods, we investigated the computational overhead of the GA fitness function, as configured in this study using the EPANET-MATLAB Toolkit. This is a cumbersome process, and one which directly impacts the number of fitness evaluations performed by the GA. By comparison, CMS requires much less overhead from loading external software and computes hydraulic states via an efficient null space solver. Table 3 compares the number of EPANET2.2 solver calls (for the 12-hour time limit) with the null space solver calls in CMS. The comparison uses select design-for-control experiments for the BWFLnet case study network. EPANET2.2 solver Null space solver 1 -n v = 1, n f = 0 11,064 4,208 3 -n v = 1, n f = 2 10,304 69,995 5 -n v = 2, n f = 0 11,064 61,274 7 -n v = 2, n f = 2 10,304 177,820 9 -n v = 3, n f = 0 10,304 87,976 11 -n v = 3, n f = 2 10,304 165,671 With the exception of Experiment 1, Table 3 reports the number of EPANET2.2 solver calls to be significantly less than the hydraulic simulations recorded in CMS. This identifies the hydraulic analysis as a bottleneck in the overall GA search procedure, particularly since its CPU times are over twice that reported using CMS. The comparison demonstrates the advantages of the proposed CMS algorithm in its ability to efficiently perform hydraulic simulations, which in turn yields more opportunities for CMS to compute solutions to the SCC problem than the GA over a given period.

Conclusions
Novel operational strategies are needed to minimize the severity and frequency of discolouration incidents in WDNs. In this paper, we propose control and design-for-control strategies for maximizing the self-cleaning capacity (SCC) of a network. This is facilitated through the operational framework of dynamically adaptive networks by extending its control capabilities to enhance water quality.
Building on previous work, we formulate an optimization problem to improve SCC through the optimal placement and operation of pressure control and automatic flushing valves. This results in a nonconvex MINLP problem. Since global solvers become intractable for large problem instances, we propose a tailored convex multi-start heuristic (referred to as CMS) to compute feasible solutions to the SCC design-for-control problem. CMS comprises a convex relaxation and randomized valve sampling procedure, which is demonstrated to effectively handle integer variables modelling valve placement. Given the high nonlinearity of the SCC objective function, CMS also implements a multi-start solver to compute feasible solutions to the SCC control problem, whilst mitigating the occurrence of poor local optima. Compared to an off-the-shelf GA implementation, which employs the EPANET2.2 hydraulic solver, CMS is shown to consistently yield good-quality feasible solutions to the SCC design-for-control problem. We demonstrate the robustness of CMS using networks with varying size and hydraulic complexity, including a large-scale operational network in the UK. Although the GA computed comparable solutions to CMS for most numerical experiments, it yielded worse-quality and, in some instances, infeasible solutions for larger and more complex conditions. Furthermore, we show the proposed multi-start solver to compute fast and scalable solutions to the SCC control problem, highlighting its suitability for application in near real-time control strategies. Altogether, the proposed convex multi-start algorithm represents progress towards the optimal design and control of smart water networks.
Moving forward, the proposed solution for maximizing SCC needs to be integrated within a wider framework for the design and operation of WDNs with dynamically adaptive control. For example, the objective to minimize average zone pressure and pressure variations should include the periodic transition to a selfcleaning mode, as formulated in this manuscript through the operation of pressure control and automatic flushing valves. From an operational perspective, it is also important that this control framework considers the management of unsteady hydraulics arising from such dynamic controls. In relation to the SCC objective function, additional features should also be investigated to improve self-cleaning performance. These may include the impacts of flow reversals as well as a tailored approach to focus on specific areas of the network (e.g. historical complaints). These considerations, combined with advances in technologies for the operation of dynamically adaptive networks, will enable better control strategies for responding to environmental, financial, and regulatory challenges.

A. Convex relaxations
This section presents a detailed derivation of the convex relaxations formulated for the nonconvex selfcleaning capacity (SCC) objective function and nonconvex HW head loss model.

A.1 Nonconvex SCC objective function
Recalling Section 3.1, we define piecewise linear relaxations for the following nonconvex sigmoidal inequality constraints where ψ + (·) and ψ − (·) model the positive and negative components of the SCC objective function f SCC , respectively; σ t := [σ + t σ − t ] T is the vector of auxiliary variables introduced to reformulate the SCC objective function; and q t ∈ R np is the vector of link flows. As presented in Section 3.1, the resulting relaxations are written as the following linear inequality constraint: and q U t and the sigmoid function parameters. Before proceeding with the derivation of (A.23), we note the following. First, we omit time indices t and link indices j for the sake of clarity. Second, we drop the sign of the sigmoid function ψ(·) and the vector of auxiliary variables σ ∈ R np ; apart from swapping the upper and lower flow bounds, the mathematical derivation does not change for the positive or negative ψ(·) functions. However, where deemed necessary to clearly explain the formulated relaxations, we use the positive sigmoid function. Finally, we use flow velocities as the domain for ψ(·), which are denoted by u := q A . It follows that the lower and upper domain bounds are denoted by u L := q L A and u U := q U A , respectively. We begin by considering the set (u, σ) σ = ψ(u), ∀u ∈ u L u U , (A.24) and its concave envelope given by whereψ(·) is a piecewise linear relaxation that depends on bounds u L and u U and the tangent point created with the sigmoid function. Following the examples presented in Boyd (2014, 2016), we discuss the steps to constructψ(·) for various domain scenarios below. First, we look for the point w such that the line from (u L , ψ(u L )) to (w, ψ(w)) is tangent to ψ(·) at (w, ψ(w)). The bisection method, an iterative numerical algorithm (Miller, 2014), is used to compute the tangent point w as no analytical expression exists for the equations of intersection in consideration. This follows the steps detailed in Algorithm A2. Since the SCC problem is mostly comprised of linear constraints, Algorithm A2 Bisection method for nonconvex sigmoid objective function 1: Input: ψ(·) model parameters, SCC threshold velocity u min , and flow velocity bounds u L and u U 2: Output: unique point w where the line from (y, ψ(y)) to (w, ψ(w)) is tangent to ψ(·) at (w, ψ(w)) 3: Define: for ψ + , y := u L ; for ψ − , y := u U ; equation of intersection f (x) := ψ (x)(x−y)+ψ(y)−ψ(x), where x = w at the root of f (x) 4: Initialize: (i) initial lower and upper bounds x 0 and x 1 ; and (ii) error tolerance tol 5: if u L ≥ u min then 6: w ← u L 7: else 8: x 1 ← x 2 12: x 2 ← x0+x1 w ← x 2 19: end if it is also convenient to construct a linear upper bound to the concave segment of ψ(·). It follows that point k represents the intersection of linear relaxations, at which point the maximum error ofψ(·) − ψ(·) occurs. The constructed concave envelopesψ(·) are described below for the four conceivable domain scenarios, with their corresponding plots shown in Figure A.9.
If u L < w ≤ u U , the concave envelopeψ(u) is defined as the following piecewise linear function (see Figure A.9a):ψ When u L < u U ≤ w, the concave envelopeψ(u) is defined a the following linear function (see Figure A.9b): When u L = w and w ≤ u U , the concave envelopeψ(u) is defined as the following piecewise linear function (see Figure A.9c):ψ Lastly, if u L = w = u U , the concave envelopeψ(u) is defined as the point (see Figure A.9d): Sigmoid functions ψ + and ψ − in (A.22) are then replaced with the appropriate concave envelopesψ(·) described in (A.26) -(A.29). After rearranging terms, and substituting flow velocity u with q A to align with the SCC problem variables, the linear inequality constraint shown in (A.23) is formulated.

A2. Nonconvex HW head loss constraints
Building on the material presented in Pecci et al. (2019) (Appendix 1), we formulate polyhedral relaxations for the HW head loss model. Recall the energy conservation equality constraint where φ(·) models the head loss materialized across network links j ∈ {1, . . . , n p }; θ t ∈ R np is the vector of auxiliary variables introduced to isolate the nonconvex head loss term; and q t ∈ R np is the vector of link flows. The resulting polyhedral relaxations are written as the following linear inequality constraint: where matrices R t and E t and vector r t depend on flow bounds q L t and q U t and the HW model parameters. We follow a similar procedure to that detailed in Liberti and Pantelides (2003) for monomials of odd degree and later adapted in Pecci et al. (2019) for the quadratic approximation head loss model. The contribution of this work concerns the derivation of polyhedral relaxations for the HW head loss model, which has a fractional exponent dictating the relationship between q and φ(·).
Omitting indices t and j, we begin by considering the set and its convex relaxation given by (q, θ) φ(q) ≤ θ ≤φ(q), ∀q ∈ q L q U , (A.33) where φ(·) andφ(·) are piecewise linear functions which formulate lower and upper relaxations of φ(·), respectively. These functions depend on flow bounds q L and q U and the tangent points created with φ(·). In the following, we present the derivations of φ(·) andφ(·) for the complete range of domain scenarios. We start by looking for the point z such that the line from (q L , φ(q L )) to (z, φ(z)) is tangent to φ(·) at (z, φ(z)). In contrast to the quadratic approximation method derived in Pecci et al. (2019), there is no analytical expression for finding root z to the intersection between φ(·) and the tangent line at (z, φ(z)). Thus, the bisection method is used to find z within an appropriate numerical tolerance. While similar to Algorithm A2, we present the exact steps taken to compute z in Algorithm A3 for the sake of completeness. Analogously, the unique pointz is computed to find the line from (q U , φ(q U )) to (z, φ(z)). The resulting set of linear relaxations for φ(·) are described below, with the corresponding plots presented in Figure A.10. Note that we define the set of polyhedral relaxations asφ(·) := {φ ∪φ}. Moreover, let k andk denote the points of intersection between the linear functions in φ(·) andφ(·), respectively.

B. Optimization-based bound tightening algorithm
In Section 3.1, we implement an optimization-based bound tightening (OBBT) algorithm to strengthen the convex relaxations formulated for Subproblem (LP). Following the supplementary material from Pecci et al. (2022), we tighten polyhedral relaxations for the nonconvex constraints by reducing the flow variable domain for link indices corresponding to the set of core links C. For each OBBT iteration, we solve two optimization problems to independently maximize and minimize flow q j,t for all links j ∈ C and hydraulic time steps t ∈ {1, . . . , n t }. These are represented by the linear program described in Subproblem (LP), with objective functions set to both minimize and maximize q j,t . We therefore solve 2n t |C| linear programs for each OBBT iteration. The algorithm is terminated once bound tightening progress diminishes between successive iterations, which is defined by a set tolerance on the maximum flow domain diameter. Moreover, we exclude forest links from the optimization scheme since their flows are defined a priori by the aggregate of downstream demands and the maximum flow defined for automatic flushing valves. Psuedocode for the OBBT algorithm is described in Algorithm B1.

17:
k ← k + 1 18: end while C. Strictly feasible sequential convex programming algorithm In Section 3.3.1, we implement a strictly feasible sequential convex programming (SFSCP) solver to compute feasible solutions to the nonlinear programming (NLP) control problem. This follows the implementations in Wright et al. (2015) and Abraham et al. (2016), with the addition of the continuous decision variable α to model the operation of flushing valves. We present the SFSCP solver psuedocode in Algorithm C1, which is adapted from Wright et al. (2015). The variables used in Algorithm C1 are defined as follows: f denotes the continuous sum of sigmoids SCC objective function f SCC ; η and α model operational settings for pressure control and flushing valves, respectively; dη and dα denote step directions for the aforementioned continuous variables; q and h represent hydraulic states computed for each control configuration; and F is the set defining the hydraulically feasible region.
Algorithm C1 SFSCP solver 1: Input: network data, fixed valve configuration, and feasible starting point x 0 2: Output: f * and corresponding operational control settings η * and α * 3: Initialize tol , k max , and set k = 1 4: Set x k ← x 0 5: Compute initial hydraulic states q k , h k , and objective function f k 6: while |f k −f k+1 | |f k | > tol and k ≤ k max do 7: η k+1 , α k+1 ← NLP control problem with first-order Taylor approximations of nonlinear SCC objective function and energy conservation constraints 8: dη k , dα k ← η k+1 − η k , α k+1 − α k 9: β = 1 10: q k+1 , h k+1 ← hydraulic simulation with controls η k + βdη k , α k + βdα k 11: f k+1 ← recompute objective function 12: while f k > f k+1 or x k+1 / ∈ F do 13: β ← β/2 14: η k+1 , α k+1 ← η k + βdη k , α k + βdα k 15: q k+1 , h k+1 ← hydraulic simulation with controls η k + βdη k , α k + βdα k f k ← f k+1 21: end while D. NLP solver performance comparison We use the SFSCP solver in this work for multiple reasons. First, it guarantees strict hydraulic feasibility for each optimization step and shows fast convergence properties in previously studied WDN control problems, as demonstrated in Wright et al. (2015). Second, since state-of-the-art nonlinear optimization solvers (e.g. IPOPT) require second-order derivatives, it avoids the complications associated with the nonsmooth HW head loss model, which has an unbounded hessian at the origin. Potential complications are exacerbated by the highly nonlinear SCC objective function used in this study. In order to justify application of the SFSCP solver, though, we conduct a performance comparison with solutions computed by IPOPT. Following Dolan and Moré (2002), we derive performance profiles to objectively compare the two solvers. Let X := {1, . . . , n x } and S := s SFSCP , s IPOPT be the set of test experiments and solvers, respectively, where n x is the number of experiments tested. Since we are interested in comparing the SCC objective function values, we set the performance metric f x,s equal to f SCC for each experiment x produced by solver s. We then compare solver performance for all experiments x ∈ X by normalizing f x,s to the best computed solution for a particular experiment by any solver s ∈ S. This comparison is expressed as follows: where r x,s is the performance ratio of solver s for experiment x. Note that, if the solver produces an error or fails to find a feasible solution to the SCC control problem, we assign f x,s = +∞. An overall performance assessment of solver s ∈ S is then obtained by deriving the cumulative distribution function where ρ s (τ ) represents the percentage of experiments within a factor τ of the best solution for solver s ∈ S. Figure D.1 presents the SFSCP and IPOPT solver performance profiles for n x = 36 design-for-control experiments. In addition to existing unidirectional pressure reducing valves (PRVs), these experiments involved a combination of n v = 1, . . . , 3 bidirectional dynamic boundary valves (DBVs) and n f = 0, . . . , 3 automatic flushing valves (AFVs), as described in Section 4.1 of the manuscript. From Figure D.1, we first observe that the SFSCP solver finds feasible solutions for all experiments, while IPOPT returns infeasibility errors for just over 20% of the tested cases. This is illustrated by the distance along the x-axis from the origin to the vertical segment of the IPOPT solver (blue) plot. When looking at individual experiments, we identified a few features which may be contributing to IPOPT's infeasibility errors. First, infeasibility occurs almost entirely for the larger case study network (e.g. BWFLnet), of which also has a greater range in operational demands and control features. Second, since good performance was recorded for the feasibility restoration problem (see Section 3.3.2 of the manuscript), we believe the highly nonlinear SCC objective function, in combination with the unbounded hessian of the HW model, to be a key factor in the encountered infeasibility errors. In any case, Figure D.1 shows the SFSCP solver to be within 10% of the best solution for all tested experiments, which is comparable to the IPOPT solver for experiments when no feasibility issues were present. This suggests that the SFSCP solver is appropriate for the current SCC design-for-problem study.