S Pellerin
email: s.pellerin@limsi.fr

B Podvin

Y Fraigneau

G Bonnavion

O Cadot

Near wake of an Ahmed body: analysis of switches and base pressure-velocity relationship

Keywords: Aerodynamics, Ahmed body, direct numerical simulation, experiments, POD

We perform the flow behind a squareback Ahmed body using DNS. The dynamics of the flow is analyzed by Proper Orthogonal Decomposition to extract its main characteristics. We focus on the relationship between the base pressure and the near-wake velocity field. The dynamics of the switches due to bimodal behavior of the flow is also considered. The linear stochastic estimation technique [4] is applied to estimate the base pressure modes from the wake flow field and vive-versa. We show that the velocity-based estimator can recover dominant base pressure fluctuations and that the pressure-based estimator can detect velocity deviation and switch mode but fails to capture shedding modes.

Introduction

In the context of drag reduction associated with environmental issues, we focus on the flow behind an Ahmed body with a square base, mimicking an SUV type profile. It has been observed experimentally [START_REF] Grandemange | Turbulent wake past a three-dimensional blunt body. part 1 -global modes and bi-stability[END_REF] that the near wake of this flow is characterized by a bimodal behavior with two symmetrical quasistationary states and switches on either side in the spanwise direction. This phenomenon is important to control because it leads to an increase of the drag associated with the base pressure evolution [START_REF] Bonnavion | Unstable wake dynamics of rectangular flat-backed bluff bodies with inclination and ground proximity[END_REF][START_REF] Haffner | Mechanics of bluff body drag reduction during transient near-wake reversals[END_REF]. The flow is calculated by DNS at Reynolds number Re = 40000, based on the incoming fluid velocity and the body height. We first investigate the dynamics of the flow due to these switches. The main topic of this study is to estimate the velocity field from the base pressure data and, conversely, to estimate the base pressure modes from the velocity field.

Summary of results

The Direct Numerical Simulation code SUNFLUIDH has been used to calculate the flow behind the squareback Ahmed body at Re= 40000. The Proper Orthogonal Decomposition is applied on numerical data which permits fist an fine analysis of the dynamics of the flow from the main energetical modes and highlights deviation, switches and shedding modes characteristics. In previous studies, POD gave a good comparison between numerical and experimental data in the near-wake [START_REF] Podvin | Proper Orthogonal Decomposition analysis and modelling of the wake deviation behind a squareback Ahmed body[END_REF]. The behavior of the switches between the two quasi-stationary states has been also studied and reduced models have been established to reproduce the switches over long periods of time [START_REF] Podvin | Low-order modelling of the wake dynamics of an Ahmed body[END_REF].

The figure 1 shows the evolution of deviation and switch modes of the velocity field, over the calculated period with DNS; the switch period is colored. The switch velocity mode (right) is positive during the switch, evidenced on the deviation mode evolution (left). We investigate here what is the footprint of the wake flow on the base pressure. Note that there are two very different kind of data: the pressure field is characterized by one component and two dimensions (1C-2D) while the velocity field is 3C-3D. We apply the POD approach to the pressure field DNS data. We first extract the POD dominant modes of the base pressure (figure 2 left) and obtain a good comparison with experimental data, on spatial modes and temporal evolution of the amplitudes. As we have seen studying the velocity field, we find analysing the pressure field that the base pressure also increases during the switches. This increase can be related to edge deformation of separation region. then apply the linear stochastic estimation technique of Podvin et al. [START_REF] Podvin | On combining linear stochastic estimation and proper orthogonal decomposition for flow reconstruction[END_REF] to the both velocity and base pressure fields. Each corresponding transfer function is then calculated. We can then estimate the pressure amplitudes from a velocity-based estimator. The comparison between the true field and the estimation is very correct and recover the dominant base pressure fluctuations (figure 2 right). We then apply the opposite approach to estimate the velocity field amplitudes from a pressure-based estimator. This estimation detects the velocity deviation and the switch mode. However, the shedding modes are not captured by this approach.

Figure 1 :

 1 Figure 1: Amplitudes evolution of the velocity modes, POD from DNS data: (left) deviation mode and (right) switch mode.

Figure 2 :

 2 Figure 2: (Left) Mode 2 of the pressure field on the base, POD from DNS data and (right) amplitudes evolution of the pressure mode 2, (top) true field and (bottom) velocity-based estimation field.


Sylvan Elhay
email: sylvan.elhay@adelaide.edu.au.

Olivier Piller
email: olivier.piller@inrae.fr

Jochen W Deuerlein
email: deuerlein@3sconsult.de

Angus R Simpson
email: angus.simpson@adelaide.edu.au

M Asce

An interior point method applied to flow constraints in a pressure dependent water distribution system

Keywords: Networks, convex programming, interior point methods, active set methods, flow control

In this paper, the fusion of an active set method (ASM) and an interior point method (IPM) is used, with a Newton method, to solve for the steady-state flows, heads and outflows of a pressure dependent water distribution system. The outflow constraints which arise from the pressure dependency are handled by an ASM and the linkflow constraints are handled by an IPM. The authors believe that this is the first time an ASM and an IPM have been used together in this way to solve a real world optimization problem. Including flow constraints in a network model allows a variety of flow control devices (flow control valves, check valves, pumps) to be modelled e ciently. The new method does not require damping. The separate treatment methods for the two constraint sets means that the linear inequality constraint qualification condition cannot be violated during iteration, unlike the case where all the constraints are handled by an ASM. The method is shown to quickly converge on nine case study networks the largest of which has more than 157,000 links and 150,000 nodes and 6,000 linkflow constraints. When tested on those same nine networks, the popular interior point optimizer package IPOPT is shown to take between about 2-4 times as long as the new method for 5 of the networks, 7 and 9 times as long for 2 networks and about equal time on the 2 remaining networks. For the largest network IPOPT takes 34 minutes while PHYAI takes 4 minutes. The generality of the method means that it has application to network design and management, capacity analysis and self-cleaning networks.

INTRODUCTION Water Distribution Network Modelling Background

The traditional demand dependent models (DDMs) ( [START_REF] Cross | Analysis of flow in networks of conduits or conductors[END_REF], [START_REF] Martin | The application of Newton's method to network analysis by digital computer[END_REF], (Wood & Charles 1972), (Todini & Pilati 1988)) of water distribution systems (WDSs) deliver the nominal nodal demands regardless of the pressures in the system.

This can lead to model solutions which are mathematically correct but which are not physically realizable. For example, if the pressure at a node falls below a certain level in a real network, then the nominal demand cannot be realized. As a consequence, attention has turned to pressure dependent models (PDMs) of WDSs and much research e↵ort has been directed at finding powerful and robust solvers to quickly find the steady-state linkflows, outflows and heads of real PDM WDSs.

PDM models require a pressure outflow relationship (POR ) which describes the outflow at each demand node as a function of available pressure: nodes at which the pressure is below a preset minimum level, h m , have zero outflow while those at which the available pressure is above a maximum preset service pressure, h s , will have the nominal demand delivered. The POR describes the proportion of the nominal demand that is delivered for pressures between h m and h s . The interested reader is directed to [START_REF] Sayyed | Noniterative application of epanet for pressure dependent modelling of water distribution systems[END_REF] or [START_REF] Deuerlein | A content-based active set method for the pressure dependent model of water distribution systems[END_REF] for discussions of some commonly used PORs .

Both [START_REF] Bhave | Node flow analysis distribution systems[END_REF] and Tabesh (1998) proposed solving the PDM problem by using a two-step iterative procedure called Node Flow Iteration. Here the problem is first solved as a DDM and then the demands are corrected according to a chosen PDM relationship.

More recently, [START_REF] Jun | Iterative methodology of pressure-dependent demand based on EPANET for pressure-deficient water distribution analysis[END_REF] proposed a solution technique in which the PDM problem is solved by repeatedly solving the corresponding DDM problem by the Global Gradient Algorithm (GGA) of Todini & Pilati (1988) in an extension to the widely-used WDS DDM solver, EPANET [START_REF] Rossman | EPANET 2 Users Manual, Water Supply and Water Resources Division[END_REF] and adjusting the demands after each solution. [START_REF] Lippai | Demand constructs for risk analysis[END_REF] introduced artificial reservoirs, artificial check valves and artificial flow control valves to prevent the outflows from lying outside a prescribed interval. Elsewhere, Piller & van Zyl (2014a) used the power equation or the Fixed and Variable Area Discharge (FAVAD) pressure-dependent leakage equation at nodes with leakage to model the dependence of flow on pressure. A review of many of the available approaches for simulating or emulating PORs can be found in Suribabu et al. (2019).

The PDM problem can alternatively be viewed as an optimization problem with linear constraints on the outflows at the nodes and the authors used this approach to solve the unconstrained primal-dual and dual problems by a Newton method with line search to find the heads, linkflows and outflows of a PDM WDS [START_REF] Elhay | A robust, rapidly convergent method that solves the water distribution equations for pressure dependent models[END_REF]. The authors later developed [START_REF] Deuerlein | A content-based active set method for the pressure dependent model of water distribution systems[END_REF]) a content-based active set method (ASM) for the PDM convex optimization problem min q,c f (q, c)

