Sylvan Elhay

Olivier Piller

Jochen Deuerlein

Angus R Simpson

A Robust, Rapidly Convergent Method That Solves the Water Distribution Equations For Pressure-Dependent Models

Keywords: pressure dependent models, consumption functions, water distribution systems, cocontent, least squares residuals, Goldstein algorithm

In the past, pressure dependent models (PDM) have suffered from convergence difficulties. In this paper conditions are established for the existence and uniqueness of solutions to the PDM problem posed as two optimization problems, one based on weighted least squares (WLS) and the other based on the co-content function. A damping scheme based on Goldstein's algorithm is used and has been found to be both reliable and robust. A critical contribution of this paper is that the Goldstein theorem conditions guarantee convergence of our new method. The new methods have been applied to a set of eight challenging case study networks, the largest of which has nearly 20,000 pipes and 18,000 nodes, and are shown to have convergence behaviour that mirrors that of the Global Gradient Algorithm on demand dependent model problems. A line search scheme based on the WLS optimization problem is proposed as the preferred option because of its smaller computational cost. Additionally, various consumption functions, including the Regularized Wagner function, are considered and four starting value schemes for the heads are proposed and compared. The wide range of challenging case study problems which the new methods quickly solve suggests that the methods proposed in this paper are likely to be suitable for a wide range of PDM problems.

INTRODUCTION

Water engineers are frequently required to find the hydraulic steady-state pipe flows and nodal heads of a water distribution system (WDS) model by solving a set of non-linear equations. In practice, the water demand components arise as a combination of various sources (such as showers, washing machines, toilets and garden use). The Demand Dependent Model (DDM) requires the delivery of the prescribed demands regardless of the available pressure or head. This requirement can lead to solutions that are mathematically correct but not physically realizable. For example, if the pressure at a node drops below a certain level, then the demand required at that node cannot be delivered. These failures are characterized by a mismatch between the demand and the available pressure at a node and they led to the development of the pressure dependent model (PDM). In PDMs there is a pressure-outflow relationship (POR) which determines the flow or delivery at a node.

There is a wide variety of approaches that have been tested in the search for suitable PDMs and fast, reliable methods to solve the resulting model equations. Early attempts to include pressure dependence in WDS analysis [START_REF] Bhave | Node flow analysis distribution systems[END_REF] modelled the dependence of flow on pressure by the (discontinuous) Heaviside function: the set demand, d, is delivered if the pressure is greater than a prescribed service pressure head, h s or it is zero if the available pressure head is below h s . [START_REF] Wagner | Water distribution reliability: simulation methods[END_REF], and later [START_REF] Chandapillai | Realistic simulation of water distribution system[END_REF], avoided the discontinuities in Bhave's model by proposing a continuously varying model in which the flow delivery is proportional to the square root of the pressure. The choice of the square root curve was based on a flow model that applies to a single, circular aperture. Both [START_REF] Bhave | Node flow analysis distribution systems[END_REF] and [START_REF] Tabesh | Implications of the Pressure Dependency of Outflows on Data Management, Mathematical Modelling and Reliability Assesment of Water Distribution Systems[END_REF] proposed solving the PDM problem by using a two-step iterative procedure.

Here the problem is repeatedly (i) solved as a DDM model and then (ii) the demands are corrected according to the DDM solution heads and a chosen PDM relationship. In another development, the computer modelling package EPANET 2 [START_REF] Rossman | EPANET 2 Users Manual, Water Supply and Water Resources Division[END_REF]) allowed users to model leakages with the power equation by defining emitters at nodes. This introduced some degree of pressure-dependent modelling by having the solver add artificial reservoirs, the elevations of which are used to calculate the emitter outflows using the DDM solver.

The practical work on pressure dependent modelling up to this point was provided with a firm theoretical underpinning when [START_REF] Deuerlein | On the hydraulic system analysis of water supply networks[END_REF] showed that for almost all of the relevant elements of a WDS model (including PDM nodes) a strictly monotone subdifferential mapping between flow and head loss can be identified, ensuring that the corresponding content and co-content functions are strictly convex and thereby guaranteeing uniqueness of the solution. Existence of the solution was established by showing that the feasible set, which is described by a system of linear equalities and inequalities, is not empty. [START_REF] Todini | A more realistic approach to the "extended period simulation" of water distribution networks[END_REF] proposed a PDM technique that does not require the introduction of a POR. The procedure uses three steps: (i) the DDM solution is determined, (ii) the pressure at any node with DDM pressure less than the service pressure is fixed and the maximum demand compatible with this constraint is calculated and (iii) the results of the second step are used to build a PDM solution which is similar to the Heaviside function. The author provided an example for which this method finds a solution while the EPANET emitters model approach for this problem fails.