| A T q + c = o, g i (c)  0, i = 1, 2, . . . , m g , (1)
where the content function, f , is strictly convex, A is a node-arc incidence matrix, the g i define m g linear constraints and q and c are the vectors of decision variables representing the network link flows and nodal outflows, respectively. The resulting ASM is robust and rapid and well-suited for PDM WDSs. Indeed, the Newton-type methods used in the ASM are particularly suitable for WDS models because even though the equations involved are non-linear, they are very structured and sparse and so even models that have hundreds of thousands of degrees of freedom can be quickly and conveniently solved on modern desktop computers. In addition, using the inverse POR , h(c), as the authors did, brings a computational advantage: the Darcy-Weisbach head loss model is quadratic and pseudohead loss function, h(c), for the Wagner (and any similar POR ) is also superlinear. The consequence is that, although in similar circumstances damping or regularization might be necessary, in this case they are not because the convergence ball for the problem is larger and as a result good convergence is more easily achieved.

Addition of Controls

However, when combinations of throttle control valves (TCVs), check valves (CHVs), flow control valves (FCVs), pressure sustaining valves (PSVs), pressure reducing valves (PRVs) or pumps are added to the models, new di culties arise [START_REF] Deuerlein | Preprocessing of water distribution systems to assess connectivity and solvability in the presence of flow control devices[END_REF]. In some cases the methods fail to converge (see [START_REF] Simpson | Modeling of pressure regulating devices a major problem yet to be satisfactorily solved in hydraulic simulation[END_REF] and [START_REF] Deuerlein | The never-ending story of modeling control devices in hydraulic systems analysis[END_REF]) or even converge to the wrong solution [START_REF] Gorev | Non-unique steady states in water distribution networks with flow control valves[END_REF]. These di culties are compounded by the fact that in many cases there is no indication that the computed solution is in fact wrong or not physically realizable.

Control devices are increasingly used in WDSs. Limited resources during hot summers, contamination of ground water and the need for more flexible and resilient systems are only a few of the drivers for this development. There is no doubt that the planning and operation of modern systems using smart devices require robust and e cient hydraulic modelling software.

An FCV constrains the flow in a link to be below some chosen value, q max , that is independent of the di↵erence between the heads at the link's initial and final nodes. If the flow is below q max , the valve is opened fully and behaves like a minor loss element.

If conditions are such that the unrestricted flow would be higher than q max , the FCV acts to restrict the flow. It is natural then to model the action of such a device by a mathematical constraint. Applying linkflow constraints to the minimization of the content associated with a WDS is then a constrained optimization problem and this formulation of the problem provides some interesting returns: the characterisation of the conditions under which solutions exist and are unique can be studied and a wide range of optimization techniques, including ASMs and interior point methods (IPMs), can be brought to bear on the problem.

Formulation of the Problem Addressed

A number of di↵erent approaches to the solution of flow constrained problems have been reported. [START_REF] Piller | Modeling control valves in water distribution systems using a continuous state formulation[END_REF] used quadratic external penalty methods applied to the valves' head loss equations. The penalty terms are interpreted as a minor head loss. Later, [START_REF] Alvarruiz | E cient modeling of active control valves in water distribution systems using the loop method[END_REF] and [START_REF] Rossman | [END_REF] used linear external penalty methods in the loop method to model flow regulating devices. The external penalty terms in these methods generate very large resistance factors and, when these exist alongside normal resistance factors in the models, the resulting matrices can have very large condition numbers. This leads, in many cases, to the need for damping or regularization schemes to be incorporated into the methods in order to avoid divergence.

In [START_REF] Piller | A content-based active set method for the pressure dependent model of water distribution systems with flow controls[END_REF] the authors presented a new content-based active set method, ASMFC, for the determination of linkflows, nodal heads and outflows in a PDM WDS which has linkflow constraints. This time the problem addressed was min q,c f (q, c) |

> <

> :

A T q + c = o, g i (c)  0, i = 1, 2, . . . , m g , s j (q)  0, j = 1, 2, . . . , m s . ( 2
)
The necessary and su cient Karush-Kuhn-Tucker (KKT) conditions for the unique minimum of the content function, f , optimization problem were given and a Newton method based on this formulation led to the ASMFC. The method was demonstrated by applying it to eight case study networks, N 1 -N 8 , with between 934 and 19,647 links and between 848 and 17,971 nodes. These case study networks each had 60 cotree linkflow constraints which either (i) limit the maximum flow in a link to be between 0 and ±10% of its unconstrained value or (ii) prescribe a link flow equality constraint (LFEC), i.e. the flow in link j is set by requiring q max,j = q min,j . The method rapidly found the solutions to all eight case study networks. An advantage of the method is that FCVs and pumps can be modelled and heuristics are not needed to determine the status of control devices in the system: their states are found as part of the solution. Furthermore, there is no risk of isolated demands in this case since in PDM problems the nodal outflow can reduce to zero if the pressure is insu cient.

In some settings flow velocities are required to be above a certain lower limit. Selfcleaning networks [START_REF] Abraham | Decreasing the discoloration risk of drinking water distribution systems through optimized topological changes and optimal flow velocity control[END_REF], for example, limit the flow capacity in all the network's cotree (by closing some valves) to ensure higher velocities in the spanning tree links. Interest in this area led the authors to conduct investigations into the performance of the ASMFC on problems like these. In these new tests the eight case study networks had constraints applied to (almost) every link in each network's cotree (rather than limiting the number of constraints to 60). Confining the constraints to the network cotrees ensured that no parts of the networks were isolated from sources. This led, for the largest network, to ⇡ 1700 constraints being imposed. It was found that the ASMFC converged rapidly for six of the eight networks but for networks N 4 and N 7 , convergence with the ASMFC remained elusive. Now, it is known (see [START_REF] Piller | A content-based active set method for the pressure dependent model of water distribution systems with flow controls[END_REF] for details) that, provided certain conditions which are known to hold in this case are met, the solution to this problem always exists and is unique. The problem to finding these solutions, it emerged, was poor starting values -a common problem in the numerical solution of nonlinear systems. If, on these more challenging problems, the ASMFC is started at a point close to the true solution it reliably converges quite rapidly. But without good starting values stalling and cycling (failure to converge) were observed.

Damping, a technique used by the authors in [START_REF] Elhay | A robust, rapidly convergent method that solves the water distribution equations for pressure dependent models[END_REF] to overcome convergence failures in PDM problems, has been helpful where initial values are far from the solution or where link characteristics are modelled by sublinear functions or even when very large and very small resistance factors in the same models lead to poor matrix conditioning. It has, however, the disadvantage that the co-energy functions (content, co-content or Lagrangian) need to be frequently evaluated and these add an undesirable computational burden. This led the authors to investigate alternative methods and resulted in the new hybrid method of this paper, HYAI.

The Solution

The di↵erence between the problems defined in ( 1) and ( 2) turns, of course, on the addition of the linear inequalities s j (q)  0 (the linkflow constraint set) and seeing that those constraints are not always well-handled by the ASM, an alternative way of handling them was found. The constraints s j (q)  0 in (2) could theoretically be implemented by omitting those constraints from the definition of the problem and adding an indicator function term I (s j (q)) to the objective function. Here, I (•) is defined by

I (x) def = ⇢ 0 x  0 1 x > 0 . ( 3
)
Log barrier IPMs replace the function I (•) in (3) by a logarithmic penalty function, 1 t (q), which drives the optimization away from the constraint boundaries while keeping them in the interiors of the constraint intervals. The key element in the HYAI method, and the reason it is said to be hybrid, is that it uses a log barrier penalty to enforce the linkflow constraints while continuing to use the ASM approach to handling the POR constraints (via the inverse POR function) on the nodes. The new method is denoted by HYAI to suggest the hybridization of an ASM and an interior point method (IPM). The method borrows from the theory of IPMs for dealing with the linkflow constraints but cannot properly be said to be an IPM in view of its use of the ASM. It will be seen that implementing HYAI requires only a small modification of the ASMFC method.

The penalty terms in HYAI are such that the barrier terms more closely approach I (•) as t increases. Thus, larger values of t give solutions which are closer to the exact solution or the ASMFC solution. But finding the minimum of f (x) can be di cult with Newton's method when t is large because the Hessian of f varies rapidly near the boundaries of the feasible set (Boyd & Vandenberghe 2009, 564). This di culty can be overcome by first finding a solution for a small value of t and then solving a sequence of problems in which the solution for each value of t is used to start the Newton method for a larger value of t.

This process of continuation can then be used to arrive at a solution for su ciently large t.

When a linkflow constraint in an ASM saturates, that e↵ectively causes the link to be removed from the system. Removing a link from the system can, in some instances, cause disconnections in the system's main graph. By comparison, in the hybrid IPM of this paper the full arc-node incidence matrix (ANIM) is always used and some terms are added to the diagonal matrix of head losses. These added terms can assist with convergence.

Unlike external penalty methods which sometimes need damping, in this log barrier IPM, damping is avoided by projecting linkflows back to the interior of the constraint interval.

The uniqueness and existence of the problem flows is assured provided the linearly independent constraint qualification (LICQ) is not violated at the solution [START_REF] Deuerlein | On the solvability of the pressure driven hydraulic steady-state equations considering feedbackcontrol devices[END_REF]) because if the LICQ holds then unique Lagrange multipliers exist and a feasible solution exists. The LICQ requires that the matrix of equality constraints (the mass balance constraints together with the saturated or active constraints) has full rank.

For the ASMFC method, the linkflow constraints are part of the LICQ matrix. The violation of the LICQ is equivalent to the isolation of some demand nodes and this can occur if, for example, the linkflow constraints are not confined to the cotree. So, where there are constraints on links in the spanning tree, one should establish (although it is usually impractical to do so) at each iteration that the LICQ is not violated.

The LICQ matrix for the HYAI method is comprised only of (i) the network node-arc incidence matrix, and (ii) the matrix of active outflow constraints, since the linkflow constraints are handled by the IPM, not the ASM. In the Appendix it is shown that the HYAI method LICQ matrix has full rank provided the NAIM has full rank and importantly, that the LICQ matrix and the ANIM of the augmented graph of the network, a key matrix in the method, have the same rank. This fact eliminates a point of failure in the Newton method which is at the basis of the new method. Whether or not a solution exists in the feasible set can be determined by a linear program (see Boyd & Vandenberghe (2009, 579) for details).

Contribution of the New Method

Network designers seek to balance conflicting requirements: (i) the need to satisfy the demands and (ii) the cost of the network construction and maintenance. The search for this balance leads to the need to minimise the head losses, which depend on link diameters and the squares of the velocities (thereby limiting the e↵ects of abrasion at singularities, elbows and shrinkages). A PDM solver with linkflow capacity constraints can be used to help find the optimal solution which balances these requirements. An important problem concerns modelling network capacity: the study of networks for which every link has flow constraints. Capacity problems arise, for example, in transport (communications networks and urban tra c flows) and power systems. A WDS capacity analysis which finds flows that are such that the corresponding velocities are bounded above and below, v min  v  v max can assist network designers in two important ways. The first is that the lower fluid velocity limit can be set su ciently high to ensure that the links so constrained are self-cleaning i.e. have flow velocities high enough to prevent the particle sedimentation that would otherwise occur. For example, in some settings, the fluid velocities in a WDS are required to be limited to no more than 1.5 m/s (sometimes the linkflow constraints are limited to every cotree link and are set to create higher velocities in just the spanning tree links).

The second way that capacity analysis assists the designers stems from the fact that (Boyd & Vandenberghe 2009, 252) the Lagrange multipliers measure the first order sensitivity of the objective function's (in this case the content's) optimal value with respect to perturbations of the constraint. Thus, an examination of the steady-state Lagrange multipliers for the lower velocity bound could indicate how much pumping is necessary in a system and how much extra energy is required to achieve that level of pumping or reduction of head loss.

The authors believe that the HYAI method is the first in which the constraints for some of the problem's objective function variables are handled by using an active set method while others are handled by an interior point method in this way. Its suitability to this approach stems from the fact that the optimization variables naturally separate into two sets: those for the nodes and those for the links.

In this paper the HYAI method is developed and its e↵ectiveness is evaluated and compared to that of ASMFC by its application to the eight case study networks previously used in the evaluation of the ASMFC together with a ninth, much larger network. The new method is shown to solve all the networks considered in this study (some of which previously could not be solved with ASMFC) and its reliability is established by showing the high level of agreement between its solutions and those of ASMFC.

The method is tested on nine networks the largest of which has about 6,000 linkflow constraints: a number which is unrealistic for network control but which is entirely reasonable in the context of design and post-analysis of the resulting outflows and saturation capacities. Moreover, the authors believe that the intensive study of the mathematical background (i) is of value in its own right and (ii) helps provide robust solutions to practical problems. Existing solvers often struggle to find the correct solution even in rather simple cases. It is the objective of the authors to provide a mathematically-based method that is able to find the correct solution for all cases where such a solution exists.

It is shown that, when tested on the nine case study networks, the popular interior point package IPOPT (Wachter & Biegler 2006) takes about 2-4 times as much wall-clock time as PHYAI on 5 networks, 7 and 9 times as long on 2 networks and about the same time on the remaining two networks. In the case of the largest network, IPOPT took about 34 minutes to solve the problem while the new method took about 4 minutes.

The solution of a small example problem is described in detail in the Appendix to help clarify the application of the new method.

The usefulness of the new method, either in its own right or as a method that can be paired with ASMFC, is discussed. Some summary statistics for four of the nine case study networks which are in the public domain and summary statistics for the corresponding solution data are provided as supplementary material to the paper.

DEFINITIONS AND NOTATION

Consider a WDS whose network graph has n p links, or arcs, and n j + n f nodes, or vertices: n j is the number of nodes at which the heads are unknown and n f 1 is the number of source nodes with fixed heads. The links of the network include control valves, pumps and pipes. Denote by q = (q 1 , q 2 , . . . , q np ) T 2 R np the vector of unknown flows in the system, h = (h 1 , h 2 , . . . , h n j ) T 2 R n j the unknown heads at the nodes in the system, u = (u 1 , u 2 , . . . , u n j ) T 2 R n j the vector of node elevations and r(q) = (r 1 , r 2 , . . . , r np ) T the vector of link resistance factors. Let A denote the n p ⇥ n j , full rank, unknown-head ANIM: the ji element of A is (i) 1 if node i is at the end of arc j, (ii) 0 if arc j does not connect to the node i, and (iii) 1 if arc j starts at node i. Let A 0 denote the ANIM, with a similar definition, for the fixed-head nodes. Let h 0 denote the vector of elevations of the n f fixed-head nodes. Denote a = A 0 h 0 . Denote by ⌘ the exponent used in the head loss formula: ⌘ = 2 for the Darcy-Weisbach model and ⌘ = 1.852 for the Hazen-Williams model. Furthermore, denote by G(q) 2 R np⇥np the diagonal matrix whose diagonal

elements are defined as [G(q)] jj = r j |q j | ⌘ 1 . Then, ⇠(q) = (⇠ 1 (q 1 ), ⇠ 2 (q 2 ), . . . , ⇠ np (q np )) T =
G(q)q is the vector whose elements model the head losses of the links in the system. In general, (e.g. for the Darcy-Weisbach formula) r = r(q) but for the Hazen-Williams formula r is independent of q. Denote the vector of the nominal demands at the nodes what follows, the symbol O denotes a zero matrix and o denotes a zero column vector of appropriate dimension for the particular case. The shorthand notation x + a, where x is a vector and a is a scalar, will be used to denote the case where every component of x has a added to it. Furthermore, it will be assumed that any matrix inverses which are shown do, in fact, exist. In what follows all the values of t will be powers of 10: t = 10 e .

Turning now to PDM problems, assume, for simplicity and without loss of generality, that every node has the same minimum pressure head, h m , and the same service pressure head, h s . Individualized minimum and service pressure heads can be implemented by replacing h m by h mi and h s by h si throughout but presents no further di culty. This modification does not change the method and only complicates data management and notation.

Define the pressure fraction, z(h

) def = (h (u + h m ))/(h s h m ). Suppose that (•) is a bounded, smooth, monotonically increasing function which maps the interval [0, 1] ! [0, 1].
The POR at a node is defined by

c(h i ) = 8 > < > : 0 i fz(h i )  0 d i (z(h i )) if 0 < z(h i ) < 1 d i if z(h i ) 1 . ( 4
)
The inverse function of the POR , the head, h i (c) expressed as a function of outflow c at node i, will be required for the development of the method. But, the function h i (c) is not in general everywhere di↵erentiable and so in its place a multivalued, sub-di↵erential mapping is considered:

h(c i ) = 8 > > > > > < > > > > > : ; if c i < 0 ( 1, h m + u i ] i f c i = 0 (h s h m ) 1 ⇣ c i d i ⌘ + h m + u i if 0 < c i < d i [h s + u i , +1) i f c i = d i ; if c i > d i . ( 5
)
Note that, in what follows, careful distinction is made between the scalar variables, h, which represent the heads and the multivalued, sub-di↵erential mapping, h(c) which represents the inverse POR .

METHOD

The problem addressed in this paper is: find the heads, linkflows and nodal outflows in a PDM WDS which has (possibly a large number of) linkflow constraints. The PDM imposes constraints to the nodal outflows and linkflow constraints are introduced in the modelling of flow controls. The constraints are expressed as

✓ w 1 (c) w 2 (c) ◆ def = n j c n j c d !  o, ✓ w 3 (q) w 4 (q) ◆ def = n p q q max n p q min q !  o. (6)
As a formalism, any links which do not have finite constraints are treated as having 1 = q min,j  q j  q max,j = 1 (in practice only the links which have finite constraints will be treated in these vectors). Define the vector (q) whose components are the log barrier terms (Boyd & Vandenberghe 2009, §11.2.1) which constrain the link flows