The system of equations for the PDM problem can be formulated in a way which runs parallel to the DDM problem formulation by including a POR element in the continuity equation. This fact was observed by [START_REF] Cheung | Extension of EPANET for pressure driven demand modeling in water distribution system[END_REF] and [START_REF] Wu | Extended global-gradient algorithm for pressure-dependent water distribution analysis[END_REF]. In an attempt to avoid using a POR, [START_REF] Ang | Solution for water distribution systems under pressure-deficient conditions[END_REF] progressively introduced a set of artificial reservoirs into the network to initiate nodal outflows. These outflows are adjusted to lie between zero and the design demand, d. Even so, this heuristic method is very time-consuming and it is, in fact, equivalent to using the Heaviside POR. Some authors (e.g. [START_REF] Lippai | Demand constructs for risk analysis[END_REF] introduced artificial check valves and artificial flow control valves to address reverse flows associated with artificial reservoirs. This approach has several shortcomings, not least of which is the fact that it involves a change in the network topology and typically increases the dimension of the problem. The consequent increase in computation time for large networks constitutes a serious disadvantage [START_REF] Wu | Extended global-gradient algorithm for pressure-dependent water distribution analysis[END_REF]. [START_REF] Giustolisi | Pressure-driven demand and leakage simulation for water distribution networks[END_REF] (and later Siew & Tanyimboh (2012) among others) recognized that adding a POR function introduced convergence problems not seen in the Global Gradient Algorithm (GGA) of [START_REF] Todini | A gradient algorithm for the analysis of pipe networks[END_REF] applied to the DDM problem. In an attempt to avoid cycling, [START_REF] Giustolisi | Pressure-driven demand and leakage simulation for water distribution networks[END_REF] used an overrelaxation parameter to correct both pipe flow and nodal head iterates where the heuristic used the L 1 norm to choose a step length. [START_REF] Siew | Pressure-dependent EPANET extension[END_REF] proposed a backtracking and line search heuristic but they corrected only the heads and not the flows. [START_REF] Giustolisi | Demand components in water distribution network analysis[END_REF] published a comprehensive study for the classification of demands in a WDS. They identified four major groups of demands (human based, volume controlled, uncontrolled orifices and leakage) and considered demand models, each type of which has its own special pressuredemand relationship. In addition, they discussed the effect of steady-state assumptions (and extended period simulations) for realistic stochastically pulsed demands and they introduced a pipe leakage model dependent on the average pipe pressure using a Fixed and Variable Area Discharge (FAVAD) technique. 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 which, in form at least, comes from the approach of [START_REF] Bhave | Node flow analysis distribution systems[END_REF]: the PDM problem is attacked by repeatedly solving the corresponding DDM problem with the GGA and adjusting the demands after each solution. They implemented their method as an extension to EPANET and then used it to compare the effects, on the solutions, of using each of four different consumption functions. However, [START_REF] Jun | Iterative methodology of pressure-dependent demand based on EPANET for pressure-deficient water distribution analysis[END_REF] made no recommendations about which of the consumption functions should be used. Some authors [START_REF] Piller | Incorporating the FAVAD leakage equation into water distribution system analysis[END_REF] used the power equation or the FAVAD pressure-dependent leakage equation at nodes with leakage to model the dependence of flow on pressure. [START_REF] Muranho | Pressure-dependent demand and leakage modelling with an EPANET extension WaterNetGen[END_REF] discussed the package WaterNetGen in which the reference pressure head of each node is set as a user-defined function. They reported that "the embedding of POR into the hydraulic solver creates some difficulties for convergence". The fact that there are so many different approaches to the PDM problem underlines the fact that existing algorithms for the PDM problem have some important limitations.

In this paper, a model in which the continuity equation includes a POR component is solved by a variation of the PDM counterpart of the GGA for the DDM problem. A Newton method is used in which, at each iteration, a linear system is solved for the heads and then the flow rates are updated using the equations for energy conservation. This method is sometimes referred to as the PDM extension of GGA. Some regularization of the POR function may be required to ensure the continuity of its first derivatives. The documented poor convergence, or even divergence, of the undamped PDM counterpart to the GGA for DDM problems is illustrated on a small network. It is shown that a new (fourth) formulation of the PDM problem, the Weighted Least Squares (WLS) optimization formulation, is equivalent to three known (equivalent) PDM formulations. The conditions for the existence and uniqueness for the WLS formulation follow. Two of the four equivalent optimization problems, the co-content (CC) and WLS versions, satisfy the conditions of a theorem due to [START_REF] Goldstein | Constructive Real Analysis[END_REF] and Gauss-Newton methods with Goldstein's line search algorithm based on those two formulations are then proposed. An important development is that using Goldstein's algorithm on the CC and WLS formulations of the optimization problems mathematically guarantees convergence.

The new methods are both robust and rapidly convergent. The effectiveness of the WLS and CC methods are demonstrated on eight benchmark water distribution network problems, the largest of which has almost 20,000 pipes and 18,000 nodes. The damped Gauss-Newton method with Goldstein's line search is shown to have convergence behavior that mirrors that of the GGA applied to DDM problems. Two modelling choices associated with the PDM are also discussed in this paper: (i) the POR or consumption function, (ii) the starting values that are needed when solving PDM problems. A weighting scheme that is necessary to ensure numerical balance between heads and flows used in the objective function is proposed and evaluated. A cubic polynomial consumption function, first introduced by [START_REF] Fujiwara | Reliability assessment of water supply systems with storage and distribution networks[END_REF], is considered and its effect is compared with that of the Regularized Wagner consumption function of [START_REF] Piller | Incorporating the FAVAD leakage equation into water distribution system analysis[END_REF].

DEFINITIONS AND NOTATION

Consider a water distribution system (WDS) that has n p pipes and n j nodes at which the heads are unknown. Denote by q = (q 1 , q 2 , . . . , q np ) T ∈ R np the vector of unknown flows in the systems and by h = (h 1 , h 2 , . . . , h n j ) T ∈ R n j the unknown heads at the nodes in the system. Let n f ≥ 1 denote the number of reservoirs or fixed-head nodes in the system, let A 1 denote the n p × n j , full rank, unknown-head node-arc incidence matrix, let A 2 denote the node-arc incidence matrix for the fixed-head nodes and let e denote the water surface elevations of the fixed-head nodes. Furthermore, denote by G(q) ∈ R np×np the diagonal matrix whose diagonal elements are such that the components, δh j (q j ), of the vector G(q)q are monotonic and of class C 1 and which represent the pipe head losses in the system (often modeled by the Hazen-Williams or Darcy-Weisbach formulae). Denote the vector [pdmjr3Figs-singl.tex: 18:46, January 19, 2016] of the desired demands at the nodes with unknown-head by d

= (d 1 , d 2 , . . . , d n j ) T ∈ R n j .
The PDM is constructed in such a way that the flow delivered at a node is determined by the pressure head at that node. Denote by h m the minimum service head (which is the sum of the minimum pressure head and the elevation head), and denote by h s the service head (which is the sum of the service pressure head and the elevation head). Suppose that γ(h) is a bounded, smooth, monotonically increasing function which maps the interval [h m , h s ] → [0, d]. The consumption function, c(h), is a function that maps the pressure head to delivery:

c(h) = 0 if h ≤ h m γ(h) if h m < h < h s d if h ≥ h s
Thus, if the pressure at a node lies between h m and h s , then the flow, or delivery, at that node lies somewhere between 0 and the set demand, d. Nodes at which the pressure head is h m or less have zero flow and those at which the pressure head is h s or greater get full delivery, d. Unlike the DDM, the PDM delivers only the flow that the solution pressure heads can provide, a feature that has spurred considerable interest in modelling pressure dependence. Denote by c(h) ∈ R n j the vector whose elements are the consumption functions at the n j nodes of the system. It is assumed in this study, and without loss of generality, that all nodes have the same values of h m and h s and the same consumption curve, γ(h). Any nodes at which the delivery is zero are said to be in failure mode. Nodes at which the delivery is between zero and d are said to be in partial delivery mode and nodes which have full delivery are said to be in normal mode.

WDS PDM EQUATIONS

The steady-state flows and heads in a WDS with PDM are usually found as the zeros of the nonlinear system of the n p + n j equations

f (q, h) = G(q)q -A 1 h -a -A T 1 q -c(h) = o, (1)
where a = A 2 e . A natural way to approach the solution of (1) is to use a Newton iteration based on the Jacobian of f ,

J (q, h) = F (q) -A 1 -A T 1 -E(h) , (2)
where F (q) and E(h) are diagonal matrices which are such that (i) the terms on the diagonal of F (q) are the q-derivatives of the corresponding terms in G(q)q and (ii) the terms on the diagonal of E are the h-derivatives of the corresponding terms in c(h). It is assumed in what follows that the diagonal terms of F and E are non-negative. Denote the energy and continuity residuals of (1) by

ρ e = G(q)q -A 1 h -a, ρ c = -A T 1 q -c(h). (3)
The Newton iteration for (1) proceeds by taking given starting values q (0) , h (0) and repeatedly computing, for m = 0, 1, 2, . . ., the iterates q (m+1) and h (m+1) from

F (q (m) ) -A 1 -A T 1 -E(h (m) ) q (m+1) -q (m) h (m+1) -h (m) = - ρ (m) e ρ (m) c
until, if the iteration converges, the relative difference between successive iterates is sufficiently small. In what follows the Jacobian J (m) will be denoted simply by J where there is no ambiguity. The January 19,2016] iterative scheme is then formally (but not computationally)

q (m+1) h (m+1) = q (m) h (m) -J -1 f (m) (4) provided J is invertible. Once the vector c (m+1) q c (m+1) h
T is found as the solution of

J (m) c (m+1) q c (m+1) h = ρ (m) e ρ (m) c , (5)
the new iterates can be computed using (4). Now, the block equations for (5) are, simplifying the notation again,

F c q -A 1 c h = ρ e (6)
and

-A T 1 c q -Ec h = ρ c . (7)
Multiplying ( 6) on the left by

A T 1 F -1 gives A T 1 c q -A T 1 F -1 A 1 c h = A T 1 F -1 ρ e (8)
and adding ( 8) to ( 7) gives

E + A T 1 F -1 A 1 c h = -A T 1 F -1 ρ e + ρ c . (9)
Once c h is determined from this equation, the term c q can be obtained from the following rearrangement of ( 6):

c q = F -1 (A 1 c h + ρ e ). (10)
Equations ( 9) and ( 10) are the PDM counterpart of the GGA method for the DDM problem.

The GGA has been widely used in the solution of the equations for DDM WDSs. For the most part it solves the DDM problems very well provided there are no zero flows when the Hazen-Williams head loss model is used or if the attainable accuracy [START_REF] Dahlquist | Numerical methods[END_REF] for the problem does not inhibit convergence. In the case of zero flows one can apply the regularizations of [START_REF] Elhay | Dealing with zero flows in solving the non-linear equations for water distribution systems[END_REF][START_REF] Piller | Modeling the behavior of a network -Hydraulic analysis and a sampling procedure for estimating the parameters[END_REF] or [START_REF] Carpentier | Water network equilibrium, variational formulation and comparison of numerical algorithms[END_REF]. The problem of low attainable accuracy remains a more significant challenge, probably addressable only with higher precision computing.

The GGA and the Cotree Flows Method (CTM) for the DDM problem are equivalent (?) in the sense that they both solve exactly the same Newton iteration equations for the same WDS. In fact, they produce exactly the same iterates for the same starting values. In a very real sense both methods only have to solve for the heads and flows which satisfy the energy equations because in the GGA the continuity equations are satisfied in every iteration after the first and in the CTM they are satisfied at every iteration. This is because the continuity equations are independent of the heads in the DDM problem. Now, the head loss formulae depend quadratically on the flow rate for the Darcy-Weisbach model and almost quadratically for the Hazen-Williams model and so the region of convergence for the Newton method applied to such a system is very large. This fact explains the very good convergence properties associated with the GGA and the CTM or their variants. But the continuity equations for the PDM problem depend on both the heads and flows. As a consequence, initial values for both flows and heads must be found and these will, in general, not satisfy the PDM continuity or energy equations. Moreover, the PDM continuity equations cannot be satisfied independently of the energy equations as in the DDM case.

It is the experience of the authors and it has been reported elsewhere (see, for example, [START_REF] Siew | Pressure-dependent EPANET extension[END_REF]) that the Newton method defined by ( 9) and ( 10) for the PDM problem exhibits January 19,2016] convergence difficulties. A small example illustrates these difficulties. The network shown in Fig. 1 has the parameters shown in Table 1. The demands shown in Table 1 were magnified by a factor of five (as were the demands in all the networks reported in this paper) to make the problem into a PDM, rather than DDM, problem. The Newton method of ( 9) and ( 10) was applied to this network with each of the four starting value schemes described later in this paper. It failed to converge in 150 iterations after many repetitions of the starting schemes in which there is a pseudo-random element or for one application of the deterministic starting scheme.

The behavior exhibited in this illustrative example is typical of the experience that the authors encountered in applying the simple Newton method of ( 9) and ( 10) to problems of this type. By contrast, a damped version of the Newton method in ( 9) and ( 10) was found to be very reliable and fast, provided suitable step size control measures are used. The [START_REF] Goldstein | Constructive Real Analysis[END_REF] step size selection algorithm, which is discussed later, was found to provide very suitable damping for the Gauss-Newton method for PDM problems. Indeed, all applications of the damped Gauss-Newton scheme with step size selection based on the Goldstein algorithm converged rapidly (usually in about seven iterations but always fewer than 14) for all repetitions of all four starting schemes on this small illustrative network.

DAMPING SCHEMES AND THE EXISTENCE AND UNIQUENESS OF SOLUTIONS

In order to address the issue of damping, four optimization problems are introduced, each of which leads to the system (1). The different formulations are useful because they lead to different metrics for the line search strategies which are used to achieve convergence of the Newton method.