(q j ) = ln( w j ) (7)
where j ranges over the components of w 3 and w 4 . Denote the content term

(c i ) = ( (h s h m ) R c i 0 1 ⇣ x d i ⌘ dx, if d i > 0, 0 i f d i = 0 denote (c) = ( (c 1 ), (c 2 ), . . . , (c n j )) T and denote ✓(c) = r c (c) = (h s h m ) 1 (c 1 /d 1 ), 1 (c 2 /d 2 ), . . . , 1 (c n j /d n j ) T . Of central inter-
est here is the restriction of the content function C(q, c), where the components of c are defined only on the intervals 0

 c i  d i . C(q, c) = np X j=1 Z q j 0 ⇠ j (s)ds a T q + c T (u + h m ) + 1 T (c) (8)
More formally, the problem addressed in this paper is:

Problem 0.1 Find the linkflows, q, and outflows, c, which minimize the objective function C(q, c) subject to the mass balance (equality ) constraints

A T q c = o (9)
and the (inequality) constraints ( 6).

The problem, restated using a log barrier technique to restrict the linkflows to the interior of their constraint intervals, is to find

min q,c ⇢ C(q, c) + 1 t (q T )1 subject to A T q c = o & o  c  d. (10)
(Note that the scalar function applied to the row vector q T here produces a row vector result.) The Lagrangian of ( 10) is

L(q, c, h, , µ) = P np j=1 R q j 0 ⇠ j (s)ds a T q + 1 T (c) + c T (u + h m ) h T (A T q + c) +µ T (c d) T c + 1 t (q T )1 subject to o, µ o. (11)
Here, h, the vector of heads, represents the Lagrange multipliers for the mass balance equality constraint and and µ are the vectors of non-negative Lagrange multipliers for the inequality constraints on the outflows, c.

Denote R = diag 1 , 2 , . . . , n j and T = diag µ 1 , µ 2 , . . . , µ n j . In view of the problem's convexity, the necessary and su cient conditions for a solution to Problem 0.1 are given by the KKT conditions:

⇠(q ⇤ ) Ah ⇤ a + 1 t r q (q ⇤ ) = o (12) h(c ⇤ ) h ⇤ L ⇤T ⇤ + U ⇤T µ ⇤ = o (13) A T q ⇤ c ⇤ = o (14) L ⇤ R ⇤ c ⇤ = o (15) U ⇤ T ⇤ (c ⇤ d) = o (16) ⇤ 0, µ ⇤ 0 . ( 17
)
where h(c ⇤ ) = ✓(c ⇤ ) + u + h m , L ⇤ is a matrix made up of rows of the identity whose indices correspond to those of the nodes at which the lower constraint is saturated (or binding) and U ⇤ is a matrix made up of rows of the identity whose indices correspond to those of the nodes at which the upper constraint is saturated. Relations ( 15) and ( 16)

represent the complementary slackness conditions, Eq. ( 12) is the conservation of energy equation and ( 13) defines h according to (5).

The Newton method for the ASM system without flow control (see Deuerlein et al.

(2019) for details) reduces (essentially) to

0 B B @ n p n j n b n p F (m) A O n j A T O S (m) b T n b e O S (m) b M b (m) 1 C C A 0 B @ q (m+1) q (m) h (m+1) h (m) c (m+1) b c (m) b 1 C A = 0 B @ ⇢ (m) e ⇢ (m) m h(c (m) b ) h (m) b 1 C A (18)
where M = r c ✓(c) is a diagonal matrix, the diagonal elements of which hold the c derivatives of the inverse POR function h(c); the matrix F denotes the diagonal matrix whose diagonal elements are q-derivatives of the corresponding terms in G(q)q; the energy and mass residuals are defined by

⇢ (m) e = ⇠(q (m) ) Ah (m) a and ⇢ (m) m = A T q (m) c (m) ;
S b is a selection matrix for the junction nodes associated with constraints that are not saturated and the subscript b indicates that the quantity in question is confined to those components which are associated with an outflow constraint that is not saturated. The The gradient of the log barrier term in ( 12) and the Hessian of ( 7) (which will be incorporated into the equivalent of ( 18)) are now discussed.

The log barrier terms

From (7), r q = P j 1 w j r q w j so r q = v = (v 1 , v 2 , . . . , v np ) T where

v j = ⇢ 1 (q j q max,j ) 1 (q min,j q j ) = 1 w 3,j 1 w 4,j if j 2 I q , 0 otherwise. ( 19
)
In view of the linearity of the constraints, r 2 q w j = 0 for w 3 and w 4 so

r 2 q = X j ✓ 1 w j r q w j ◆ ✓ 1 w j r q w j ◆ T . (20)
The j-th term in the sum (20) is

1 (q j q max,j ) 2 + 1 (q min,j q j ) 2 = 1 w 2 3,j + 1 w 2 4,j (21)
so if we denote C 3 = diag {q q max }, and C 4 = diag {q min q} then def = X j r 2 q (q j q max,j ) + r 2 q (q min,j q j ) =

0 B @ C 2 3 + C 2 4 O O O O O O O O 1 C A . (22)
A typical log barrier function (q) and its first two derivatives, r q (q) and r 2 q (q), are shown in Fig. 1. Note that r 2 q (q) > 0 and, importantly, this establishes that the IPM preserves strong convexity for this problem and this in turn guarantees that the uniqueness of the solution is preserved.

Setting the gradient of the log barrier Lagrangian (11) to zero

r B = 0 B @ ⇠(q) Ah a 1 t v A T q c h(c) h L T + U T µ 1 C A = 0 B @ o o o 1 C A (23)
gives the system of equations to be solved by the Newton method. It is evident from the form of the first block equation in ( 23) that the log barrier term can be considered a minor head loss. The Jacobian of this system is

V = 0 B @ F + 1 t A O A T O I O I M 1 C A so the system for the Newton update is now V (m) ⇣ q (m+1) T , h (m+1) T , c (m+1) T ⌘ T = r (m) B .
Denote, dropping superscripts where the context is clear,

b F = F + 1 t b ⇢ e = G(q)q Ah a 1 t v ⇢ m = A T q c ⇢ c = h(c) h L T + U T µ
It is worth noting here that the log barrier terms 1 t v and 1 t can be viewed as minor head loss terms and the vector 1 t v represents the surrogate Lagrange multipliers for the linkflow constraints at steady state (Boyd & Vandenberghe 2009, 567). As well, it is clear that only a minimal e↵ort is required to modify code to solve a PDM WDS by ASM in order to handle linkflow constraints via this technique. Figure 2 shows a typical function made up of a link head loss, ⇠(q), and the corresponding log barrier term v(q)/t. The surrogate Lagrange multipliers, ⌫ and  for head losses h 1 and h 2 are also shown.

The system to be solved is now (cf. ( 18))

0 B @ n p n j n j n p b F A O n j A T O I n j O I M 1 C A 0 B @ q h c 1 C A = 0 B @ b ⇢ e ⇢ m ⇢ c 1 C A . ( 24
)
Note that in line with the ASM, the outflow updates are applied to only those nodes in (24) can, noting that n b  n j and that S b S T b = I n b , be rewritten more compactly as

I
0 B @ n p n j n b n p b F A O n j A T O S T b n b O S b M b 1 C A 0 B @ q h c b 1 C A = 0 B @ b ⇢ e ⇢ m ⇢ c b 1 C A ( 25
)
where

M b 2 R
F q A h = b ⇢ e (26) A T q S T b c b = ⇢ m S b h + M b c b = ⇢ c b (27)
or

A T q A T b F 1 A h = A T b F 1 b ⇢ e A T q S T b c b = ⇢ m S T b M 1 b S b h + S T b c b = S T b M 1 b ⇢ c b
and adding these three equations gives

⇣ A T b F 1 A + S T b M 1 b S b ⌘ h = A T b F 1 b ⇢ e + ⇢ m + S T b M 1 b ⇢ c b . ( 28
)
Once h is found from the solution of this system, (26) gives

q = b F 1 (A h b ⇢ e ) ( 29
)
and then (27) gives

c b = M 1 b S b h ⇢ c b . ( 30
)
The updated unknowns are found as

h (m+1) = h (m) + h (m) , q (m+1) = q (m) + q (m) , c (m+1) b = c (m) b + c (m) b (31)
and the updated Lagrange multipliers are found as

(m+1) = ⇣ h (m+1) h m u ⌘ , µ (m+1) = h (m+1) h s u (32)
Relations ( 28)-( 32) form the basis of the Newton iterative method. The iterations are run until either (i) the relative di↵erence between the norms of successive iterates is su ciently small, (ii) too many iterations have been executed, (iii) there has been no reduction in the di↵erence between successive iterates after a prescribed number of steps (stalling), or (iv) The equation residuals are too small to warrant further computation because the correction terms have become too small.

The key matrix in this system is the Schur complement in Eq. ( 28),

A T b F 1 A + S T b M 1 b S b .
The term, S b , represents the ANIM for the system's pseudo links, the links which connect the virtual nodes to a ground node. In the iterative process this term acts like a Levenberg-Marquardt regularization term which helps with convergence. It is revealing to note that the Schur complement can also be written

( A T S T b ) ✓ b F 1 M 1 b ◆ ✓ A S b ◆ (33)
showing that the right-hand-side matrix is the ANIM for the extended graph.

The linkflows are required to always be in the interior of their constraint intervals so, after the linkflows have been updated, any updates which lie on the boundary of, or outside, their constraint intervals must be projected back into the interiors of their constraint intervals. Setting of the distance, ⌧ j , from the boundary to the linkflow projection point inside the constraint interval, the projection tolerance, requires some delicacy and is discussed later. The tolerance, ⌧ j , should be adjusted for the value of t but should never, as a consequence of the precision being used, be so small that the projection is onto the constraint boundary since the log barrier terms in ( 19) and ( 21) do not then exist. This is the reason that the HYAI method does not accommodate equality constraints, q min,j = q max,j . Where a linkflow is required, moreover, to take a set value this should be handled by treating that flow as known data and not handled as a tightly constrained unknown. In addition, the mass balance should be corrected to take this into account and the corresponding Lagrange multiplier can be found from the energy residuals.

Implementation Details

Starting values

The HYAI method for PDM problems itself requires starting guesses for the linkflows, heads and outflows. The initial HYAI linkflows for unconstrained links or for links with semi-infinite constraints were set to correspond to a velocity of 1/3 m/s. The finitely constrained linkflows were set to the flow rate at the midpoint of the constraint interval.

The outflows were set to one half of the nominal demand d i /2. The heads were set to

h i (d i /2
), the value of the inverse POR for an outflow which is half the nominal demand.

Projection

In IEEE double precision arithmetic, the precision used in all the calculations reported here, ✏ m ⇡ 10 16 . Taking a projection tolerance which depends only on t, or on t and the length of the constraint interval, for example choosing ⌧ = (q max q min )/t, can lead to undesirable consequences. An example illustrates this point.

Example 1 Suppose t = 10 20 , q max = 10 6 , q min = 0. Then q max q min = 10 6 and ⌧ = (q max q min )/t = 10 26 and the projection point will be computed as q max ⌧ = 10 6 10 26 = 10 6 in IEEE Double Precision arithmetic. Thus, this linkflow is projected exactly onto the constraint boundary. Using ⌧ = 1/t does not help because again for any q max larger than about 10 4 , q max ⌧ = q max in IEEE DP. ⇤

If ✏ p is a small multiple of the computer arithmetic's machine epsilon, ✏ m , then projecting to q min,j + ⌧ min,j on the left boundary and to q max,j ⌧ max,j on the right boundary with

⌧ min,j = max ✓ 1 t , ✏ p |q min,j | ◆ , ⌧ max,j = max ✓ 1 t , ✏ p |q max,j | ◆ (34)
ensures that the iterates are always projected to a point close to the boundary (and consistent with the current value of t) but still inside the constraint interval.

Continuation

It was mentioned earlier that IPMs sometimes diverge for the large values of t that are necessary to get su cient accuracy. Where this happens it may be possible to find a smaller value of t which delivers convergence and then use the solution for that t as the starting value for a larger value of t. Very few iterations are then typically required for the larger values of t. By this process of continuation it is possible to get convergence for large values of t that might diverge with the usual starting values. The stopping test for the HYAI solutions which start the ASMFC method need not use the same small stopping tolerance which is ultimately required in the ASMFC. In the large network examples discussed later, a HYAI stopping tolerance of ✏ s = 10 4 was found su cient to guarantee convergence of the ASMFC. This suggests that if the HYAI method is being used on its own and continuation is necessary, then the smaller values of t can be handled with larger stopping tolerances. In the examples discussed below the authors used a sequence of reducing stopping tolerances defined by ✏ s = 10 10 /t to progress from solutions for t = 10 12 up by factors of 10 to t = 10 20 .

An illustrative example of continuation is given in the next section. It is seen there that the first value of the penalty parameter, t = 10 12 , required six iterations to satisfy the stopping test with ✏ s = 10 10 /t = 10 2 . For the next value of the penalty parameter t = 10 13 the stopping tolerance was reset to ✏ s = 10 3 and convergence took just one iteration. Each of the subsequent values of t up to t = 10 20 , required only one iteration.

Slightly more iterations were typically required for the continuation stages with larger networks. Network N 4 , for example, took 35, 5, 1, 2, 5, 1 and 5 iterations for the values of exponent e increasing from 14 to 20, respectively, giving a total of 54 iterations.

There are some practical factors which need to be considered when using continuation:

(i) it may need some experimentation to determine the largest value of t for which direct convergence (i.e. convergence without using continuation) is possible (taking t smaller than this value necessarily leads to more iterations than are necessary) and (ii) once that value of t has been obtained it may take some more experimentation to determine the largest factor by which to increase t (again, increasing t by too small a factor un-necessarily increases the number of iterations). Finding a balance between the competing factors here is necessary. The interested reader is referred to Boyd & Vandenberghe (2009, 570 ↵.)

Accuracy A solution, b x(t), t > 0, to the problem min x ⇢ f (x) + 1 t (x) subject to g i (x)  0, i = 1, 2, . . . , m
is called a central point and every central point yields a dual feasible pair which provide a lower bound of the objective function's optimal value v ⇤ (see Boyd & Vandenberghe (2009, 565) for details). From this it follows that f (b x(t)) v ⇤  m/t and that the point b x(t) is m/t-suboptimal. Thus, to find a solution with guaranteed accuracy of the objective function ✏ f requires only that b x(t) is a solution for t = m/✏ f . In this study the authors were able to calibrate the accuracy of the HYAI method head, linkflow and outflow solutions by comparing them with those of the ASMFC. Thus, for network N 1 , using exponents e = 13-19, the linkflows and heads agree to about 5 digits and the outflows agree relatively to about 10 decimals in infinity norm. At e = 20 the linkflows agree to more than 8 decimals, the heads to more than 11 and the outflows to more than 13. These last numbers are about what one would expect with iteration stopping tolerances of ✏ s = 10 10 and led to the choice of e = 20 for all the results reported here on test networks.

Selective Updating and Restarting

Convergence can, in some cases, be assisted by the following simple restarting technique. If the linkflow constraints were being handled by an active set method, three index sets, I q l , I q b and I qu would be defined as in [START_REF] Piller | A content-based active set method for the pressure dependent model of water distribution systems with flow controls[END_REF]. Any linkflows at a constraint boundary for which the Lagrange multipliers are negative would be assigned to the set I q b , the set of linkflows which are updated during the iteration (linkflows in I q l and I qu which have non-negative Lagrange multipliers need not be updated since they already satisfy the KKT complementary slackness conditions). Using a strategy which parallels this was found useful to assist and accelerate convergence. A description of that strategy follows.

Recall that here all linkflows are in the interiors of their constraint intervals. Every link for which the flow is su ciently close to one of its constraint boundaries is identified.

Those close to the lower constraint bound are associated with the index set I q l and those close to the upper boundary are associated with the index set I qu . Their corresponding energy imbalances are then computed:  = Gq Ah a for the flows in I q l and ⌫ = (Gq Ah a) for those in I qu . The imbalances,  and ⌫, are approximate, or surrogate, Lagrange multipliers. All links in I q l with  j < 0 and links in I qu with ⌫ j < 0 then have their flows reset to the midpoints of their constraint intervals. This strategy can be considered a restarting of the iterative process and was found to assist in attaining rapid convergence. For this purpose the linkflow q j was considered to be su ciently close to a constraint boundary if its distance to the boundary was smaller than the tolerance, ✏ s (q max,j q min,j ), where ✏ s is a modest multiple of ✏ p : it was found that ✏ s needs to be at least one order of magnitude larger than the projection tolerance ✏ p in order to avoid inconsistent set assignments. The number of iterations required was approximately halved by using restarting in those case study networks where direct convergence was achieved (N 1 , N 2 , N 3 , N 5 and N 6 ) but where continuation was necessary (N 4 , N 7 , N 8 and N 9 ) the number of iterations was about the same. This restarting strategy is therefore recommended.

THE HYAI ALGORITHM

In this section the basic HYAI algorithm is described.

HYAI Algorithm Summary

Compute (i) ⌧ min,j , ⌧ max,j , for all j and (ii) the initial values q (0) , h (0) and c (0) .

The iteration loop

For m = 0, 1, 2, . . . repeat steps (a) to (j) until the stopping test is satisfied (a) Let n b denote the number of elements in

I c b . Compute the matrix M (m) b 2 R n b ⇥n b . (b) ⇢ (m) m A T q (m) c (m)
(c) Solve (28) for h (m+1) and update the heads h (m+1) using ( 31).

(d) Use ( 29) to get q (m+1) and, excluding any constrained links with (i) linkflows which are su ciently close to a boundary and (ii) have positive Lagrange multipliers, update the pipe flows q (m+1) using ( 31) .

(e) Use (30) to get c (m+1) and use (31) to update the nodal outflows c (m+1) i for which

i 2 I c b . (f) (m+1) ⇣ h (m+1) h m u ⌘ (g) µ (m+1) h (m+1) h s u (h) Update the index sets I c b , I c l , I cu . (i) c (m+1)
max(min(c (m+1) , d), o) (project any exterior outflows onto the constraint boundaries)

(j) q (m+1) max(min(q (m+1) , q max ⌧ ), q min + ⌧ ) (project any linkflows not in the interior of the constraint interval back into the interior)

It is evident that the HYAI method can be easily implemented as a simplification of the ASMFC. The log barrier terms in ( 21) are added to the first block component of the right-hand-side of the system and the terms shown in ( 19) are added to the diagonal F matrix. No damping is required and there is no need to be concerned about violating the LICQ.

Algorithmic Variations

The new method can be used in at least three di↵erent ways. If direct convergence is achieved for su ciently large t nothing more is necessary. If, however, direct convergence is not possible for large enough t, then a quickly obtained low precision solution, that can be used in one of two ways, is obtained for the largest value of t that gives direct convergence. That solution can now be used to start the HYAI method for larger t and repeating this process of continuation leads to a solution for the value of t required. Few iterations are required after the first stage of this continuation process. Alternatively, the quickly obtained low precision solution can be used to start the ASMFC method.

Again, few iterations are then required by the ASMFC. This second technique, which will be designated PHYAI to suggest a pairing of the methods, turns out to be very e↵ective, particularly on networks that challenge ASMFC. These variations are illustrated and compared in the next section.

EXAMPLES AND DISCUSSION

The new HYAI method and the PHYAI methods were applied to the following problems:

(i) a small illustrative network with n p = 10 links and n j = 8 nodes, (iv) the small network in Fig. 3. This example is explained in detail in the Appendix to help clarify the process.

The networks in (i) and (ii) were used by the authors in [START_REF] Piller | A content-based active set method for the pressure dependent model of water distribution systems with flow controls[END_REF] and their EPANET .inp files can be found in the Supplementary Materials to that paper.

Denote the measures of the errors between successive q iterates by

(m+1) q = max j |q (m+1) j q (m) j | 1 + |q (m+1) j | , ( 35
)
with corresponding notation for the errors of the c and h iterates. In all cases (i) the nominal demands were magnified by a factor of 5 to illustrate the behaviour of the methods on PDM problems (ii) the 1-side regularized Wagner POR was used (see [START_REF] Deuerlein | A content-based active set method for the pressure dependent model of water distribution systems[END_REF] for details)

(iii) the mid-interval starting scheme described above was used and

(iv) the iterations were run until the relative di↵erences between successive iterates 36) were smaller than the prescribed stopping tolerance ✏ s = 10 10 . This stopping tolerance (even though it is smaller than would normally be used in practice) was chosen to ensure that the quadratic convergence normally associated with Newton's method was evident.

(m+1) q , (m+1) c , (m+1) h (
The iterations were stopped when (i) the stopping test was satisfied, (ii) the residuals were too small to make any further computation practicable because the iterative correction terms vanished, or (iii) the iteration stalled but left residuals that are su ciently small.

The iteration counts reported below are certainly conservative: practical stopping tolerances, which would be considerably larger than 10 10 , would require many fewer iterations than are tabled below. A single iteration of the ASMFC method takes roughly the same wall-clock time as one iteration of the HYAI method.

Except where otherwise stated, all the calculations reported in this paper were performed using the authors' own codes written for Matlab (2020a) which uses IEEE standard double precision floating-point arithmetic. The tests for the case study networks N 1 -N 8

were run on a PC with an i7-4700MQ CPU and those for N 9 were run on a PC with an i9-10880H CPU.

Small Illustrative Network Example

Figure 4 shows the small illustrative network and its steady-state solution when all links are essentially unconstrained (none of the constraints are saturated) while Fig. 5 shows the same network when three constraints, all of which are active at steady-state, are applied. The constraints on links 3 and 9 are acting as FCVs while the flow in link 10 is at the lower link flow constraint boundary and, since the corresponding multiplier in this case can be interpreted as a negative head loss, this constraint can be seen to be acting as a pump. The continuation technique is illustrated in Table 1. The initial value of t, which corresponds to an exponent e = 12, takes six iterations to satisfy the stopping test and for each subsequent value of t, each increasing by one order of magnitude, takes only one iteration. For the last entry in the table, where t = 10 20 , HYAI and ASMFC solutions agree relatively to about 14 decimal digits but for smaller values of t they match to fewer digits.

HYAI and PHYAI Methods on Larger Networks

The results of applying both the HYAI method and the PHYAI method (the paired HYAI and ASMFC methods) to the nine case study networks are now discussed. In the case of the PHYAI method there are two cases to consider: (i) direct convergence for the required t or (ii) a continuation solution for the required t. In the case of direct convergence, the HYAI iterations were stopped if either (i) the (more relaxed) stopping test with ✏ s = 10 4 was satisfied or (ii) more than 12 HYAI iterations had been executed. The solution thus produced was then used to start the ASMFC with ✏ s = 10 10 . Where continuation was required, the graded tolerances scheme which relates the stopping tolerance to the value of t was used to terminate the HYAI stage of the PHYAI. As reported earlier, the HYAI solution for the larger values of t required few iterations.

Columns two to four of Table 2 show the number of links, nodes and sources. Every cotree link for which the unconstrained flow was positive was constrained to [0, +10%] of that flow while if the flow negative it was constrained to [ 10%, 0] of its unconstrained value. Any constraints for which the width of the constraint interval was smaller than 0.001 L/s was culled from the collection. The number of constraints that remained is shown in column five. Column 6 shows the total nodal delivery as a percentage of the total nominal demand.

Column seven gives the number of HYAI iterations required to satisfy the stopping test (36) with ✏ s = 10 10 and columns 8-10 show the three quantities h q,h,c def = log 10 q,h,c

for the HYAI method. They represent the number of decimal digits of agreement between the last two HYAI iterates as measured by the relative error in (36).

Five of the networks N 1,2,3,5,6 , converged directly with between 12 and 23 iterations for e = 20. Networks N 4,7,8,9 , on the other hand, all required continuation starting from e = 12 (to achieve an e = 20 solution) because, if started with e = 20, they had not satisfied the stopping test after more than 70 iterations.

Both HYAI and PHYAI methods took about the same number of iterations for networks N 1,2,3,5,6 but for networks N 4,7,8,9 the HYAI method took about twice as many iterations as PHYAI because of continuation. Column 11 shows the number, k p h , of iterations required by the HYAI stage of the PHYAI scheme and column 12 shows the number, k p a , of iterations required for the ASMFC stage of the PHYAI scheme. Most networks required 3-7 ASMFC iterations, one took 10 and the largest network, N 9 , took 18. Column 13 shows the sum, k p h + k p a and the following three columns show the numbers a q,h,c which represent the number of decimal digits of agreement between the last two ASMFC PHYAI iterates as measured by the relative error in (36).

Column 17 shows the quantity ↵ q which is defined by

↵ q def = log 10 ⇢ max j |q h,j q a,j | 1 + |q a,j |
and columns 18 and 19 show the corresponding quantities, ↵ h for h and ↵ c for c. These numbers indicate, in decimal digits, the degree of agreement between the solutions of the HYAI method and those of the PHYAI and they are entirely consistent with a stopping tolerance of 10 10 .

Thus, the PHYAI is the recommended method where speed is an issue in view of the smaller number of iterations it required. But, one could start by using the HYAI method and, in the event that continuation is required, switch to ASMFC to achieve convergence.

IPOPT AND EPANET ON THESE PROBLEMS

Exactly the same problems used to demonstrate the PHYAI and HYAI methods on networks N 1 -N 9 were solved using the Matlab Toolbox written by P. Carbonetto and updated by [START_REF] Bertolazzi | ebertolazzi/mexipopt[END_REF] for the o↵-the-shelf package, IPOPT (Wachter & Biegler 2006). The objective function ( 8) subject to the equality constraints ( 9) and the variable bounds (6) was optimized using IPOPT. The comparison with the authors' Matlab codes was performed to provide some indication of the e ciency of the new method.

The 'mumps' linear solver and the 'adaptive' 'mu strategy' options were used (see https://coin-or.github.io/Ipopt/OPTIONS.html for details). In all cases the IPOPT termination status for the stopping tolerance 10 10 (the same as was used on the PHYAI and HYAI codes) was shown as "Optimal solution found". The results of these tests are shown in Table 3. Column 2 of Table 3 shows, k I , the number of IPOPT iterations required to satisfy the termination test and column 3 shows the corresponding PHYAI number, k P .

Column 3 shows the optimal value of the objective function, C(q, c), column 5 shows the final delivery fraction, ⇣ as a percentage. Column 6 shows, s I , the number of wall-clock seconds for IPOPT to satisfy the stopping test and column 7, shows the number of PHYAI wall-clock seconds to satisfy the stopping test. Column 8 shows the ratio of these times.

The PHYAI codes took fewer iterations than the IPOPT codes in all cases. But an iteration count comparison in a case such as this does not necessarily tell the whole story.

Here, wall-clock times are more revealing. The IPOPT method takes about the same wall-clock time as the PHYAI method in two cases, about 2-4 times as long in five others and 7 and 9 times as long in the remaining two cases. For the largest network the savings in real time are significant: about 4 minutes for PHYAI as opposed to about 34 minutes for the IPOPT codes.

The nine test networks were also run with the latest release of EPANET (version 2.2 from July, 23rd, 2020) which includes the simulation of pressure dependent demands. Flow control is, to the knowledge of the authors, treated by valve-state heuristics as in previous versions. In general, the results are promising for the smaller test networks (N 1 -N 6 ).

For networks N 7 , N 8 and N 9 EPANET fails to converge after 500 iterations. This may be because of the heuristics used for tackling flow control devices. The status report file shows a huge number of state switches over the iterations. Another problem is that EPANET cannot handle positive lower bounds on flows which are required in the design and control of pumps. It is restricted to FCV upper bounds and check valves (one-directional flows, q > 0).

CONCLUSIONS

The authors' ASMFC, which solves for the heads, linkflows and outflows of a PDM WDS, uses an active set method to enforce both the PDM constraints on the outflows and the modelling constraints on the linkflows. The method has been shown to be rapid and robust when applied to eight case study networks, N 1 -N 8 , with up to 60 linkflow constraints, some of which are equality constraints. Deeper investigations revealed that when every one of the (up to ⇡1700) links in the networks' co-trees were constrained to ±10% of their unconstrained values, six of the networks converged rapidly and without di culty with the ASMFC but two of the networks, N 4 and N 7 , failed to converge.

A new hybrid method, HYAI, which combines active set principles and log barrier IPM principles has been presented which overcomes these convergence di culties: the constraints on outflows inherent in PDM WDSs are handled by an ASM and the linkflow constraints required to model flow control devices are handled by a log barrier IPM. Implementation details concerning projection, continuation, restarting and selective updating have been outlined which improve the performance of the method. The following points are noted.

(a) Unlike external penalty methods HYAI does not require the computational burden of damping but it may require continuation which typically increases the number of iterations.

(b) Whereas with external penalty methods it can be very hard to cause iterates which move far away from the constraint interval back to approach the constraint interval, the HYAI linkflow iterates are always within the constraint intervals.

(c) Implementing HYAI requires only a very simple modification of the ASMFC method.

The new HYAI method and the continuation strategy were demonstrated on a small illustrative network and then used to find the steady-state heads, linkflows and outflows of all eight previously used case study networks and one other, much larger network. The HYAI method achieved direct convergence with five of the nine networks but the other four required continuation. Thus, the HYAI method succeeded in finding steady-state solutions where the ASMFC had not. But on those networks which challenged the ASMFC and some others, the HYAI method on its own took about twice as many iterations as the PHYAI. Pairing the HYAI and ASMFC methods, the PHYAI method, which uses low precision HYAI solutions to start the ASMFC was found to solve the challenging systems in fewer iterations overall. Three di↵erent approaches are available: (i) attempt a solution with HYAI seeking direct convergence. If that fails (ii) use a low precision HYAI solution to start the ASMFC or (iii) use the low precision solution to initiate a continuation process that leads to a HYAI solution for the required t. An alternative is to begin by finding a low precision HYAI solution and using that to start ASMFC. This strategy will reduce the average number of iterations required over a range of networks. The new method was also illustrated in detail on a small example network to help clarify the method for readers.

The same problems that were used to illustrate the new method were solved with the freely available interior point optimization package IPOPT. The IPOPT method takes about the same wall-clock time as the PHYAI method in 2 cases, about 2-4 times as long in 5 others and 7 and 9 times as long in the remaining 2 cases. For the largest network the savings in real time are significant: about 4 minutes for PHYAI as opposed to about 34 minutes for the IPOPT codes.

The authors would, in future work, like to extend this technique to handle pressure control problems and to more deeply investgate the use of capacity analysis in WDS design and in optimal network control.

SUPPLEMENTARY MATERIAL

Network parameter summary statistics and summary solution statistics for those networks which are in the public domain (N 1 , N 3 , N 4 and N 7 ) are avaiable from the ASCE as material supplementary to this paper.

DATA AVAILABILITY STATEMENT

Some or all data, models, or code generated or used during the study are proprietary or confidential in nature and may only be provided with restrictions (e.g. anonymized data). The four networks N 2 , N 5 , N 6 and N 8 are not available either because they are proprietary or because of security concerns.

Some or all data, models, or code used during the study were provided by a third party. (N 9 ). Direct request for these materials may be made to the provider as indicated in the Acknowledgments.

EPANET .inp files for the network shown in Fig. 4 and networks N 1 , N 3 , N 4 and N 7

are available from the ASCE library (www.ascelibrary.org) as material that is supplementary to [START_REF] Deuerlein | A content-based active set method for the pressure dependent model of water distribution systems[END_REF].

In addition, EPANET .inp files for the network shown in Fig. 3

if ✓ A S b ◆ x = ✓ o o ◆ ) x = o
and this is immediate since the first block equation is Ax = o and the full rank of A ensures that this is so. In fact the matric b C and the LICQ matrix C have the same rank, as can be seen by the following argument. Recall that the HYAI LICQ matrix is C are both determined by the rank of S a A T : the rows of A T that are selected by S a must form a linearly independent set. This will always be true if A has full rank but it must hold if A has less than full rank: the set of active POR constraints must be such that the corresponding rows of A T form a linearly independent set. Thus, the LICQ is never violated in the ASM solution of a PDM problem if the NAIM has full rank and there are no linkflow constraints. If LICQ is not violated then unique Lagrange multipliers exist and a feasible solution exists.

C = 0 B @ n p n j n j A T I n c l O S c l n cu O S cu 1 C A Denote S T a = ( S T c l S T cu ). Then, we can rewrite C = ✓ A T I O S a ◆ C has full rank if ✓ A O I S T a ◆ ✓ x 1 x 2 ◆ = ✓ o o ◆ ) x 1 = o, x 2 = o.