Four equivalent optimization problems

The first optimization problem is couched in terms of the determination of the set of unknown flows. Denote by e j the j th column of an identity matrix of appropriate dimension.

Problem 1.1 Define the content function

C(q) = np i=1 q i 0 δh i (u)du -a T q + n j j=1 -e T j A T 1 q 0 c -1 j (v)dv. ( 11
) Denote U = {q ∈ R np |o ≤ -A T 1 q ≤ d}. Find min q∈U C(q).
The content and the co-content functions, which are co-energy, appear to have been first introduced by [START_REF] Cherry | Cxvii. some general theorems for non-linear systems possessing reactance[END_REF] and [START_REF] Millar | Cxvi. some general theorems for non-linear systems possessing resistance[END_REF] to solve electrical network equations. They proved that solving the network equations for power systems is equivalent to minimizing a co-energy function.

Using the identity [START_REF] Parker | Integrals of inverse functions[END_REF])

y 0 f -1 (v)dv = yf -1 (y) - f -1 (y) f -1 (0) f (w)dw (12)
the last term in (11) may be rewritten to give

min q∈U np i=1 q i 0 δh i (u)du -a T q -c -1 (-A T 1 q) T A T 1 q - n j j=1 c -1 j (-e T j A T 1 q) (hm) j c j (w)dw (13)
[pdmjr3Figs-singl.tex: 18:46, January 19, 2016]

The Lagrangian of this problem is, denoting by α ≥ 0 the Lagrange multiplier vector for the lower bound constraint on A T 1 q and denoting by β ≥ 0 the Lagrange multiplier vector for its upper bound constraint, 14) can be rewritten, showing its dependency on ζ, as

L(q, α, β) = np i=1 q i 0 δhi(u)du -a T q -c -1 (-A T 1 q) T A T 1 q - n j j=1 c -1 j (-e T j A T 1 q) (hm) j cj(w)dw + α T A T 1 q -β T (A T 1 q + d). ( 14
) Denote ζ by ζ = c -1 (-A T 1 q) + β -α. But then (
L(q, ζ, β) = np i=1 q i 0 δh i (u)du -a T q -ζ T A T 1 q -β T d - n j j=1 c -1 j (-e T j A T 1 q) (hm) j c j (w)dw, (15)
whence, provided that the definition of c is extended so that c j = 0 , if h j ≤ (h m ) j , and

c j = d j if h j ≥ (h s ) j , L(q, ζ) = np i=1 q i 0 δh i (u)du -a T q -ζ T A T 1 q - n j j=1 ζ j (hm) j c j (w)dw (16)
and this leads to the equivalent problem of finding min q max ζ L(q, ζ). Importantly, the gradient of L is f (q, ζ) indicating that Eq. ( 1) is a necessary optimality condition and a saddle-point equation. This suggests that ζ and h are identical and we can replace ζ by h to get

Problem 1.2 Find min q max h L(q, h).
The Lagrangian or primal-dual problem is unconstrained.

The fact that δh is a monotonic, C 1 -differentiable function means that it is possible to express q as a function of h, using the first block-equation of (1), as

q(h) = δh -1 (A 1 h + a) (17)
with δh -1 being the function inverse of the head loss model δh. It is possible, by analogy with the approach of [START_REF] Collins | Solving the pipe network analysis problem using optimization techniques[END_REF], to arrive at a Co-content optimization formulation of the PDM problem. Write

Z(h) = L(q(h), h) = np i=1 q i (h) δh -1 (0) δh i (u)du -a T q -h T A T 1 q - n j j=1 h j (hm) j c j (w)dw (18)
and use (12) to get the following formulation for the dual function:

Z(h) = - np i=1 e T i (A 1 h+a) δh -1 (0) δh -1 i (u)du - n j j=1 h j (hm) j c j (w)dw. ( 19
)
The optimization problem associated with this formulation is then

Problem 1.3 Find max h Z(h).
[pdmjr3Figs-singl.tex: 18:46, January 19, 2016]

Solving Problem 1.3 will be referred to as using the Co-Content (CC) approach.

The fourth optimization problem considered here uses the energy and continuity residuals of (3). Denote by W ∈ R (np+n j )×(np+n j ) a diagonal matrix of positive weights and define

θ(q, h) = 1 2 W 1 2 f (q, h) 2 2 = 1 2 f T W f . (20)
The problem considered now is

Problem 1.4 Find min q,h θ(q, h). (21)
Solving Problem 1.4 will be referred to as using the Weighted Least Squares (WLS) approach.

Existence and uniqueness of solutions [START_REF] Piller | Least action principles appropriate to pressure driven models of pipe networks[END_REF] proved that the solutions to Problem 1.1, Problem 1.2 and Problem 1.3 exist provided that the set U is not empty, is closed and that C(q) is continuous and norm-coercive. They proved that there is a unique solution provided that U is convex and C is strictly convex. But, an optimization which reduces the value of the objective function θ(q, h) of ( 20) to zero clearly solves (1). Since the solutions to Problem 1.1, Problem 1.2 and Problem 1.3 are also the solutions to (1) then it follows that the solution to Problem 1.4 always exists, is unique and is the same as the solutions of Problem 1.1, Problem 1.2 and Problem 1.3.

The existence and uniqueness of the DDM solutions are not guaranteed for networks in which unsourced subnetworks are disconnected from their main networks. Then (the equivalent of) U is empty and there is no DDM solution if the subnetwork has any non-zero demands (see [START_REF] Deuerlein | Preprocessing of water distribution systems to assess connectivity and solvability in the presence of flow control devices[END_REF] for more details). However, the PDM problem always has a solution because U is always non-empty.

The existence and uniqueness of solutions to the PDM WDS problems under modest conditions motivates the search for robust and reliable methods to find them. One of the main aims of this paper is to demonstrate the effectiveness of two versions of the damped Gauss-Newton method on the PDM WDS problem: the WLS and the CC approaches. The damping or step-size control algorithms are based on the methods of [START_REF] Goldstein | Constructive Real Analysis[END_REF] and are proven on a set of eight case study networks with between 932 and 19,647 pipes and between 848 and 17,971 nodes. These case study networks were previously used in [START_REF] Simpson | A forest-core partitioning algorithm for speeding up the analysis of water distribution systems[END_REF] and ?.

The method of (4) can also be viewed as the Gauss-Newton method [START_REF] Gratton | Approximate Gauss Newton methods for nonlinear least squares problems[END_REF] for the WLS formulation given in Problem 1.4. This can be seen from the following argument. Recalling the definitions of f (q, h) in (1), denoting

∇ x = ( ∂ ∂x 1 , ∂ ∂x 2 , . . . , ∂ ∂xn ) and noting that J in (2) is symmetric, ∇ T q,h θ(q, h) = J W f (22)
and so the Hessian H ∈ R (np+n j )×(np+n j ) for the objective function θ q (m) , h (m) can be found as

H = ∇ q,h (J W f ) = J W J + Q
where Q involves the second-order terms. The term Q is ignored in the Gauss-Newton method and so the resulting iteration scheme is

q (m+1) h (m+1) = q (m) h (m) -(J W J ) -1 J W f (m) = q (m) h (m) -J -1 f (m) . ( 23
)
[pdmjr3Figs-singl.tex: 18:46, January 19, 2016] and this is just (4). Importantly, the term Q involves the system residuals for least squares problems and if the problem has zero residuals at the solution (as in the present case) then the quadratic convergence of the full-Hessian Newton method obtains in the Gauss-Newton variation [START_REF] Gratton | Approximate Gauss Newton methods for nonlinear least squares problems[END_REF].

DAMPED NEWTON METHOD FOR THE SYSTEM IN EQ. (1)

The damped Newton method for (4) is

q (m+1) h (m+1) = q (m) h (m) -σ (m+1) J -1 f (m) . ( 24
)
for some choice of step-size, σ (m+1) . Thus, when the terms, c

, and c

(m+1) h
, of ( 9) and ( 10) have been found, the new iterate can be computed as

q (m+1) h (m+1) = q (m) h (m) -σ (m+1) c (m+1) q c (m+1) h . ( 25
)
In the next section the step size selection algorithm of [START_REF] Goldstein | Constructive Real Analysis[END_REF] is briefly described. Only the WLS or CC optimization problem objective functions can be used in this approach because Problem 1.1 is a constrained problem and Problem 1.2 is a saddle-point problem.

The Goldstein criteria for step size selection in a minimization problem

Denote by -φ (m) = -φ(q (m) , h (m) ) the descent direction chosen for the m-th step and suppose that the proposed step length is σ (m) . It is assumed that ∇θ (m) φ (m) ≥ 0 since otherwise -φ does not represent a descent direction. If ∇θ (m) φ (m) = 0 then the current point is an extremum or saddle point and no further iteration is justified.

Let 0 < µ 1 ≤ µ 2 < 1 be chosen parameters. Define the (scalar) Goldstein index by

g θ(q (m) , h (m) ), σ (m) = θ q (m) , h (m) -θ q (m+1) , h (m+1)
σ (m) ∇θ (m) φ (m) (26)
where

q (m+1) h (m+1) = q (m) h (m) -σ (m) φ (m) . (27)
If µ 1 ≤ g ≤ µ 2 then the step size σ (m) is accepted. If g > µ 2 then the step length 3σ (m) /2 is proposed.

Otherwise, the step length σ (m) /2 is proposed. Problem 1.1 is a constrained problem and Problem 1.2 is a saddle-point problem. Only the equivalent minimization problems of Problem 1.3 and Problem 1.4 can be used in (26). Indeed, for the case of Problem 1.4, the denominator of ( 26) is

σ (m) ∇θ (m) φ (m) = σ (m) J W f (m) T J -1 f (m) = σ (m) f (m) T W f (m) = 2σ (m) θ (m)
and so (26) simplifies, for this case, to

g θ (m) , σ (m) = θ (m) -θ q (m+1) , h (m+1) 2σ (m) θ (m) . ( 28
)
It is important to note that Goldstein's algorithm is not heuristic. If the conditions of Goldstein's theorem [START_REF] Goldstein | Constructive Real Analysis[END_REF]) are met, then convergence is mathematically guaranteed. However, there are choices that can be made for some of the parameters in the algorithm and different choices of these parameters may affect the speed of convergence. The important point is that the existence and uniqueness of the solution for the WLS and CC formulations is proved and that therefore the conditions of Goldstein's theorem can be met, mathematically guaranteeing convergence for any choice of the parameters within the range specified by the theorem. One advantage of using the WLS formulation of the problem is that the objective function θ of ( 20) is, unlike the corresponding L 1 function in [START_REF] Giustolisi | Pressure-driven demand and leakage simulation for water distribution networks[END_REF] differentiable, something that is required in order to satisfy the conditions of the Goldstein Theorem.

Summary of the algorithm

The algorithm takes input starting values q (0) , h (0) and an objective function, ψ which, in this context, is either the weighted least squares function θ(q, h) of Problem 1.4 or -Z(h), the negative of the co-content function of Problem 1.3.

(a) Compute φ (m) as the solution of J φ (m) = -f (m) and compute ∇ψ (m) (b) If ∇ψ (m) φ (m) = 0 then no descent possible, no further iteration is justified.

(i) set σ (m) = 0 and (ii) set q (m+1) = q (m) , and h (m+1) = h (m) .

(iii) Exit.

(c) If ∇ψ (m) φ (m) > 0 then, set σ (m) = 1, choose 0 < µ 1 ≤ µ 2 < 1 and proceed as follows:

(i) Compute q (m+1) h (m+1) = q (m) h (m) -σ (m) φ (m)
and the Goldstein index

g ψ (m) , σ (m) = ψ (m) -ψ q (m+1) , h (m+1)
σ (m) ∇ψ (m) φ (m) (ii) If µ 1 ≤ g ψ (m) , σ (m) ≤ µ 2 then set h (m+1) = h (m+1)
and, in the case of the WLS formulation of the problem set q (m+1) = q (m+1) , and then increment m and go to step (a).

(iii) ElseIf g ψ (m) , σ (m) > µ 2 then increase step length: set σ (m) = 3σ (m) /2 and go to step (c)(i) (iv) Else decrease step length: set σ (m) = σ (m) /2 and go to step (c)(i) m) represents an ascent, not descent, direction and this indicates an error condition. Exit.

(d) If ∇ψ (m) φ (m) < 0 then -φ (
It is worth noting that when using the WLS formulation of the problem, both the heads and the flows are updated at each step whereas in the CC formulation only the heads are updated. Before presenting results which illustrate the effectiveness of the methods, some preliminary issues are addressed.

MODELLING CHOICES

[pdmjr3Figs-singl.tex: 18:46, January 19, 2016]

In moving from a DDM to a PDM there are two important model choices to be made: (i) the consumption function model and (ii) the starting values to be used in the iteration. Some particular choices for these models and the consequences of their use are discussed in the sections following.

The consumption function

The consumption function describes what is sometimes called the nodal hydraulic availability or nodal delivery in a system (for a useful review of four consumption function models in the context of reliability assessment and analysis see [START_REF] Jun | Iterative methodology of pressure-dependent demand based on EPANET for pressure-deficient water distribution analysis[END_REF]). The flow, q, at an aperture has usually been modelled by a relationship in which the flow is proportional to a power n of the pressure head h, q ∝ h n and where n has been variously estimated [START_REF] Van Zyl | The effect of pressure on leakage in water distribution systems[END_REF][START_REF] Cheung | Extension of EPANET for pressure driven demand modeling in water distribution system[END_REF] to lie in the interval n ∈ [0.5, 2.79]. [START_REF] Tanyimboh | A new nodal outflow function for water distribution networks[END_REF] proposed a consumption function whose form around h j = (h m ) j more closely resembles the choice of exponent n = 2:

c T (h j ) = d j e α+βh j 1 + e α+βh j
, all h j , and where the parameters α and β can be derived empirically or, in the absence of empirical data, with a formula provided by the authors. Yet another variation, which uses sinusoidal functions, was proposed by [START_REF] Tucciarelli | Leak analysis in pipeline systems by means of optimal valve regulation[END_REF]. [START_REF] Wagner | Water distribution reliability: simulation methods[END_REF] proposed a consumption function whose form is based on the exponent choice n = 0.5. Let h j denote the head at node j. Denote also

z(h j ) = h j -(h m ) j (h s ) j -(h m ) j . ( 29
)
The Wagner consumption function is defined by

c W (h j ) =    0 if z(h j ) ≤ 0 d j z(h j ) if 0 < z(h j ) < 1 d j if z(h j ) ≥ 1 (30)
where d j denotes the demand at the j-th node. The Wagner consumption function has a discontinuous derivative at h j = (h m ) j and its value at h j = (h s ) j is less than d j and these properties sometimes have undesirable effects on the convergence behaviour of the iterative methods (see e.g. [START_REF] Ackley | Head driven analysis of water distribution systems[END_REF], [START_REF] Giustolisi | Water distribution network pressure-driven analysis using the enhanced global gradient algorithm (egga)[END_REF] and [START_REF] Muranho | Pressure-dependent demand and leakage modelling with an EPANET extension WaterNetGen[END_REF]). Because of these effects, [START_REF] Piller | Least action principles appropriate to pressure driven models of pipe networks[END_REF] proposed regularizing the function by smoothing it with a cubic interpolating polynomial which matches function and derivative values either side of the points h j = (h m ) j and h j = (h s ) j . Thus, the Regularized Wagner consumption function, denoted here by c R (h), is continuous and has a continuous first derivative.

The choice of n = 0.5 in the design of the Wagner function is based on a model that applies to a single, circular aperture and it describes the instantaneous flow for given pressure heads. The nodal demands in a network model that is not an all-pipes model are frequently, in practice, derived by measuring total water usage for a group of 50-100 houses over a period of some months and then calculating an average daily use for the whole collection of houses represented by that single node. Clearly, the delivery at empirically derived demands such as these are not faithfully modelled by c W (h).

Even where an all-pipes model is used, a formula based on the flow at a single outlet is unlikely to faithfully model water consumption in a setting where showers, toilets, irrigation systems and taps are all used.

A C 1 cubic consumption function, c C (h j ), was studied in the context of reliability analysis in [START_REF] Fujiwara | Reliability assessment of water supply systems with storage and distribution networks[END_REF], where it was first proposed, and in [START_REF] Fujiwara | Reliability analysis of water distribution networks in consideration of equity, redistribution, and pressure-dependent demand[END_REF]). It bears [pdmjr3Figs-singl.tex: 18:46, January 19, 2016] in (b) was found to provide starting values that lead to convergence but not as often as the schemes in (a)(i) and (a)(ii).

ILLUSTRATION OF THE WLS AND CC METHODS

In what follows, the results of applying the WLS and CC methods to eight case study networks are reported in order to illustrate the viability of the methods on a variety of quite different, and challenging, networks. Firstly, the case study networks are described and some implementation details are given. Secondly, the convergence behaviours of the two methods are described and a comparison is made of how that behaviour is affected by which of the consumption functions, c C (h j ) and c R (h j ), is used. Thirdly, the differences between the solutions which result from using c C (h j ) and c R (h j ) are reported.

Implementation and the details of the case studies

All the calculations reported in this paper were done using codes specially written for Matlab 2012b and 2013a (The Mathworks 2012, 2013) and which exploit the sparse matrix arithmetic facilities available in that package. Matlab implements arithmetic that conforms to the IEEE Double Precision Standard and so machine epsilon for all these calculations was 2.2 × 10 -16 .

Columns 2, 3 and 4 of Table 2 show the numbers of pipes, n p , the numbers of nodes, n j , and the numbers of sources, n f , for the eight case study networks used for testing. All the networks use the Darcy-Weisbach head loss model. These networks, apart from some necessary changes, are those used previously in [START_REF] Simpson | A forest-core partitioning algorithm for speeding up the analysis of water distribution systems[END_REF] and ?. Four of the networks used in this paper are available as supplemental data. In all cases the demands of the network were magnified by multiplying them by a factor of five to ensure that the problem was actually a PDM problem and not a DDM problem.

In all tests reported here the minimum pressure head and service pressure head were set to h m = 0 m and h s = 20 m, respectively. The iteration stopping test was

h (m+1) -h (m) ∞ h (m+1) ∞ ≤ , and q (m+1) -q (m) ∞ q (m+1)
∞ ≤ with = 10 -6 . This tolerance was used here to confirm the quadratic convergence of Newton's method even though such a small tolerance is unlikely to be required in practical applications.