AN EXAMPLE IN DETAIL

The flows are, in what follows, quoted in L/s except where otherwise stated. The network shown in Fig. 3 The outflow constraints (6) for this case are c 1,2  0, c 1  10, c 2  15 and the linkflow constraints are q 1 1  0 and q 1  0. Thus, w 1 and w 2 each have two components and w 3 and w 4 each have one. At start all of the outflows are at the midpoints of their constraint intervals so none of their constraints are saturated: c 1 = 5, c 2 = 7.5. Thus, S

b in ( 18) is an identity matrix. The initial flow in link 1 is q 1 = 0.5 (the midpoint of its constraint interval) so its constraint is also not saturated at start. The other two links have initial flows, q j , j = 1, 2, given by q j = ⇡ 12 D 2 j = 16.4 which correspond to velocities of 1/3 m/s. The initial point is therefore inside the problem's feasible region. The KKT conditions ( 12) to (17) now read

⇠ 1 h 1 + h 2 1 t v 1 = 0 ⇠ 2 + h 1 15 = 0 ⇠ 3 + h 2 15 = 0 h(c 1 ) h 1 1 + µ 1 = 0 h(c 2 ) h 2 2 + µ 2 = 0 q 1 + q 2 c 1 = 0 q 1 + q 3 c 2 = 0 1 c 1 = 0 2 c 2 = 0 µ 1 (c 1 d 1 ) = 0 µ 2 (c 2 d 2 ) = 0 1,2 0, µ 1,2 0
The system (18) is, dropping iteration subscripts and noting that since no constraints are saturated the Lagrange multipliers 1,2 = µ 1,2 = 0,

0 B B B B B B B B B B B @ f 1 -1 1 0 0 f 2 1 0 0 0 f 3 0 1 0 0 -1 1 0 0 0 -1 0 1 0 1 0 0 0 -1 0 0 0 -1 0 h 0 (c 1 ) 0 0 0 0 -1 h 0 (c 2 ) 1 C C C C C C C C C C C A 0 B B B B B B B B B B @ q 1 q 2 q 3 c 1 c 2 h 1 h 2 1 C C C C C C C C C C A = 0 B B B B B B B B B B @ ⇠ 1 h 1 + h 2 ⇠ 2 + h 1 15 ⇠ 3 + h 2 15 q 1 + q 2 c 1 q 1 + q 3 c 2 h(c 1 ) h 1 h(c 2 ) h 2 1 C C C C C C C C C C A
where the f i represent the frictional head loss derivatives on the diagonal of the matrix F in (18). We remind the reader that h(c) here is the inverse POR defined in (5).

The log barrier terms of ( 19) and ( 21) now need to be added to this system to get the system which is to be solved. Thus, f 1 is to be replaced by b

f 1 = f 1 + 1 t z 1 = f 1 + 1 t(q 1 1) 2 + 1 tq 2 1
and the first component of the right-hand-side vector now has 1 t v 1 = 1 t(q 1 1) + 1 tq 1 subtracted from it, giving the system to be solved at this step as