The starting scheme described in Section (a)(ii) was used for all the tests and the same seed was used to start the pseudo-random number generators for all runs. The Goldstein index limits were set, as a result of testing, to µ 1 = 1/10 and µ 2 = 1 -µ 1 in all the testing reported here. The Goldstein line search for suitable damping requires the calculation, at each iteration, of the Goldstein index (26). For the WLS scheme the expression in (26) reduces to (28) and so each subiteration during the line search requires one calculation of the expression θ q (m+1) , h

. This involves recomputing the right hand side of ( 27) with the proposed value of σ (m) and then forming θ q (m+1) , h

, a relatively fast computation. Computing the Goldstein index for the CC line search requires the computation of the expression in (18), part of which involves evaluating the inverse head loss function δh -1 i of (17) to get q i (h). This inversion is a simple matter if the head loss is modeled by the Hazen-Williams formula but it is more challenging when the head loss is modeled by the Darcy-Weisbach formula which takes quite different forms for laminar, transitional and turbulent Reynolds numbers. Given the difficulty (or perhaps the impossibility) of finding closed-form expressions for the inverse function in that case, this inversion was performed using the Matlab function fsolve in the calculations for this report.

[pdmjr3Figs-singl.tex: 18:46, January 19, 2016]

The integrals were evaluated using the Matlab function integral. The impact of these differences between the WLS and CC line search, or subiteration, calculations is discussed later.

The residuals in the objective function for Problem 1.4 should be weighted to account for significant differences in scale of the heads and flows data. Denote the inverse, diagonal, weighting matrix for the energy residuals by M and the inverse, diagonal, weighting matrix for the continuity residuals by N . Then,W = diag M -1 , N -1 . In this study the weights used were based on demands and fixedhead node elevations: the energy residuals are each divided by the maximum head among the fixedhead nodes and the continuity residuals are weighted by dividing all components by the maximum demand (it is assumed that not all demands are zero). Thus, for this case M = (max h f ) 2 I and N = (max d i ) 2 I. This weighting scheme proved satisfactory but least squares schemes in which the residuals were unweighted frequently led to convergence difficulties.

Convergence behaviour

Columns 5-12 of Table 2 show the numbers of iterations and subiterations, or line search steps, that were required to solve the eight case study networks by both the WLS and CC methods and for the two consumption functions c C (h j ) and c R (h j ). Both the WLS and CC schemes converged in quite modest numbers of iterations with both consumption functions for all the networks. The WLS scheme required many fewer iterations than the CC scheme and, in all but one case, required many fewer subiterations than the CC scheme. The main iterations in both cases require comparable computation but, as was pointed out earlier, there is some difference between the two methods in the computation required for subiterations. On one hand, each subiteration of the WLS scheme requires one evaluation of objective function θ of (20), a simple and rapid calculation. On the other hand, each subiteration of the CC scheme requires one evaluation of objective function Z of (18). The second integral in ( 18) is simple to compute explicitly for both of the consumption functions c C (h j ) or c R (h j ). But the first integral in (18) involves the inversion of the function δ(h j ) and, while this inversion for the Hazen-Williams head loss model can be written in closed form, it requires significant computation if the head loss is modelled by the Darcy-Weisbach formula.

The authors believe that WLS approach provides the preferred choice: it is easier to implement and, although no carefully designed timings tests have been conducted to compare the WLS and CC methods, it appears to be faster than the CC method. The difficulties associated with the CC line search when head loss is modelled by the Darcy-Weisbach formula make the CC method less attractive. In any case, both have been demonstrated to converge rapidly on a wide range of network types.

Consumption function effects

The choice of consumption function can, in some cases, have a noticeable effect on the solution heads and flows of a PDM problem. Nodes in the network which have positive demand will be referred to as demand nodes. Recall that demand nodes in a network for which the PDM solution has zero delivery (c(h j ) = 0) are said to be in failure mode, demand nodes for which the delivery falls between the minimum and the service level (0 < c(h j ) < d) are said to provide partial delivery and demand nodes which deliver the full demand (c(h j ) = d) are said to give full delivery. In what follows the numbers of demand nodes in these three categories are reported for the PDM solutions of the eight case study networks. The number of nodes for which the solution has negative pressures (i.e. for which the delivery is zero) is also reported.

Node counts for failure, partial delivery and full delivery

Table 3 compares various aspects of the solutions for the two consumption functions c R (h j ) and c C (h j ). Columns 2 and 3 show the total deliveries as percentages of the total initial demands. Column 4 shows the numbers of demand nodes. Columns 5 and 6 show the numbers of demand nodes in failure mode, Columns 7 and 8 show the numbers of demand nodes in partial delivery mode and Columns 9 and 10 show the numbers of demand nodes in full delivery mode. The last two columns show the numbers of nodes in the solutions for which the pressure is below zero. Although in most cases the numbers of demand nodes in the different modes are similar, there are some quite marked differences. Networks N 4 , N 5 and N 7 show quite large numbers of nodes in failure mode when c R (h j ) is used but not when c C (h j ) is used.

Head differences

Frequency distributions of the differences between the heads and flows of the two solutions obtained using the two consumption functions with each network were produced in order to to better understand the effect that the choice of consumption function can have on the solutions obtained. Fig. 3 shows the frequency distributions of the differences in the heads (m) between the solutions for the Regularized Wagner consumption function, c R (h j ), and the cubic consumption function c C (h j ) for Network N 1 . Although most heads there are very similar, some 100 of the 848 heads in that case differ by as much as 2 m. The variation in differences between solution heads for the two consumption functions across the other case study networks is quite marked. Fig. 4 shows the corresponding frequency distributions for the flows (L/s) and shows greater agreement between the two solutions than for the heads.

Another characterization of the differences between the solution heads and flows for the two consumption functions c R (h j ) and c C (h j ) can be seen in Table 4. There, Columns 2 and 3 show the intervals containing most of the differences of the heads and flows, respectively. Thus, for N 4 almost all the solution heads differ by more than 3 m but less than 5 m and the solution flows for N 8 differ by no more than 0.5 L/s. Columns 4 and 5 show, respectively, the means of head and absolute flow differences. The scale of the differences between the solution heads for the two consumption functions c R (h j ) and c C (h j ) suggests that more research is necessary to find and calibrate appropriate models of consumption at demand nodes.

CONCLUSIONS

The Newton method PDM counterpart of the GGA for DDM problems is shown, by a small example, to exhibit failure to converge if no damping is used. This behaviour has been reported elsewhere. It has been shown that a new (fourth) formulation of the PDM problem, the WLS optimization formulation, is equivalent to three known (equivalent) PDM problem formulations. The conditions for the existence and uniqueness of the solution to the WLS formulation follow and two of the four equivalent optimization problems, the CC and WLS versions, are used as the bases for Gauss-Newton methods with Goldstein step selection. The damped method is proved, on a challenging set of eight case study networks, to have convergence behaviour that mirrors that of the GGA on DDM problems. The line search scheme based on the WLS optimization problem is shown to be significantly more economical than that based on the CC optimization. Thus, the PDM counterpart to the GGA for DDM problems is seen to be solvable robustly and rapidly provided the recommended modifications to the Newton method are employed.

The cubic consumption function, c C (h), of [START_REF] Fujiwara | Reliability assessment of water supply systems with storage and distribution networks[END_REF] is compared with the Regularized Wagner function of [START_REF] Piller | Least action principles appropriate to pressure driven models of pipe networks[END_REF], c R (h). In particular, the number of iterations required for solution and the differences in heads and flows between solutions obtained were compared. The steady-state solution heads for c C (h) differed from those for c W (h) by as much as 5 m for some nodes. The reasons for these differences were not investigated and more work is needed in order to better understand the effects that the consumption function choice has on the solutions.

Four starting value schemes for the heads in the system (unnecessary to initiate the DDM problem but necessary for the PDM problem) were proposed and compared. The two which use equal flow velocities and pseudo-random heads were found to be very effective and another, based on using the invertibility of the matrix A T 1 A 1 was found to be unreliable. The scheme based on the DDM solution of the problem was found to be less reliable than the two best schemes but sometimes effective. The WLS PDM solution method reliably finds the solution in roughly the same number of iterations as are required to find a solution to the corresponding DDM problem for the same network. Given the small number of iterations required by the new method, it would be hard to recommend a starting scheme in which the number of iterations to find the starting values is the about same as the number of iterations to find the PDM solution.

A residual weighting scheme based on maximum fixed-head elevation and maximum nodal demand was proposed and the authors' experiments suggest that the proposed scheme is quite suitable and that unweighted schemes can present convergence difficulties. Furthermore, the wide range of delivery fractions and PDM node fractions together with the small number of iterations required to solve these challenging case study networks of quite different scales suggests that the methods proposed in this paper are likely to be suitable for a wide range of PDM problems.