0 B B B B B B B B B B B @ f 1 + 1 t z 1 -1 1 0 0 f 2 1 0 0 0 f 3 0 1 0 0 -1 1 0 0 0 -1 0 1 0 1 0 0 0 -1 0 0 0 -1 0 h 0 (c 1 ) 0 0 0 0 -1 h 0 (c 2 ) 1 C C C C C C C C C C C A 0 B B B B B B B B B B @ q 1 q 2 q 3 c 1 c 2 h 1 h 2 1 C C C C C C C C C C A = 0 B B B B B B B B B B @ ⇠ 1 h 1 + h 2 1 t v 1 ⇠ 2 + h 1 15 ⇠ 3 + h 2 15 q 1 + q 2 c 1 q 1 + q 3 c 2 h(c 1 ) h 1 h(c 2 ) h 2 1 C C C C C C C C C C A
The modification required to the system in (18) to get the system (25) for solution is thus seen to be minimal.

Rather than solve this system as a whole the smaller dimension system of ( 28) is solved for the heads, h 1 , h 2 , and then ( 29) and ( 30) are used to find the corresponding linkflows, q 1 , and outflows, c 1 , c 2 . Once the updated unknowns are found using (31), the Lagrange multipliers, 1,2 and µ 1,2 are found using (32). Of course, there are no Lagrange multipliers for the linkflow constraints.

The HYAI and ASMFC phases each took 4 iterations to satisfy the stopping test with a stopping tolerance of 10 14 when t = 10 10 . The final heads were h 1 = 14.92 m and h 2 = 14.88 m. The final outflows were c 1 = 8.64 L/s and 12.94 L/s. The final linkflows q 1,2,3 were 1.00, 9.64, 11.94 L/s, respectively so the flow constraint on link 1 was saturated at the upper bound. Its final Lagrange multiplier was ⌫ 1 = 0.0378. The set assignments did not change after the first iteration.

TABLES

  with unknown-head by d = (d 1 , d 2 , . . . , d n j ) T 2 R n j and denote by n d the number of nodes with non-zero demands. Denote by c(h, d) 2 R n j the vector whose elements are the POR function values at the n j nodes of the system and denote 1 = (1, 1, . . . , 1) T . Throughout

  index set for these components is, following the notation in[START_REF] Piller | A content-based active set method for the pressure dependent model of water distribution systems with flow controls[END_REF], denoted I c b and it has n b elements. Similarly, denote by I q the set of indices which represent the links with finite constraints and denote by n Iq the number of linkflow constraints.

  c b . Denote by S b 2 R n b ⇥n j the selection matrix made up of the rows of an n j identity which correspond to those nodes with positive nominal demand and indices in I c b . Eq.

  n b ⇥n b is the diagonal matrix with the derivatives of the inverse POR functions for nodes in I c b on its diagonal and c b 2 R n b is the vector of unsaturated outflows with d i > 0 and ⇢ c b is the outflow residual computed only for those outflows which are not subject to active constraints. The system for the hybrid method is now seen to be obtained from the ASM in Deuerlein et al. (2019) by simply adding the log barrier penalty term 1 t to the matrix F on the left and the log barrier term 1 t v to ⇢ e on the right. The block equations are b

(

  ii) eight case study networks N 1 -N 8 [four of which are available online in the ASCE Library (www.ascelibrary.org)], and which have between 934 and 19,647 links and between 848 and 17,971 nodes and (iii) a larger network with n p = 157, 044, n j = 150, 630, and four sources, n f = 4 (Sitzenfrei et al. 2020).

  The block equations are Ax 1 = o and x 1 + S T a x 2 = o. Assume these equations hold. Multiplying the second block equation on the left by A gives AS T a x 2 = o and x 2 = o i↵ AS T a has full rank. Now, S a b C = ( S a A T O ) by the orthogonality of S a ,and S b . So the ranks of C and b

  was solved using the PHYAI method to illustrate the workings of the method in detail. It has ANIM all links have length 500 m, diameters 250 mm, roughnesses 0.03 mm, (b) nodes 1 and 2 have zero elevation and the tank at node 3 has an elevation of 15 m, (c) node 1 has demand 10 L/s and node 2 has demand 15 L/s, (d) the flow in L/s through link 1 when no links are constrained is q 1 = 2.0023 L/s, (e) the flow in link 1, the network's only co-tree link, is constrained to lie between zero and 1 L/s, (f) the Wagner POR was chosen for this example(Wagner et al. 1988), (g) the minimum pressure head, h m = 0, and the service pressure head, h s = 20 m.

  Figure captions.

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure2: A representational plot of a function made up of a link head loss, ⇠(q), together with the corresponding log barrier term, v(q)/t.

  and networks N 1 , N 3 , N 4 and N 7 with the flow constraints used in the tests reported here (flow control valves and check valves) are available from the ASCE library as material that is supplementary to this paper. Suribabu, C. R., Renganathan, N. T., Perumal, S. & Paez, D. (2019), 'Analysis of water distribution network under pressure-deficient conditions through emitter setting',

	Drink. Water Eng. Sci. 12, 1-13.
	Tabesh, M. (1998), Implications of the Pressure Dependency of Outflows on Data Man-
	agement, Mathematical Modelling and Reliability Assesment of Water Distribution
	Systems, PhD, Department of Civil Engineering, University of Liverpool, UK.
	Todini, E. & Pilati, S. (1988), A gradient algorithm for the analysis of pipe networks.,
	John Wiley and Sons, London, pp. 1-20.
	Wachter, A. & Biegler, L. T. (2006), 'On the implementation of a primal-dual interior
	point filter line search algorithm for large-scale nonlinear programming', Mathematical
	Programming 106(1), 25-57.
	Wagner, J., Shamir, U. & Marks, D. (1988), 'Water distribution reliability: sim-
	ulation methods', J. Water Resour. Plann. Manage. 114(3), 276-294. DOI:
	10.1061/(ASCE)0733-9496(1988)114:3(276).
	APPENDIX
	UNIQUNESS OF THE SOLUTION AND THE LICQ
	The matrix b C = ( A T	S T b ) of (33) has full rank if A has full rank. This fact is
	established

Wood, D. & Charles, C. (1972)

, 'Hydraulic network analysis using linear theory', Journal of the Hydraulics Division 98(7), 1157-1170.

Table 1 :

 1 Convergence data for the HYAI method with continuation applied to the network shown in Fig.4.

	k	e	(k) h	(k) q	(k) c
	1 12 6.298900e 01 4.506700e 01 2.839100e 01
	2 12 5.112000e 01 8.036900e 01 9.616300e 02
	3 12 1.953100e 01 5.618000e + 00 1.234900e 01
	4 12 9.981000e 02 6.948300e 01 4.075300e 02
	5 12 1.918700e 02 5.451400e 03 7.863300e 03
	6 12 7.781600e 04 2.985700e 04 7.923300e 04
	1 13 7.309900e 05 2.684000e 05 9.161300e 05
	1 14 3.206800e 06 1.181400e 06 3.846000e 06
	1 15 6.226300e 09 2.295200e 09 7.469500e 09
	1 16 5.038500e 14 5.592700e 13 6.243100e 13
	1 17 5.012200e 15 5.590300e 14 6.255000e 14
	1 18 6.415700e 16 5.358600e 15 5.957200e 15
	1 19 6.490100e 15 7.314800e 15 7.148600e 15
	1 20 6.497500e 14 7.312400e 14 7.386900e 14

Table 2 :

 2 The parameters for the nine case study networks, the convergence results for the HYAI and PHYAI methods and the agreement between the solutions for the two methods.

	HYAI	PHYAI

†

required continuation to achieve a solution for e = 20.

Table 3 :

 3 PHYAI and IPOPT iteration numbers and wall clock times compared.

	ID k I k P	C ⇤ (q, c)	⇣%	s I	s P	s I /s P
	N 1 33 14	8.931 69.8	1.2	0.5	2.4
	N 2 34 15	0.229 15.3	0.7	0.3	2.3
	N 3 51 20	67.064 51.8	2.9	0.4	7.2
	N 4 39 25	169.188	7.8	0.9	0.8	1.1
	N 5 16 14	1.023 26.4	0.4	0.3	1.3
	N 6 34 18	1.690 20.0	4.3	1.2	3.6
	N 7 50 21	43.346 12.3	11.8	3.0	3.9
	N 8 39 22	27.563 25.9	17.5	4.9	3.6
	N 9 80 33	123.043	5.1 2016.0 221.7	9.1

ACKNOWLEDGMENTS

The authors thank Univ.-Prof. Dipl.-Ing. Dr. Robert Sitzenfrei who provided the authors with the virtRom network, N 9 . This network is available from https://www.uibk.ac.at/umwelttechnik/softwareanddatasets/. One author gratefully thanks Dr J. Kautsky for useful discussions. One author of this paper was supported in part by the German Ministry for Education and Research (BMBF Project W-Net 4.0 02WIK1477C).

Figure 1: Representational plots of the log barrier function (q) of ( 7) and its first two