The robust solution algorithm introduced in this paper is able to deal with, amongst other conditions, insufficient pressures and excessive demands. Networks N 1 , N 2 , N 5 and N 6 were derived from real world networks by removing pumps and control devices. The extension of this work to systems which have pumps and control devices would be a useful contribution to the field as would the investigation of this technique applied to extended period simulations and rigid water column modeling. There is also a great need for improved mathematical methods that successfully deal with ill-posed systems and other situations where existing modelling techniques reach the limits of their theoretical bases. Thus, future work could aim to develop hydraulic models suitable for extreme operational conditions (which can have a significant impact on the hydraulic performance of control devices and pumping stations) or even extreme event situations like natural disasters, terrorist attacks or electrical power blackouts. The stable and robust calculation of WDS hydraulics in such anomalous situations is a basic requirement for all model-based decision systems. Existing simulation techniques cannot handle these critical events adequately and often fail because of the lack of convergence.

Indeed, in the case where the hydraulic simulations run online, the robustness of the solver is particularly important: the operational data are transferred from the supervisory control and data acquisition system which automatically updates the states of valves, pumps, etc. and that data is fed directly to an online solver. Network operations and catastrophic events sometimes cause parts of a network to suffer from insufficient pressure or sometimes segment a network into components which have inadequate connections to sources or perhaps have no connection at all to a source. In such a case, the resulting system can be underdetermined and existing solvers often fail to converge, converge to the wrong solutions or even worse, cease executing. This is not acceptable for practical online simulation. Developing techniques to handle such conditions in PDM systems remains a challenge for researchers in this field.

SUPPLEMENTAL DATA

The data for case study networks N1, N3, N4 and N7, which are modifications of networks in the public domain, are available online in the ASCE Library (www.ascelibrary.org). The other four networks considered in this paper are not available because of security concerns. [pdmjr3Figs-singl.tex: 18:46, January 19, 2016]

[

  pdmjr3Figs-singl.tex: 18:46, January19, 2016]

[

  pdmjr3Figs-singl.tex: 18:46, January19, 2016]

Figure 2 :

 2 Figure 2: A family of curves showing flow q as a function proportional to h n for a range of exponents n together with the cubic consumption function of Fujiwara & Ganesharajah (1993), c C (h), The exponential consumption function of Tanyimboh & Templeman (2004), c T (h), and the the Regularized Wagner function of Piller (1995), c R (h), (which is based on a value of n = 0.5).

Figure 3 :Figure 4 :

 34 Figure 3: Frequency distributions of the differences in the heads (m) between the solutions for the Regularized Wagner consumption function, c R (h), and the cubic consumption function c C (h) for Network N 1 .

Table 1 :

 1 Pipe and node data for the network shown in Fig.1. The network has a single reservoir, Node 1, with an water surface elevation of 100 m. The demands that are shown above were magnified by a factor five to cause the problem to be a PDM, rather than DDM, problem. τ i τ si τ i τ si τ i τ si τ i τ si

	TABLES
				Pipes			Nodes
	ID From To L(m) D(mm)	(mm) Elev (m) d (L/s)
	1	1	2 1000	100	0.3	-	-
	2	1	4	400	300	0.3	10	50
	3	2	3	400	200	0.3	10	30
	4	2	5	100	300	0.3	19	20
	5	3	6	500	200	0.3	10	30
	6	4	5	700	300	0.3	5	0
	7	4	7	700	200	0.3	9	80
	8	5	6	400	300	0.3	5	90
	9	5	8	400	250	0.3	0	90
	10	6	9	100	300	0.3	-	-
	11	7	8	900	300	0.3	-	-
	12	8	9	500	300	0.3	-	-

[pdmjr3Figs-singl.tex: 18:46, January 19, 2016] f

Table 2 :

 2 Number of pipes, n p , nodes, n j , sources, n f , iterations, τ i , and subiterations, τ si , required to solve the eight case study networks by WLS and CC schemes for the two consumption functions c C (h) and c R (h).

		% Delivery	d > 0	c(h) = 0	0 < c(h) < d	c(h) = d	p < 0
	ID c C (h) c R (h) c N 1 86.9 89.0 474	11	13	135	142	328	319	50	54
	N 2	52.5	65.7	661	34	38	503	507	124	116	42	47
	N 3	92.1	93.9	1770	34	34	221	227	1515	1509	34	34
	N 4	26.8	27.4	1609	21	347	1521	1211	67	51	21	380
	N 5	49.2	51.3	1241	35	195	1168	1023	38	23	80	421
	N 6	68.6	70.6	3173	37	48	2683	2733	453	392	122	145
	N 7	56.5	59.6 10552	74	457	9505	9313	973	782	85	534
	N 8	97.2	97.7 15332	0	0	3119	3206 12213 12126	1	1

C (h) c R (h) c C (h) c R (h) c C (h) c R (h) c C (h) c R (h)

Table 3 :

 3 Comparison of the deliveries, numbers of demand nodes, nodes in failure mode, partial delivery mode and full delivery mode, and nodes with negative pressure for the cubic consumption function, c C (h), and the Regularized Wagner consumption function, c R (h).

[pdmjr3Figs-singl.tex: 18:46, January

19, 2016]

Visiting Research Fellow, School of Computer Science, University of Adelaide, South Australia, 5005, sylvan.elhay@adelaide.edu.au.

Senior Research Scientist, Irstea, Water Department, Bordeaux Regional Centre, Cestas F-33612, France.

Senior Researcher, 3S Consult GmbH, Karlsruhe, Germany & Adjunct Senior Lecturer, School of Civil, Environmental and Mining Engineering, University of Adelaide, South Australia, 5005.

 4 Professor, School of Civil, Environmental and Mining Engineering, University of Adelaide, South

Australia, 5005. [pdmjr3Figs-singl.tex: 18:46, January 19, 2016]

Author-produced version of the article published in J. Water Res. Plan. & Man., 2016, 142(2), p. 04015047-1-04015047-12 The original publication is available at http://ascelibrary.org/doi/10.1061/%28ASCE%29WR.1943-5452.0000578 ¡http://ascelibrary.org/doi/10.1061/(ASCE) WR.1943-5452.0000578¿ doi : 10.1061/(ASCE) WR.1943-5452.0000578 some resemblance to c T (h j ) but, unlike c T (h j ), attains the values 0 and d j at the left and right endpoints of the interval and has zero derivatives at those two endpoints. This function is well integrated into a PDM solver and it was used in this investigation. Its form and properties are now briefly reviewed and its effect is examined in what follows.

Denote r(t) = t 2 (3 -2t), t the independent variable. The cubic consumption function, c C (h j ), is defined by

The consumption functions c R (h j ), c T (h j ) and c C (h j ) are shown in Fig. 2 along with a family of curves which show consumption curves proportional to h n with various values of n ∈ [0.5, 2.79].

A natural question concerns what effect, if any, choosing two different consumption functions would have on the the solution process and the solutions. In particular, would one of the consumption functions require more computation than the other for the same problem? And would the solutions so obtained differ by much between the two cases? These questions are addressed in a later section by comparing the results of using the consumption functions c R (h j ) and c C (h j ). Some investigations by other authors have used different consumption function models for different nodes. It is assumed in this investigation that all nodes in the WDS have the same consumption function in order that the comparison of the effects of the two consumption functions considered are made more apparent.

Starting values for the heads

The PDM problem requires values for both the initial flows, q (0) , and heads, h (0) . The following schemes were investigated.

(a) All flow velocities equal, pseudo-random heads: q (0) consistent with a velocity of 0.3048 m/s (= 1 ft/s) in each pipe and the heads chosen to be

r, where r is the sampled value of a pseudo-random variable uniformly distributed in [0, 1], or (ii) h j = r where r is the sampled value of pseudo-random variable uniformly distributed in [h m , h s ], h m the minimum pressure head and h s the service pressure head.

(b) Initial flows and heads from the GGA solution of the DDM problem:

Here the DDM solution is used as the starting point for the PDM problem. Any negative head in the solution is replaced by either the formula in (a)(i) or that in (a)(ii).

Another scheme, in which the initial flows q (0) are set to match a velocity of 0.3048 m/s (=1 ft/s) and h (0) is found as the solution to A T 1 A 1 h (0) = A T 1 (Gq (0) -a) (which is easily derived from the first block equation of (1)), was also trialed. The matrix A T 1 A 1 is guaranteed, by the full rank of A 1 , to be invertible. This scheme proved to be unreliable. In fact, the schemes in (a)(i) and (a)(ii) were found, by the authors, to provide the most reliably successful starting values. The scheme described Figure 1: The small illustrative network described in Table 1 and used to demonstrate the failure of the undamped Newton method to converge.

[pdmjr3Figs-singl.tex: 18:46, January 19, 2016]


Peter Benner
email: benner@mpi-magdeburg.mpg.de

Jens Saak
email: saak@mpi-magdeburg.mpg.de

Martin Stoll
email: stoll@mpi-magdeburg.mpg.de

Heiko K Weichelt
email: weichelt@mpi-magdeburg.mpg.de

Efficient solvers for large-scale saddle point systems arising in feedback stabilization of multi-field flow problems

Keywords: Coupled flow control, large-scale saddle point systems, preconditioned GMRES, least-squares commutator approach, algebraic multigrid, Chebyshev-Semi-Iteration

. Combining well known approximation methods like a least-squares commutator approach for the Navier-Stokes Schur complement, an algebraic multigrid method, and a Chebyshev-Semi-Iteration, an efficient preconditioner is derived and tested for different parameter sets by using a simplified reactor model that describes the spread concentration of a reactive species forced by an incompressible velocity field.

Introduction

In this paper we investigate the solution of large-scale saddle point systems arising in control problems for coupled partial differential equations (PDEs). The starting points are recent publications concerning the boundary feedback stabilization of non-coupled flows like the linear Stokes flow in [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF] and the non-linear Navier-Stokes flow in [START_REF] Bänsch | Riccati-based boundary feedback stabilization of incompressible Navier-Stokes flow[END_REF]. The analytic approach for this feedback stabilization is given by Raymond in, e.g., [START_REF] Raymond | Feedback boundary stabilization of the two-dimensional Navier-Stokes equations[END_REF].

Using the projection idea proposed by Heinkenschloss et al. [START_REF] Heinkenschloss | Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations[END_REF], Benner et al. [START_REF] Bänsch | Riccati-based boundary feedback stabilization of incompressible Navier-Stokes flow[END_REF][START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF] show that the solution of certain saddle point systems is the key ingredient to ensure that the numerical solution lies on the correct solution manifold, i.e., the space of discretely divergence-free velocity fields, without performing an explicit projection.

Applying these ideas to a coupled flow problem, namely the Navier-Stokes equations combined with a diffusion-convection equation, leads to saddle point systems with a more complicated block structure [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF]. Solving these systems efficiently requires the use of appropriate preconditioners. This paper investigates an efficient iterative solution strategy via the use of preconditioned Krylov subspace methods based on the framework derived in [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF]. Here we consider the full feedback system for the coupled multi-field flow problem, while in [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF], only the linear Stokes case was treated without coupling to another field equation. Moreover, this paper complements [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF] in the sense that there, we have focused on presenting results on the convergence of the Newton-ADI method for solving the algebraic Riccati equation determining the stabilizing feedback control for the coupled system, where the saddle point problems in the innermost step of the Newton-ADI iteration were solved by sparse direct methods, while here, we study preconditioned iterative solvers for this step.

The remainder of this paper is organized as follows. Section 2 briefly recalls the feedback stabilization approach for multi-field flow problems from [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF] that leads to large-scale saddle point systems. Afterwards, we discuss properties of these saddle point systems to derive a suitable preconditioner in Section 3. Section 4 shows numerical results before we conclude the paper and give a short outlook to further investigations in Section 5.

Derivation of saddle point systems

The derivation of the block structured saddle point systems in [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF] starts with the linearized coupled flow problem defined for t ∈ [0, ∞) and x ∈ Ω ⊂ R 2 . The linearized Navier-Stokes equations that describe, up to first order, the difference between actual and desired velocity and pressure are given as

∂ ∂t z - 1 Re ∆z + (w • ∇)z + (z • ∇)w + ∇p = f l , div z = 0, on [0, ∞) × Ω. (1)
They are then coupled via the velocity field z(t, x) with the linearized diffusionconvection equation

∂ ∂t c z - 1 ReSc ∆c z + (w • ∇)c z + (z • ∇)c w = 0, on [0, ∞) × Ω (2)
that describes the concentration of a reactive species denoted by c z (t, x). The stationary linearization points w(x) for the velocity and c w (x) for the concentration are assumed to be given. The equations are scaled with the Reynolds number Re and the Schmidt number Sc. Using the mixed Taylor-Hood finite elements [START_REF] Hood | Navier-Stokes equations using mixed interpolation[END_REF] for the velocity and pressure in Eqs. (1) as well as linear ansatz functions for the concentration in Eq. ( 2), we end up with a system of discrete differential-algebraic equations (DAE) that can be written as the control system:

  M z 0 0 0 0 0 0 0 M c   d dt   z p c   =   A z G 0 G T 0 0 -R 0 A c     z p c   +   B z 0 0   u, (3a)
y = 0 0 C c   z p c   (3b)
with the first block row for velocity (of dimension n z ), the second row for pressure (of dimension n p ), and the third row for concentration (of dimension n c ) [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF].

The matrix pencil

  A z G 0 G T 0 0 -R 0 A c   A ,   M z 0 0 0 0 0 0 0 M c   M
is of dimension n × n with n = n z + n c + n p and has 2n p infinite eigenvalues [START_REF] Cliffe | Eigenvalues of block matrices arising from problems in fluid mechanics[END_REF].

In [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF], a linear-quadratic regulator (LQR) approach is applied to system (3) for determining the stabilizing control function u. The solution of this LQR problem is a linear feedback control u(t) = K(z(t), p(t), c(t)), determined via the solution of an algebraic Riccati equation (ARE) defined on the subspace of discretely divergence-free vector fields. The resulting ARE is then solved using a Newton-ADI algorithm. This method yields a threefold nested iteration. In the innermost loop, saddle point systems of the form

  A T z + q i M z G -R T G T 0 0 0 0 A T c + q i M c   =A T +qiM=:Fi   Λ z Λ p Λ c   Λ =   Ỹz 0 Ỹc   . Y (4)
have to be solved for certain ADI shifts q i ∈ C -and a block right hand side Y.

The whole nested iteration is given in [2, Algorithm 1] and is omitted here due to space constraints.

Preconditioned iterative solvers for block structured saddle point systems

The use of direct solvers in ( 4) is only suitable for moderate problem sizes and two-dimensional problems. Although iterative methods can handle much larger systems, their performance will deteriorate if the mesh-size decreases. To avoid this, a suitable preconditioner P i ∈ C n×n is introduced such that [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF][START_REF] Saad | Iterative methods for sparse linear systems[END_REF]). Before we derive a suitable preconditioner P i we need to describe the properties of the saddle point system and their influence on the chosen preconditioner.

P -1 i F i Λ = P -1 i Y is solved instead of (4) (see

Properties

The matrices M z , M c are symmetric and positive definite, G, R are of full rank, and the ADI shift q i ∈ C -is contained in the convex hull of the finite spectrum of (A, M). The shifted system matrix F i is indefinite ∀q i ∈ C -. Due to the different q i , the matrix F i changes in each ADI step and, therefore, the preconditioner has to be adapted in each ADI step as well. Nevertheless, for the remainder of this section we assume a fixed ADI shift q i = q to omit the index i if it is obvious.

Derivation of block preconditioner

Adapting the ideas from [3, Section 3.2] we consider

F =   F z G -R T G T 0 0 0 0 F c   = F N SE -RT 0 F c with F z := A T z + qM z , F c := A T c + qM c , R := R 0 , (5)
and F N SE as the saddle point system for the non-coupled Navier-Stokes flow as it is used in [START_REF] Bänsch | Riccati-based boundary feedback stabilization of incompressible Navier-Stokes flow[END_REF]. Using the preconditioner P N SE from [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF], we define a block preconditioner for the use with GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] to solve with the block structured saddle point system (5) as follows:

P = P N SE -RT 0 P c =   P z 0 -R T G T -P SC 0 0 0 P c   ⇒ P -1 =   P -1 z 0 P -1 z R T P -1 c P -1 SC G T P -1 z -P -1 SC P -1 SC G T P -1 z R T P -1 c 0 0 P -1 c   .
In contrast to the preconditioner derived in [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF], we cannot achieve a block lower triangular matrix due to the coupling matrix R. Applying P -1 to F yields

P -1 F = P -1 z Fz P -1 z G -P -1 z R T +P -1 z R T P -1 c Fc P -1 SC G T P -1 z Fz-P -1 SC G T P -1 SC G T P -1 z G -P -1 SC G T P -1 z R T +P -1 SC G T P -1 z R T P -1 c Fc 0 0 P -1 c Fc (6)
If one assumes P z = F z , P c = F c , and P SC = G T F -1 z G as ideal approximations in [START_REF] Elman | Block preconditioners based on approximate commutators[END_REF], this leads to

  I z F -1 z G -F -1 z R T + F -1 z R T P -1 SC G T -P -1 SC G T P -1 SC G T F -1 z G -P -1 SC G T F -1 z R T + P -1 SC G T F -1 z R T 0 0 I c   =   I z * 0 0 I p 0 0 0 I c  
and our iterative method would converge within one step. The goal is to find good approximations for P z , P c , and P SC that can be evaluated fast and still cluster the eigenvalues in a suitable way such that our iterative solver shows fast convergence [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF]. Instead of calculating the inverse P -1 to apply the preconditioner P, we consider the solution of a linear system

  P z 0 -R T G T -P SC 0 0 0 P c     x z x p x c   =   b z b p b c   (7)
that can be solved in three steps:

Step I:

x c = P -1 c b c , (8a)
Step II:

x z = P -1 z (R T x c + b z ), (8b)
Step III:

x p = P -1 SC (G T x z -b p ). ( 8c
)
In conclusion, the coupling matrix R only leads to a matrix-vector multiplication.

In steps I and II, one needs to solve with the shifted velocity and concentration system matrices as defined in [START_REF] Cliffe | Eigenvalues of block matrices arising from problems in fluid mechanics[END_REF]. For both steps, an algebraic multigrid (AMG) method can be used as it is described below. But first, we discuss the more challenging step III that is handled as follows.

Approximation methods

Schur Complement Approximation. P SC is an approximation of the Navier-Stokes Schur complement SC := G T F -1 z G ∈ R np×np . Unfortunately, the matrix SC would be a dense matrix that includes the inverse of F z . To avoid the use of this matrix, we follow the approach in [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF][START_REF] Stoll | All-at-once solution of time-dependent Stokes control[END_REF] and use a slightly modified variant of the least squares commutator approach as it is described in [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF]Section 8.2]. Namely, we consider the shifted Oseen operator in the velocity space

F z = - 1 Re ∇ 2 + w • ∇ + qI.
Note that it is common practice to omit the reaction term (z • ∇)w that appears in the linearized Navier-Stokes equations to derive preconditioners [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF]Section 8]. Similar to [7, Section 8.2] and [START_REF] Elman | Block preconditioners based on approximate commutators[END_REF], we suppose that there exists an analogous operator on the pressure space defined as

F p = (- 1 Re ∇ 2 + w • ∇ + qI) p .
The least squares commutator of the shifted Oseen operator with the gradient operator is defined as

E = (F)∇ -∇(F p )
and is supposed to become small in some sense [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF]. Using the discrete versions of the operators, we end up with

E = (M -1 z F z )M -1 z G -M -1 z G(M -1 p F p )
with M p the mass matrix and F p = A T p + qM p the shifted system matrix, both defined on the pressure space. Premultiplying this by G T F -1 z M -1 z and postmultiplying by F -1 p M p yields

G T M -1 z GF -1 p M p ≈ G T F -1 z G = SC, [3]
. The large and dense matrix G T M -1 z G cannot be used explicitly, but it is shown in [7, Section 5.5.1] that this matrix is spectrally equivalent to the Laplacian S p defined on the pressure space for an inf-sup stable discretization and an inflow-outflow problem [7, Section 8.2], as it is considered in this paper. Finally, we obtain

P SC ≈ S p F -1 p M p ⇒ P -1 SC ≈ M -1 p F p S -1 p .
In [START_REF] Benzi | Modified augmented Lagrangian preconditioners for the incompressible Navier-Stokes equations[END_REF] the authors use a similar approach for the Navier-Stokes equations. In summary, the application of P -1 SC requires to solve with S p (step IIIa), multiply with F p (step IIIb), and solve with M p (step IIIc). The step IIIa can be solved with an AMG method, similar to the steps I and II.

Algebraic Multigrid. As it is described above, the steps I (8a), II (8b) and IIIa are solved using an AMG method [START_REF] Ruge | Algebraic multigrid (amg)[END_REF]. Due to the possibly complex ADI shifts q in (8a) and (8b), we use the AGMG package developed by the group of Y. Notay [START_REF] Napov | An algebraic multigrid method with guaranteed convergence rate[END_REF][START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF][START_REF] Notay | Aggregation-based algebraic multigrid for convection-diffusion equations[END_REF]. In all three cases we use the MATLAB ® -based implementation to solve systems of the form

F x = b
with a sparse matrix F ∈ {F z , F c , S p }. Details about the used parameters for the function agmg are discussed in Subsection 4.2. For more details about the internally used methods and the implemented syntax we refer the reader to [START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF]. Although the AGMG method can handle complex arithmetic, it needs significantly more steps to converge to the desired tolerance. Additionally, we note that agmg is a non-linear function, such that one should use a flexible iterative method, e.g., FGMRES [START_REF] Saad | A flexible inner-outer preconditioned gmres algorithm[END_REF]. However, our numerical experiments do not show any drawbacks using a standard GMRES implementation.

Chebyshev-Semi-Iteration. Although the solution of step IIIc with the symmetric positive definite mass matrix M p is relatively cheap, this can still be accelerated by using the Chebyshev-Semi-Iteration as it is described, e.g., in [START_REF] Stoll | All-at-once solution of time-dependent Stokes control[END_REF]. Numerical tests showed that one needs only 4 -6 steps to obtain a suitable result for the preconditioner, which results in a speedup that is shown in Subsection 4.2.

The next section depicts selected results to show the performance of the preconditioned iterative method.

Numerical examples

To test the efficiency of the preconditioned iterative method, the same data and configurations as in [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF] are used. After refining the initial triangulation of the reactor model in Figure 1, we end up with the variable dimensions as depicted in Table 1b. Furthermore, we define five parameter sets for different combinations of Reynolds and Schmidt numbers as shown in Table 1a. We use the MATLAB implementation of GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] to solve the saddle point systems (4) for selected ADI shifts q i that appear during the Newton-ADI iteration. Each q i is used for three ADI steps with four right hand sides every time. Thereby, the number of GMRES steps and the CPU times are measured and arithmetically averaged. The preconditioner P is evaluated as a MATLAB function handle that solves the linear system (7) using the steps [START_REF] Heinkenschloss | Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations[END_REF]. The GMRES tolerance is set to 10 -10 to ensure the same convergence of the ADI iteration that a direct solve would imply [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF]. Although a few complex ADI-shifts q i appear for each parameter set during the Newton-ADI process, the pictures only show the real parts of q i . All computations were executed in MATLAB R2012a on a a 64-bit server with 2×Intel ® Xeon ® X5650 @2.67GHz, 12 Cores (6 Cores per CPU) and 48 GB main memory available.

Influence of ADI shifts and Reynolds and Schmidt numbers

The influence of the variation of the Reynolds and Schmidt numbers as given in Table 1 is depicted in Figure 2. To obtain the best approximations for the 1a.

preconditioning steps (8a) -(8c), a direct solver is used to solve with F z , F c , and S p . It can be observed that ADI shifts -10 5 < Re (q i ) < -10 1 , between 20 -25 GMRES steps are needed. As soon as the absolute value of q i gets smaller then 10 the number of steps increases. This is a natural behavior, because the influence of the mass matrices M z and M p vanishes. Nevertheless, GMRES converges within at most 40 -80 steps for all parameter configurations. An empirical test to set: q i = -10 ∀|q i | < 10, during the Newton-ADI process showed similar ADI convergence behavior as for the original shift selection, without the drawback of higher GMRES cost for certain shifts. In summary, the derived preconditioner is suitable concerning different Reynolds and Schmidt numbers, as well as different ADI shifts.

Approximations using AMG and Chebyshev-semi-iteration

As described in Subsection 3.3, the different preconditioning steps should be solved by an easy to evaluate approximation that is accurate enough to ensure the convergence of GMRES, but avoids the use of sparse factorizations of largescale sparse matrices. We exchanged the direct solver by its approximation step by step and depict the results in Figure 3. At first, we use the MATLAB based function agmg [START_REF] Notay | An aggregation-based algebraic multigrid method[END_REF] to solve with F z and F c in (8b) and (8a) with an accuracy of 10 -10 . Depending on the used ADI shift, the function agmg needed 1 -30 steps. Thus, the times to solve the whole saddle point system with the same number of GMRES steps increased a little bit compared to the direct solver. At second, we approximately solved with S p in step IIIa using agmg as a preconditioner. This was sufficient enough to achieve the GMRES accuracy and, furthermore, decreased the time. Finally, we applied a Chebyshev-Semi-Iteration to approximately solve with M p in step IIIc. The obtained speedup finally decreased the times below the time used by the direct in each step without the loss of 4) with GMRES a representative selection of real ADI shifts from the Newton-ADI iteration for different approximations of the preconditioning steps [START_REF] Heinkenschloss | Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations[END_REF].

any accuracy in GMRES. to the above addressed problems with complex ADI shifts in agmg, we restrict our comparison in Figure 3 to a selection of real ADI shifts. The selection has been performed such that the span of all ADI shifts appearing in the entire Newton-ADI process is covered. Where those shifts clustered we only chose one representative per cluster.

At the end of this section it should be noted that the suggested preconditioned GMRES method for the considered class of saddle point problems would show its full strength in comparison to a direct solver when using finer discretizations, leading to larger dimensions, and in particular when moving to 3D problems. This will be addressed in future work.

Conclusions and outlook

We have recalled the formation of block structured saddle point systems as they arise within the Riccati-based feedback stabilization approach for coupled flow problems that avoids any explicit projection [START_REF] Bänsch | Optimal control-based feedback stabilization of multi-field flow problems[END_REF]. We were able to extend the results from [START_REF] Benner | Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow[END_REF], developed for uncoupled Stokes flow, to the coupled flow described by incompressible Navier-Stokes and a diffusion-convection equation. For that reason, the least-squares commutator approach in [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF] has been modified to approximate the shifted Navier-Stokes Schur complement. Exploiting the block structure of the arising preconditioner guarantees a fast evaluation within GMRES Ėach of the blocks can either be approximated by an AMG method or a Chebyshev-Semi-Iteration. Several numerical experiments showed that the derived preconditioning method is able to solve the arising saddle point systems efficiently independent of the different parameter settings. Only the use of complex ADI shifts during the Newton-ADI process is not yet optimally covered by this approach and will be investigated in the future.

Fig. 1 :

 1 Fig. 1: Initial triangulation of the reactor model with coordinates and boundary conditions [2]. Set Re Sc I 1 1 II 1 10 III 10 1 IV 1 100 V 10 10 (a) Different parameter settings.

10 Fig. 2 :

 102 Fig.2: Average number of GMRES steps for a representative selection of ADI shifts from the Newton-ADI iteration for the configuration sets in Table1a.

Fig. 3 :

 3 Fig.3: Average time to solve Eq. (4) with GMRES a representative selection of real ADI shifts from the Newton-ADI iteration for different approximations of the preconditioning steps[START_REF] Heinkenschloss | Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations[END_REF].

Table 1 :

 1 Test parameter settings.