Ido Cohen

Guy Gilboa

Energy Dissipating Flows for Solving Nonlinear Eigenpair Problems

This work is concerned with computing nonlinear eigenpairs, which model solitary waves and various other physical phenomena. We aim at solving nonlinear eigenvalue problems of the general form T (u) = λQ(u).

In our setting T is a variational derivative of a convex functional (such as the Laplacian operator with respect to the Dirichlet energy), Q is an arbitrary bounded nonlinear operator and λ is an unknown (real) eigenvalue. We introduce a flow that numerically generates an eigenpair solution by its steady state.

Analysis for the general case is performed, showing a monotone decrease in the convex functional throughout the flow. When T is the Laplacian operator, a complete discretized version is presented and analyzed. We implement our algorithm on Korteweg Vries (KdV) and Nonlinear Schrdinger (NLS) equations in one and two dimensions. The proposed approach is very general and can be applied to a large variety of models. Moreover, it is highly robust to noise and to perturbations in the initial conditions, compared to classical Petiashvili-based methods.

Introduction

Nonlinear elliptic equations arise in various problems in physics, e.g. for stationary solutions of equations, such as Bose-Einstein condensates (BEC), Nonlinear Schrdinger (NLS) and Korteweg Vries (KdV) [START_REF] Wang | A new approach for numerically solving nonlinear eigensolution problems[END_REF]. In this work we focus on nonlinear problems of the form,

T (u) = λQ(u), (1)
where u is a function in a Banach space U, and T and Q are (possibly) nonlinear operators. More specifically, we assume T to be a subgradient of a convex, proper, lower-semi-continuous functional J,

T (u) ∈ ∂ u J(u), (2)
where ∂ u J(u) denotes the subdifferential of J(u). On the right-hand-side of (1), Q : U → U is a bounded (possibly) nonlinear operator. We refer to functions u 1 which admit (1) as eigenfunctions, with a corresponding eigenvalue λ ∈ R. Our aim is to find a pair (u, λ) ∈ U × R that admits [START_REF] Álvarez | Petviashvili type methods for traveling wave computations: I. analysis of convergence[END_REF], referred to as an eigenpair.

Eigenpairs of nonlinear operators appear in various fields of science and engineering. Their analysis can provide deeper understanding and significant insights related to nonlinear systems. Nonlinear eigenvalue analysis is an active field of research, from both a theoretical and a computational perspective. In some nonlinear problems, such as [START_REF] Effenberger | Robust solution methods for nonlinear eigenvalue problems[END_REF] the underlying operators are linear but the dependency on λ is nonlinear. A recent review of such Nonlinear Eigenvalue Problems (NEP) appear in [START_REF] Güttel | The nonlinear eigenvalue problem[END_REF]. These studies are part of a different branch of problems (not in the scope of this paper). We examine solutions which can be formulated by Eq. [START_REF] Álvarez | Petviashvili type methods for traveling wave computations: I. analysis of convergence[END_REF], where the eigenvalue is linearly dependent. We summarize below the main related studies.

Solitary Waves as Solutions of Nonlinear Eigenvalue Problems

A pioneering work in this field, which was followed by many, is that of Petviashvili [START_REF] Petviashvili | Equation of an extraordinary soliton[END_REF]. It was aimed at finding numerical approximations of stationary solutions for the Kadomtsev-Petviashvili equation with positive dispersion (KPI equation). The method, originally, was develop to obtain stationary solutions of wave equations of the form

-M u + u p = 0, (3)
where M is a positive, self-adjoint operator and p is a constant. We note that M should be invertible, as the iterative procedure is based on its inversion. Conditions for the convergence of Petviashvili's method were established in [START_REF] Pelinovsky | Convergence of petviashvili's iteration method for numerical approximation of stationary solutions of nonlinear wave equations[END_REF]. The approach of Petviashvili was later generalized and applied to a family of nonlinear problems, such as [START_REF] Petviashvili | Solitary waves in plasmas and in the atmosphere, gordon and breach pub[END_REF][START_REF] Lakoba | A generalized petviashvili iteration method for scalar and vector hamiltonian equations with arbitrary form of nonlinearity[END_REF][START_REF] Álvarez | Petviashvili type methods for traveling wave computations: I. analysis of convergence[END_REF][START_REF] Álvarez | Numerical generation of periodic traveling wave solutions of some nonlinear dispersive wave systems[END_REF]. However, all these algorithms assume M is invertible. Moreover, it is not aimed at finding eigenpairs, but at solving a more restricted problem. When casted within the formulation of (1), the eigenvalue is set to a unit value (λ = 1). Our proposed method is based on a forward flow, and hence M is not require to be invertible. Moreover, it is aimed at finding eigenpairs, of unknown λ. The resulting eigenpair is related to an initial condition, provided by the user, which can emerge from noisy experimental data, for instance. Yang and Lakoba [START_REF] Lakoba | A generalized petviashvili iteration method for scalar and vector hamiltonian equations with arbitrary form of nonlinearity[END_REF][START_REF] Yang | Universally-convergent squared-operator iteration methods for solitary waves in general nonlinear wave equations[END_REF][START_REF] Lakoba | Conjugate gradient method for finding fundamental solitary waves[END_REF]] generalized Petviashvili's iteration method, accelerated the inverse power method and used modified conjugate gradient to find solitary waves. In [START_REF] Yang | Newton-conjugate-gradient methods for solitary wave computations[END_REF] it was suggested to combine the conjugate gradient method with accelerated inverse power method into a unified algorithm, which coincides with Petviashvili's method for small enough error. This method was shown to provide fast convergence rates. In our work we compare the numerical results to this method and to a modified version of it for adaptively computing eigenpairs . The focus of this paper is on the robustness of the methods, rather than on the convergence rate. We note that our forward flow requires considerably more iterations to converge, compared to algorithms based on inversion. However, it is much more stable and general.

Variational Formulations of Eigenpair Problems

Eigenvalue problems are often analyzed theoretically and solved numerically based on energy minimization methods. Within a variational setting, it can be shown that an eigenpair is an extremum of a generalized Rayleigh quotient [START_REF] Trefethen | Numerical linear algebra[END_REF]. This extends in a natural manner the linear case, where any eigenvector v of a Hermitian matrix A, admitting Av = λv, is an extremal point of the associated Rayleigh quotient R(v) = (v * Av)/(v * v) = λ, where v * is the conjugate transpose of v. The studies of [START_REF] Bao | Computing the ground state solution of bose-einstein condensates by a normalized gradient flow[END_REF][START_REF] García-Ripoll | Optimizing schrödinger functionals using sobolev gradients: Applications to quantum mechanics and nonlinear optics[END_REF][START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF][START_REF] Caliari | A minimisation approach for computing the ground state of gross-pitaevskii systems[END_REF] aim at finding a minimum (or a local minimum) of an energy functional associated to the eigenvalue problem. In [START_REF] Bao | Computing the ground state solution of bose-einstein condensates by a normalized gradient flow[END_REF][START_REF] García-Ripoll | Optimizing schrödinger functionals using sobolev gradients: Applications to quantum mechanics and nonlinear optics[END_REF] a constrained steepest descent is used for solving ground states of BEC. Alternative approaches, such as [START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF][START_REF] Caliari | A minimisation approach for computing the ground state of gross-pitaevskii systems[END_REF], are based on constrained energy minimization techniques with a suitable Lagrange function. In the above cases, both sides of the eigenvalue problem (1) should have an associated energy. This puts some limitations on the variety of problems that can be solved. In our work this restriction is relaxed (so only the operator T is associated with an energy term).

Eigenpairs Associated with Total Variation

The total variation (TV) functional, J T V = |∇u(x)|dx, has been thoroughly investigated in recent decades [START_REF] Burger | A guide to the tv zoo[END_REF]. Since its introduction to the image processing field for denoising and deconvolution by [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], it has been used as an edge preserving regularizer for algorithms related to stereo, optical flow, segmentation and many other computer vision tasks [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF]. Eigenpairs associated with TV were investigated by [START_REF] Meyer | Oscillating Patterns in Image Processing and Nonlinear Evolution Equations[END_REF] and [START_REF] Andreu | Minimizing total variation flow[END_REF]. It was shown that convex disk-like shapes are eigenfunctions of the nonlinear eigenvalue problem

T (u) = λu, T (u) ∈ ∂ u J T V (u),
where ∂ u J T V (u) is the subdifferential of TV. For smooth, non-vanishing gradient, we have ∂ u J T V (u) = -div(∇u/|∇u|), which is the 1-Laplacian operator. In recent years, a theory of nonlinear transforms for one-homogeneous functionals has been formulated [START_REF] Gilboa | A spectral approach to total variation[END_REF][START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF][START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF]. It is based on the analysis of nonlinear eigenvalue problems. The work of [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF] proposed a flow for finding eigenfunctions of one-homogeneous functionals, and is described in more detail below. Numerical methods for finding p-Laplacian eigenpairs were proposed in [START_REF] Yao | Numerical methods for computing nonlinear eigenpairs: Part i. iso-homogeneous cases[END_REF] and [START_REF] Hor Ák | Numerical investigation of the smallest eigenvalues of the p-laplace operator on planar domains[END_REF]. Solutions of semilinear elliptic eigen-problems were presented in [START_REF] Wang | A new approach for numerically solving nonlinear eigensolution problems[END_REF]. Ground states of generalized eigenvalue problems, which may involve also a smoothing kernel, were analyzed in [START_REF] Benning | Ground states and singular vectors of convex variational regularization methods[END_REF]. An iterative algorithm for finding nonlinear eigenpairs by an extended inverse-power method was proposed by [START_REF] Hein | An inverse power method for nonlinear eigenproblems with applications in 1-spectral clustering and sparse pca[END_REF]. All of these methods are based partially or mostly on the fact that the eigenpair satisfies an extremum of the associated (generalized) Rayleigh quotient.

Our work generalizes the flow of [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF]. It goes beyond the variational setting and allows the nonlinear operator Q in Eq. ( 1) to be very general.

Main Contributions

The main contributions of this work are as follows:

1. A new flow is introduced, which generates eigenpairs admitting [START_REF] Álvarez | Petviashvili type methods for traveling wave computations: I. analysis of convergence[END_REF]. A decrease in the convex functional J associated with the operator T is shown in the continuous setting, under mild assumptions (Sec. 2). When J is a regularizing term, such as the Dirichlet energy, we get a dissipating flow which ensures robustness against noise and stability of the algorithm.

2. The need to omit residual parts perpendicular to the operator T are explained. The concept of the Q * subspace, where all eigenfunctions must reside in, is introduced (Sec. 2.4). A complementary flow which directs the function u to this subspace is presented and combined with the main flow (Sec. 2.5, Sec. 3).

3.

A discrete explicit implementation is proposed and analyzed, when T is the negative Laplacian. We formulate a time step bound that guarantees the monotonic decrease of J throughout the flow (Section 4).

4. We employ our algorithm on KdV and NLS equation and draw a comparison between our algorithm and Newton Conjugate Gradient Methods (NCGM) [START_REF] Yang | Newton-conjugate-gradient methods for solitary wave computations[END_REF] in one and two dimensions (section 5), showing excellent stability of the flow, for various initial conditions (Sec. 5).

2 Flows for Generating Eigenpairs

Preliminary Notations and Definitions

We first introduce some standard definitions and notations used in this paper.

Let Ω be a bounded subset of R N . We use the L 2 inner product, u, v = Ω uv • dx and its associated norm u = u, u . Let J : U → R be a convex functional. We denote by ∂J(u) the subdifferential set of J at u. Let us recall the subgradient inequality: For any T ∈ ∂J(u), the following inequality holds:

J(v) ≥ J(u) + T, v -u . (4)
The set of all T for which this inequality holds is the subdifferential. A necessary and sufficient condition for u to be a global minimum of J is if zero is contained in ∂J(u). We recall the "chain rule" for the differentiation of functionals (see Lemma 3.3 of Brezis [START_REF] Brezis | Ope rateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] and recent extensions e.g. in [START_REF] Apidopoulos | On a second order differential inclusion modeling the fista algorithm[END_REF]). Let J be a convex, lower semi-continuous, proper functional, τ > 0, and u ∈ W 1,2 ((0, τ ); U). Let also T ∈ L 2 ((0, τ ); U), such that T ∈ ∂J(u(t)) a.e. in (0, τ ). Then the function

J • u : [0, τ ] → R is absolutely continuous in [0, τ ] with: d dt J(u(t)) = T, u t ∀T ∈ ∂J(u(t)) a.e. in (0, τ ). ( 5
)

Nossek & Gilboa Flow

Our work generalizes the work of Nossek & Gilboa [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF], for which J is assumed to be convex, absolutely one-homogeneous functional and Q(u) = u. We briefly describe this flow and its main properties. The flow is given by,

u t = u u - T T , T ∈ ∂ u J, u(0) = f, (6)
where f admits f = 0 and has zero mean. If J is invariant to a global constant change, i.e. J(u) = J(u + c), ∀c ∈ R, the solution u (t) of this flow has the following properties:

1. The mean value of u (t) is preserved throughout the flow:

u (t) , 1 = 0.
2. The value of J(u(t)) is decreasing with time

d dt J (u (t)) ≤ 0, a.e. in (0, ∞)
where equality is reached iff u is an eigenfunction.

3. The L 2 norm of the solution is increasing with time,

d dt u (t) 2 ≥ 0,
where equality is reached iff u is an eigenfunction.

4. A steady-state, u t = 0, holds iff u is an eigenfunction.

For further details see [START_REF] Nossek | Flows generating nonlinear eigenfunctions[END_REF]. Aujol et al. [START_REF] Aujol | Theoretical analysis of flows estimating eigenfunctions of one-homogeneous functionals for segmentation and clustering[END_REF] analyzed the flow and proposed a modification for which existence and uniqueness of the continuous flow can be shown as well as convergence of the discrete iterative algorithm. Following this, we now turn to formulate the main flow which is designed to find solutions for Eq. ( 1). We will see later that an auxiliary (complementary) flow is required for this generalization to work properly, in order to allow minimal assumptions regarding the operator Q.

The Main Flow

We introduce a flow which generates eigenpairs of the form of Eq. (1). Given some initial condition u(t = 0) = f , u(t) is evolved by the following PDE,

u t (t) = M (u(t)), (7)
where

M (u) = s Q(u) Q(u) - T (u) T (u) ,
and s = sign( Q(u), T (u) ).

Theorem 1. Assuming the flow of Eq. (7) exists and is unique and J(u) satisfies Eq. (5), the solution u(t) has the following attributes:

1. The value of J is decreasing with time,

d dt J(u) ≤ 0, a.e. in (0, ∞),
where equality is reached iff u admits Eq. (1).

2. A necessary condition for steady state u t = 0 holds iff u is an eigenfunction.

Proof. We denote

T = T (u(t)), Q = Q(u(t)).
1. From ( 5) and ( 7) we deduce,

d dt J(u) = T, u t = T, s Q Q - T T = | T, Q | Q -T ≤ 0, a.e. in (0, ∞). ( 8
)
If u is an eigenfunction in the sense of Eq. ( 1), then the eigenvalue λ can be evaluated by using [START_REF] Álvarez | Petviashvili type methods for traveling wave computations: I. analysis of convergence[END_REF] and taking the inner-product with respect to either Q or T , to have,

λ = T, Q Q 2 = T 2 T, Q = sign( T, Q ) T Q . ( 9
)
This yields | T, Q | = Q • T and therefore (d/dt)J(u) = 0. Conversely, if (d/dt)J(u) = 0 then, following (8), | T, Q | = T • Q ,
i.e., Q and T are collinear and u admits Eq. ( 1), hence an eigenfunction.

2. If u admits Eq. ( 1) then we can assign T = λQ with λ evaluated by the expression on the right of Eq. ( 9) to get u t = M (u) = 0 in [START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF]. Conversely, let u t = 0 then T and Q are collinear, hence we reach Eq. ( 1).

2.4

The Space of Solutions Q * Theorem 1 guarantees certain stability of the flow, since J(u) decreases with time. In addition, the flow stops at an eigenpair, which is our aim. We will see in the analysis below that u needs to be in a certain subspace of U, otherwise, one will attain only trivial solutions of zero eigenvalue. We begin by stating some definitions and notations. Let us define two subsets of U, T ⊥ and Q * , associated with the operators T and Q.

Definition 1. Let T ⊥ be the set perpendicular to the image of the operator T :

T ⊥ = {v ∈ U : T (u), v = 0, ∀u ∈ U}.
We note that T ⊥ is a linear subspace. This can be shown by taking two functions in that set v 1 , v 2 ∈ T ⊥ and observing that the following relation holds:

T (u), α 1 v 1 + α 2 v 2 = α 1 T (u), v 1 + α 2 T (u), v 2 = 0, ∀u ∈ U, α 1 , α 2 ∈ R.
We can now denote a basis for T ⊥ . Definition 2. Let B be an orthonormal basis of T ⊥ . Lemma 1. Let T ⊥ be defined as in Def. 1, then:

1. J(u) is invariant to addition of v ∈ T ⊥ , i.e.: J(u) = J(u + v), ∀u ∈ U, ∀v ∈ T ⊥ .
2. If the following two conditions hold: (i) The operator T is self-adjoint.

(ii) For any minimizer u of J, the subdifferential set, ∂ u J(u), has a single element. Then T ⊥ = arg min{J(u)}.

Proof.

1. As J(u) is convex then, based on the subdifferential inequality (4), we have

J(u) ≥ J(u + v) -T (u + v), v = J(u + v), ∀v ∈ T ⊥ , and
J(u + v) ≥ J(u) + T (u), v = J(u), ∀v ∈ T ⊥ .
2. If u ∈ arg min{J} and v ∈ T ⊥ then:

J(u) ≥J(v) + T (v), (u -v) =J(v) + T (v), u =J(v) + v, T (u) =J(v),
but as u minimizes J we also have

J(v) ≥ J(u), therefore J(v) = J(u) and v ∈ arg min{J}. Conversely, if u ∈ arg min{J} then T (u) = 0, since T is self-adjoint we get T (u), v = u, T (v) = 0, ∀v ∈ U.
Corollary 1. Let the conditions of Lemma 1, item 2 hold and min {J(u)} = 0, then:

1. T ⊥ is the null space of J(u).

2. J(u) is invariant to addition of elements in the null space.

We would now like to define a second significant subset, denoted by Q * .

Definition 3. Let Q * be defined by

Q * = {u ∈ U : Q(u), e i = 0 ∀e i ∈ B} ,
where B is defined in Def. 2 and e i are its orthonormal elements.

We would now like to define an energy which measures how far we are from

Q * : E(u) = 1 2 ei∈B Q(u), e i 2 . ( 10
)
Remark 1. An alternative definition of Q * is the null space of the energy functional E(u).

We observe that every v ∈ T ⊥ is perpendicular to Q(u * ) where u * is an eigenfunction with λ = 0. This is a straightforward consequence of Eq. ( 1) and Def. 1. Hence, λ Q(u * ), e i = T (u * ), e i = 0. We can summarize this by the following statement: Remark 2. All eigenfunctions with nonzero eigenvalue belong to Q * . Proposition 1. Let Q * be defined as in Def. 3 and u(t) evolves according to Eq. [START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF]. If u(t) is not an eigenfunction with λ = 0, then the following properties hold for all t ∈ (0, ∞):

1. If u(t) / ∈ Q * then (d/dt)J(u(t)) < 0.
2. If u(t) ∈ Q * then the projection of u(t) on the linear space T ⊥ is preserved throughout the flow (Eq. ( 7)). In other words:

u(t), v = const, ∀v ∈ T ⊥ .
Proof.

1. This is an immediate consequence of Theorem 1 and Remark 2. According to Remark 2 and since we assume u(t) is not an eigenfunction with zero eigenvalue, if u(t) / ∈ Q * then u(t) is necessarily not an eigenfunction. Recall Theorem 1, if u(t) is not an eigenfunction then (d/dt)J(u(t)) < 0.

2. Let us compute the derivative of u(t), v with respect to time:

d dt u(t), v = u t , v = s Q(u) Q(u) - T (u) T (u) , v = s Q(u), v Q(u) - T (u), v T (u) . If v ∈ T ⊥ then T (u), v = 0 and since u ∈ Q * then Q(u), v = 0. Consequently, d dt u(t), v = 0 ⇒ u(t), v = const.
Corollary 2. Let T be a self adjoint operator and min {J(u)} = 0. If u(t), the solution of the main flow Eq. ( 7), admits u(t) ∈ Q * , then its projection on the null space of J is preserved.

This corollary is a generalization of the first attribute of the Nossek & Gilboa flow which guarantees that the mean of u is preserved. A more significant consequence of Prop. 1 is that if u / ∈ Q * the main flow can reach a steady state only at "trivial" eigenfunctions, with λ = 0. This leads us to propose a complementary flow which retains u ∈ Q * .

Complementary flow

The above discussion emphasizes the importance of retaining u(t) to be in Q * . Following Remark 1 we define a complementary flow that aims at maintaining E(u) → 0 without changing the value of J(u). Let us first compute the variational derivative of E(u):

∂ u E = i Q(u), e i • ∂Q i (u), (11)
where

∂Q i = j (-1) j d j dx j ∂Q ∂u (j) e i ,
and u (j) denotes the jth derivative of u. We would like the complementary flow not to interfere with the main flow, and specifically not to increase J. A flow that decreases E(u) on one hand and retains J(u) unaffected on the other hand is:

u t = C(u) (12)
where

C(u) = -∂ u E + ∂ u E, T (u) T (u) 2 T (u).
Let us compute the time derivatives of J and E:

d dt J(u) = T (u), u t = T (u), C(u) = T (u), -∂ u E + ∂ u E, T (u) T (u) 2 T (u) = 0. (13)
For E we have,

d dt E(u) = ∂ u E, u t = ∂ u E, C(u) = -∂ u E 2 + ∂ u E, T (u) 2 T (u) 2 ≤ 0, ( 14
)
where the last inequality follows Cauchy-Schwarz. Thus in the flow of ( 12) E(u(t)) decreases with time whereas the value of J(u(t)) is unchanged. Note that one should verify that ∂ u E and ∂ u J are not exactly parallel with opposite directions. The relations between ∂ u E and ∂ u J are farther discussed in Sec. 5.

Combined Flows

Constrained optimization problems in variational calculus are well studied, see e.g. [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Dem'yanov | Exact penalty functions and problems of variation calculus[END_REF][START_REF] Demyanov | Exact penalty functions in isoperimetric problems[END_REF]. In light of Prop. 1 we would like our flow [START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF] to be constrained to u ∈ Q * . Otherwise, we will either not reach a steady state or reach a trivial one (an eigenfunction of zero eigenvalue), which is usually not of interest. Such constrained flows are tied, though not equivalent, to exact penalty methods, as in [START_REF] Dem'yanov | Exact penalty functions and problems of variation calculus[END_REF][START_REF] Demyanov | Exact penalty functions in isoperimetric problems[END_REF]. Inspired by these works, we couple the main flow [START_REF] Bao | Ground-state solution of bose-einstein condensate by directly minimizing the energy functional[END_REF] and the complementary flow [START_REF] Caliari | A minimisation approach for computing the ground state of gross-pitaevskii systems[END_REF] as follows,

u t = M (u) + αC(u), (15)
where α ∈ R + and M (u) and C(u) are as defined in the contexts of Eqs. ( 7) and [START_REF] Caliari | A minimisation approach for computing the ground state of gross-pitaevskii systems[END_REF], respectively. This combined flow should retain the properties of the main flow, while keeping E(u) zero (for α large enough). Let us now show that (d/dt)E(u) and (d/dt)J(u) are nonpositive:

d dt J(u(t)) = T (u), u t = T (u), M (u) ≤0 +α T (u), C(u) =0 ≤ 0. d dt E(u(t)) = ∂ u E, u t = ∂ u E, M (u) + α ∂ u E, C(u) = ∂ u E, s Q(u) Q(u) - T (u) T (u) + α -∂ u E 2 + ∂ u E, T (u) 2 T (u) 2 ≤ ∂ u E - ∂ u E, T (u) T (u) + α -∂ u E 2 + ∂ u E, T (u) 2 T (u) 2 = -∂ u E + ∂ u E, T (u) T (u) ≤0 -1 + α ∂ u E + ∂ u E, T (u) T (u) .
The inequality of the third row is based on Cauchy-Schwarz inequality,

s ∂ u E, Q(u) / Q(u) ≤ ∂ u E . If α
1 and the elements ∂ u E(u) and T (u) are not collinear (in opposite directions) we can deduce that:

d dt E (u(t)) ≤ 0.

Towards a Discrete Setting

We propose in Algorithm 1 the general discrete iterations that approximate the combined flow of [START_REF] Demyanov | Exact penalty functions in isoperimetric problems[END_REF]. It is based on an explicit scheme for the time discretization. The time step parameters dt C and dt M are not resolved here. A full analysis and bounds on these time steps is given next, for the common case of T (u) = -∆u.

Here, instead of choosing a specific α, the complementary flow is evolved (for each time step of the main flow) until a certain accuracy is reached, with respect to being close enough to Q * . This is analogous to having a time dependent α(t) in [START_REF] Demyanov | Exact penalty functions in isoperimetric problems[END_REF]. An illustration of the combined flow is shown in Fig. 1a where M (u) and C(u) are invoked alternatingly. As discussed above, in the continuous setting C(u) does not change the value of the functional J(u) and is along its level lines. However, M (u) might cause u to be out of Q * . The complementary flow, C(u), "returns" u to Q * leaving J(u) unchanged. Note that in the exact penalty algorithm [START_REF] Dem'yanov | Exact penalty functions and problems of variation calculus[END_REF] the flow is run in several iterations with increasing α. The combined flow is repeated until the solution is close enough to an eigenpair in terms of the respective angle (See Fig. 1b and Algo. 1).

Algorithm 1 Eigenpair Generating Flow

1: Inputs: u 0 2: Initialize: δ, θ th stop conditions Compute θ 0 according Eq. ( 16)

3: while θ k > θ th do 4: while E (u) > δ do 5: u k,l = u k,l-1 + u k,l-1 t • dt C (complementary flow) 6:
end while 7:

u k+1 = u k + u k t • dt M (main flow) 8:
Compute θ k+1 according Eq. ( 16) 9: end while Our stopping condition is based on the absolute angle between T (u) and Q(u),

|cos(θ)| = | T (u), Q(u) | T (u) • Q(u) . ( 16
)
When the angle θ is zero, T and Q are collinear, and an exact eigenpair is attained. For numerical purposes, we seek a good approximation and stop when θ is below a certain threshold. Further discussions related to the stopping condition are in A.

Realization for T (u) = -∆u

As many phenomena are based on eigenpairs, in the sense of Eq. ( 1), with T (u) the negative Laplacian operator, we devote the rest of this paper to that particular case. In this case J is the Dirichlet energy,

J(u) = 1 2 ∇u 2 .
It is clear that J admits the conditions of Eq. (1) and Eq. ( 5). The null space of the Dirichlet energy is all functions in U of constant value. According to Cor. 2 it coincides with T ⊥ . Therefore,

T ⊥ = {v ∈ U : v = const}.
The orthonormal basis of

T ⊥ is B = 1 |Ω| .
The energy functional E(u) is defined by,

E(u) = 1 2|Ω| Q(u), 1 2 , ∂E = 1 |Ω| Q(u), 1 ∂Q,
and ∂Q is the variational derivative of Q(u), 1 . Naturally, Q * is all u ∈ U such that Ω Q(u) dx = 0. The main and complementary operators, respectively, are,

M (u) = s Q(u) Q(u) + ∆u ∆u , (17)
where s = sign( Q(u), -∆u ), and

C(u) = -∂ u E + ∂ u E, ∆u ∆u 2 ∆u. ( 18
)
In the following section we resolve the time step sizes for the discrete implementation when T is the negative Laplacian operator.

Implementation for the Laplacian operator

For the discrete implementation of the flow we follow the general algorithm presented in the previous section. However, instead of performing a full iterative process of the complementary flow for each main time step, we perform a single step. It is shown that for T = -∆, with an appropriate time step of the complementary flow, the energy E vanishes, up to a first order Taylor approximation.

We suggest the following iterations to approximate the combined flow: Initialize with u 0 , iterate until convergence,

u k+ 1 2 =u k + M (u k ) • dt M u k+1 =u k+ 1 2 + C(u k+ 1 2 ) • dt C . (19)
In the rest of this section we set bounds for dt M and dt C , to ensure stable flows. These bounds guarantee that in every cycle (from k to k + 1) we obtain J(u k ) -J(u k+1 ) ≤ 0 and E(u k ) ∼ = 0, ∀k.

Characterization of J(u) and E(u) within the Discrete Setting

We now proceed with analyzing the behavior of J(u k ) and E(u k ) using the iterations of ( 19), where M (u k ) and C(u k ) are defined by ( 17) and [START_REF] Ekeland | Convex analysis and variational problems[END_REF]. For simplicity, we use the same notations of gradient (∇) and Laplacian (∆) for the discrete operators.

Main Flow Iteration

We first analyze the influence of an iteration of the main flow (first row of Eq. ( 19)) on the functional J. We use the standard first order backward difference time-step approximation of explicit schemes. We examine the change in J at iteration k, denoted by D k J , following a time step of the main flow:

D k J :=J(u k+ 1 2 ) -J(u k ) =J u k + dt M • M (u k ) -J(u k ).
Using the identity J(u) = 1 2 -∆u, u yields:

D k J = 1 2 -∆ u k + dt M • M (u k ) , u k + dt M • M (u k ) --∆u k , u k =dt M -∆u k , M (u k ) + 1 2 dt 2 M ∇M (u k ) 2 . ( 20
)
According to Theorem 1 the expression ∆u k , M (u k ) is non-negative, thus, if

0 ≤ dt M ≤ 2 ∆u k , M (u k ) ∇M (u k ) 2 (21)
then D k J ≤ 0. Let us now examine the influence of the main flow on the functional E. We assume E(u k ) is approximately zero (up to a first order Taylor approximation), to be justified in the next paragraph. The Taylor approximation of E(u k+ 1 2 )is given by:

E(u k+ 1 2 ) = 1 2|Ω| Q(u k+ 1 2 ), 1 2 = 1 2|Ω| Q(u k + dt M M (u k )), 1 2 = 1 2|Ω| Q(u k ), 1 + dt M ∂Q(u k ), M (u k ) + O(dt 2 M ) 2 ∼ = 1 2|Ω| Q(u k ), 1 + dt M ∂Q(u k ), M (u k ) 2
Under the assumption E(u k ) is approximately zero (i.e. Q(u k ), 1 ∼ = 0) we have,

D k E :=E(u k+ 1 2 ) -E(u k ) ∼ =E(u k+ 1 2 ) ∼ =dt 2 M • ∂Q(u k ), M (u k ) 2 2|Ω| . (22)

Complementary flow iteration

We now analyze the influence of the complementary flow (second row of Eq. ( 19)) on the energies J and E. Here we begin with the analysis of E(u).

D

k+ 1 2 E :=E(u k+1 ) -E(u k+ 1 2 ) =E(u k+ 1 2 ) + dt C • ∂E(u k+ 1 2 ), C(u k+ 1 2 ) + O(dt 2 C ) -E(u k+ 1 2 ) ∼ =dtC • ∂E(u k+ 1 2 ), C(u k+ 1 2 ) .
The step size for which E(u k+1 ) vanishes, as a first order approximation, is therefore

dt C = - E(u k+ 1 2 ) ∂E(u k+ 1 2 ), C(u k+ 1 2 ) . (23)
The step size dt C is positive since the denominator is negative (according to Eq. ( 14)). Finally, we examine the influence of the complementary flow on the functional J. The change of J along one time step of the complementary flow is

D k+ 1 2 J :=J(u k+1 ) -J(u k+ 1 2 ) = 1 2 -∆u k+1 , u k+1 --∆u k+ 1 2 , u k+ 1 2 = 1 2 -∆(u k+ 1 2 + dt C • C(u k+ 1 2 )), u k+ 1 2 + dt C • C(u k+ 1 2 ) --∆u k+ 1 2 , u k+ 1 2 .

14

Note that according to Eq. ( 13) ∆u k+ 1 2 , C(u k+ 1 2 ) = 0 and therefore

D k+ 1 2 J = dt 2 C ∇C(u k+ 1 2 ) 2 .
Using Eq. ( 23),

D k+ 1 2 J = E 2 (u k+ 1 2 ) ∇C(u k+ 1 2 ) 2 ∂E(u k+ 1 2 ), C(u k+ 1 2 ) 2 . ( 24
)
Algorithm 2 Eigenpair Generation for the Laplacian Operator 1: Inputs: u 0 2: Initialize:

C(u 0 ), Eq. ( 18) dt C , Eq. ( 23)

u 1 = u 0 + C(u 0 ) • dt C Set threshold θ th θ 1 ,
Eq. ( 16) 3: while θ th < θ k do 4:

u k+ 1 2 = u k + dt M • M (u k ), Eq. (30) 5:
C(u k+ 1 2 ), Eq. ( 18)

6:
dt C , Eq. ( 23)

7:

u k+1 = u k+ 1 2 + C(u k+ 1 2 ) • dt C 8:
θ k+1 , Eq. ( 16) 9: end while 4.1.3 Conditions for the monotonic decrease of J Following the above constraints on the time step sizes of both flows we can now analyze the change in J for the entire flow, Eq. ( 19),

J(u k+1 ) -J(u k ) =D k+ 1 2 J + D k J . ( 25
)
We would like to find a bound on dt M which ensures J(u k+1 ) -J(u k ) ≤ 0.

Substituting for D

k+ 1 2 J
and D k J by Eqs. ( 24) and ( 20), respectively, we have,

J(u k+1 ) -J(u k ) =E 2 (u k+ 1 2 ) ∇C(u k+ 1 2 ) 2 ∂E(u k+ 1 2 ), C(u k+ 1 2 ) 2 + dt M -∆u k , M (u k ) + 1 2 dt 2 M ∇M (u k ) 2 .
The term E(u k+ 1 2 ) is evaluated according to the first order approximation of Eq. ( 22) to obtain,

J(u k+1 ) -J(u k ) ∼ =dt 4 M • ∂Q(u k ), M (u k ) 4 4|Ω| 2 • ∇C(u k+ 1 2 ) 2 ∂E(u k+ 1 2 ), C(u k+ 1 2 ) 2 + dt M -∆u k , M (u k ) + 1 2 dt 2 M ∇M (u k ) 2 . ( 26
)
A deeper analysis shows that when the directions of ∂E(u) and ∆u are far enough from being collinear the fourth order term is small and thus can be omitted. After neglecting the fourth order term on the right-hand-side of Eq. ( 26) the inequality coincides with Eq. ( 21). We now use the following inequality based on the operator norm of the discrete gradient:

∇u 2 ≤ 4d h 2 u 2 , ( 27
)
where d is the dimension and h is the spatial grid size. Based on the expression on the right-hand-side of Eq. ( 21), and the definition of M (u) for the Laplacian case, ( 17), we have

2 ∆u k , M (u k ) ∇M (u k ) 2 ≥ h 2 2d ∆u k , M (u k ) M (u k ) 2 = h 2 2d | ∆u k ,Q(u k ) | Q(u k ) + ∆u k 2 + 2 | ∆u k ,Q(u k ) | ∆u k • Q(u k ) = h 2 4d ∆u k .
(28) Therefore, to ensure the expression in ( 26) is negative, we set the following bound on the time step size,

dt M ≤ h 2 4d ∆u k . ( 29
)
According to Eq. ( 20), D k J depends on dt M quadratically, therefore, the optimal dt M is half of the upper bound. Consequently, we choose

dt M = h 2 8d ∆u k . ( 30
)
Typical values of dt M in our experiments were in the range [0.002, 0.6] for one dimension and [0.05, 0.13] for two dimensions. These values highly depend on the noise level. A detailed description of the algorithm, in the case of the Laplacian operator, is presented in Algorithm 2.

Interpretations

Let us revisit Eq. ( 1) and the main flow. Given u k we can evaluate the eigenvalue according to Eq. ( 9) and define the error of the solution u k as:

err(u k ) = -∆u k -λ(u k ) • Q(u k ) ( 31
)
which is zero if u k is an eigenfunction. In addition, by using the expression of dt M in Eq. ( 30), we can calculate the change of u at the kth step as,

u k+ 1 2 -u k =M (u k ) • dt M = h 2 8d ∆u k + sign( Q(u k ), -∆u k ) ∆u k Q(u k ) Q(u k ) = h 2 8d ∆u k + λ(u k ) • Q(u k ) = - h 2 8d • err(u k ). (32)
We note that [START_REF] Lakoba | A generalized petviashvili iteration method for scalar and vector hamiltonian equations with arbitrary form of nonlinearity[END_REF] suggests a general eigenvalue approximation by: λ(u) = -∆u, f (u) / Q(u), f (u) with an arbitrary operator f (u) as long as this ratio exists. In particular we can have f (u) = u. Then we can replace Line 4 in Algo. 2 with Eq. ( 32) using other eigenvalue evaluations. We can conclude that any flow in the sense of Eq. ( 32) with an eigenvalue approximation for which Theorem 1 holds, the later properties of Sec. 2 (Lemma 1 and Prop. 1) also hold. This conclusion is valid not only for T (u) = -∆u but for any operator T (u) admitting Eq. ( 1). However, different step size constraints should be found.

Applications

In this section we present several numerical experiments using the proposed method. We compare our results with NCGM introduced in [START_REF] Yang | Newton-conjugate-gradient methods for solitary wave computations[END_REF]. NCGM assumes the eigenvalue is given, therefore we supply the correct eigenvalue to the algorithm. In addition, we compare our method to an adaptive version of NCGM when the eigenvalue is not given (thus the task is to find a full eigenpair, as our flow does). The eigenvalue is approximated at each iteration based on [START_REF] Lakoba | A generalized petviashvili iteration method for scalar and vector hamiltonian equations with arbitrary form of nonlinearity[END_REF]. We denote this method by Adaptive NCGM (ANCGM).

1 Dimension

Soliton Equations

In 1895 Korteweg Vries (KdV) formulated a mathematical model of waves on shallow water surfaces which were previously described by Russell. Later studies have shown that the KdV equation is the continuum limit of a discrete nonlinear mass-spring model [START_REF] Pava | Nonlinear dispersive equations: existence and stability of solitary and periodic travelling wave solutions[END_REF]. The formulation by Kruskal and Zabusky [START_REF] Zabusky | Interaction of" solitons" in a collisionless plasma and the recurrence of initial states[END_REF] for such phenomena was:

u t + uu x + δ 2 u xxx = 0,
with δ a small real scalar. Reformulating this expression for a stationary wave yields the following equation:

-u XX = λ -cu + u 2 2 , ( 33
)
where c is the wave velocity, X = x -ct and λ = δ -2 . Naturally, λ can be understood as an eigenvalue. In this example we have Q(u) = -cu + u 2 /2. The solution to this equation models well a family of solitary waves referred to as solitons, (named by Kruskal and Zabusky). In this case an analytic solution can be formulated, which makes it feasible to compare the experimental results to the ground truth. The solution for (33) is:

u(X) = 3c • sech 2 √ c • λX 2 .
The implementation of our flow is based on Algo. 2, where for solving Eq. ( 33) we have:

Q(u) = -cu + u 2 2 , ∂Q(u) = -c + u, ∂ u E = Q(u), 1 (-c + u) .
The operators in NCGM [START_REF] Yang | Newton-conjugate-gradient methods for solitary wave computations[END_REF] are:

L = ∆ + λ (c -u) , M = ∆ + λ • c,
where λ is given. The notations L and M are taken from [START_REF] Yang | Newton-conjugate-gradient methods for solitary wave computations[END_REF] (not to confuse with M (u) of our algorithm). In the adaptive version, ANCGM, λ is approximated by λ(u) = -∆u, u / Q (u) , u .

Experimental setup

In our experiments we initialize with the a noisy version of the solution, by adding white Gaussian noise N of a certain variance to the analytic solution for some given λ. Note that we also allow some bias in the simulated measurements so the noise mean in not zero. The noise parameters are denoted by a normal distribution N ∼ N (µ, σ 2 ), with mean µ and variance σ 2 . We employ the aforementioned algorithms with a stopping condition θ th ≤ 2 • (where θ is the angle between T and Q, as defined in ( 16)). In Fig. 2 the results of the generated eigenfunctions are shown. The noise parameters are µ = 0.1 with variance values of σ 2 = 0.2 and σ 2 = 0.36, for the first and second experiments, respectively. For this relatively simple case we compare our results to ANCGM, where both algorithms try to find an eigenpair (λ is not known). We show that our algorithm is quite stable, for different levels of noise. ANCGM on the other hand, collapses to the trivial solution (-∆u = λQ(u) = 0) for high noise levels (and is less accurate also for moderate noise). We should emphasize that our algorithm runs considerably more iterations than ANCGM. However, our iterations are very fast, as they are based on explicit forward-flow computations, compared to computationally-intensive iterations of 2a we summarize the results of ANCGM and our flow when the initial function u 0 , red solid line, is an eigenfunction with additive white Gaussian noise (of mean of 0.1 and variance of 0.2). In Fig. 2b we repeat the experiment with higher level of noise (variance is 0.36).

ANCGM which involve inversion. Overall, the computations are comparable. The running time (using Matlab code) in the experiment with variance 0.2 are t = 0.336[sec] for ANCGM and t = 0.1[sec] for our algorithm. For variance 0.36, it was t = 1.7[sec] for ANCGM and t = 1.19[sec] for ours.

According to Theorem 1, the main flow reaches a steady state when an eigenpair is attained. Thus, it is expected that if u t = 0 then sign( Q(u), -∆u )Q(u)/ Q(u) = -∆u/ ∆u . In Fig. 3 we plotted sign( Q(u), -∆u )Q(u)/ Q(u) and -∆u/ ∆u . Moreover, on the top left side of these plots we present the approximated eigenvalue λ = -∆u, u / Q (u) , u and the computed pointwise ratio λ(x) = -∆u(x)/Q(u(x)) (u(x) is the result of our algorithm at convergence). One can observe these graphs coincide very well. In Fig. 4 the values of J(u(t)) and the mean of u(t) as a function of t are shown. Our theoretical analysis predicts J(u(t)) should be decreasing with time (Theorem 1 in the time continuous case and the bound on dt M of Eq. ( 29) in the discrete setting). Moreover, since u(t) is very close to Q * , the mean of u(t) should be almost constant in time (Prop. 1, item 2). It is shown that the numerical implementation approximates well the theoretical analysis.

Low SNR

An essential assumption of algorithms based on linearization, such as NCGM, is that the initial guess u 0 is close enough to the true solution. In this example we try to demonstrate what happens when this assumption is not valid. In this experiment we initialize with the eigenfunction corresponding to λ = 2 and c = 0.25 and add noise with normal distribution N ∼ N (0.25, 0.25). We compared our algorithm and NCGM, giving both of them this initial function (and additionally providing NCGM the correct λ). Fig. 5 shows the results of this experiment. In Fig. 5a the eigenfunction u * , the initial function u 0 , our summarizes the result of NCGM and our flow when the initial function is far from the eigenpair (high level of noise). Whereas NCGM collapses to the trivial solution our flow reaches a meaningful solutions. In Fig. 5b it is demonstrated that a good eigenpair approximation is attained. In Fig. 5c the evolution of J and the mean of u as a function of number of iterations are presented.

algorithm and NCGM results are plotted. Our algorithm converges to another eigenfunction whereas NCGM collapses to a trivial solution. Fig. 5b illustrates the similarity of the two opposite normalized counterparts of the flow u t , i.e. sQ(u)/ Q(u) and -∆u/ ∆u showing we indeed approximate an eigenfunction well. In Fig. 5c the value of J and the mean of u are plotted as a function of iteration.

One Dimensional Nonlinear Schrdinger (NLS) Equation

In this part we would like to emphasize the simplicity of our approach. Our algorithm is easy to adjust to other applications and the adjustment amounts to redefining Q(u) and its variational derivative. We demonstrate this by the NLS equation as discussed in [START_REF] Shabat | Exact theory of two-dimensional self-focusing and one-dimensional self-modulation of waves in nonlinear media[END_REF][START_REF] Zakharov | On the complete integrability of a nonlinear schrödinger equation[END_REF]:

i ∂u ∂t + ∂ 2 u ∂x 2 -k u 2 u = 0.
To find a stationary solution we assign u(x, t) = u(x)e iµ(u)t . This yields -u xx = -µ(u)u -ku 3 , which can be rewritten as -u xx = λ(-μ(u)u -ku 3 ) with μ(u) = 1 -3u 2 and k = 2. Finally we reach the problem formulation of (1),

-u xx = λ(u 3 -u).
Therefore, we have

Q(u) = u 3 -u, (
Q(u) = -V 0 sin 2 x + sin 2 y u + σu 3 + µu, ∂Q(u) = -V 0 sin 2 x + sin 2 y + 3σu 2 + µ.
The results of our algorithm and NCGM with u 0 are summarized in Fig. 7. In Fig. 7a the initial guess u 0 is shown. The results of NCGM and our flow are presented in Figs. 7b and7e, respectively, and the angle evolutions along iterations are presented in Figs. 7c and7f. In Fig. 7d we show the behavior of the functional J and the mean of u as a function of time.

We performed a similar experiment of the problem with noise level of N (0, 0.0125). The result of the experiment are similarly organized in Fig. 8. It is easy to see that NCGM converges approximately to linear case with very small amplitude signal.

Decomposition into Eigenfunctions

A signal may be modeled in some cases by a linear combination of nonlinear eigenfunctions. In this case we would like to decompose it into its basic elements. We suggest below a simple iterative algorithm to perform that (Algo. 3). In Fig. 9 we show an example of such a decomposition for the case of KdV. We compute the eigenfunctions φ 1 , φ 2 , φ 3 with respective eigenvalues λ 1 , λ 2 , λ 3 . We used as initial condition a combination of two eigenfunctions of KdV (with eigenvalues λ = 0.25, λ = 4). One of them is translated in space (therefore their correlation is low). In the decomposition results, the first two eigenfunctions are very similar to the original ones. The third one is of small amplitude. In that case the linear part becomes dominant and the problem reduces to finding eigenfunctions of the ordinary heat equation. Thus, we get a harmonic function.

Conclusions

In this paper we presented a new nonlinear flow which can generate a very broad family of nonlinear eigenfunctions. The flow reaches a steady state only for eigenpairs. Moreover, the process is a forward flow, which is based on the Algorithm 3 Nonlinear Eigenfunction Decomposition 1: Inputs: u 0 2: while u k-1 > th do 3:

Find an eigenpair (φ k , λ k ) using Algo. original operators and does not need their inversion. We thus obtain a very stable process and our assumptions regarding the operators are less restrictive. When T is the negative Laplacian operator, we prove that the flow decreases the Dirichlet energy. This yields a highly stable smoothing process, which is robust to noise.

We introduced the space Q * , where all nontrivial eigenfunctions reside in, and explained the requirement of the flow to be constrained to this space. A complementary flow was proposed in order to achieve this requirement. We have tested our algorithm on several nonlinear eigenvalue problems in one and two dimensions, and have shown the stability of the algorithm and its robustness to various initial conditions. Finally, it was demonstrated how potentially our process can be used iteratively, in order to decompose a signal into a linear combination of eigenpairs.

A Considerations for the stopping condition

Implementations of numerical solutions are concerned with issues related to the appropriate stopping condition. As T (u) and Q(u) are fully correlated when eigensolution of Eq. ( 1) is attained, the most natural criteria is the angle between them (see Eq. ( 16)). Other criteria worth mentioning are (d/dt)J(u) Eq. ( 8) and M (u) which are shown to be zero when an eigenpair is reached. The norm of the main flow can be calculated as:

M (u) = √ 2 1 - | Q(u), T (u) | Q(u) • T (u) . ( 34
)
If u / ∈ Q * then ∃e ∈ B s.t. Q(u), e = 0. Let us denote the projection of Q(u) on T ⊥ as:

Q ⊥ (u) = ei∈B e i • Q(u), e i = 0, and the "residue" of Q(u) as:

Q (u) = Q(u) -Q ⊥ (u).
As E(u) is positive we get the following inequalities:

|cos(θ)| = | T (u), Q(u) | T (u) • Q(u) = T (u), Q (u) T (u) • Q(u) ≤ Q (u) Q(u) = Q (u) Q (u) 2 + 2E(u)
and similarly

M (u) = √ 2 1 - | Q(u), T (u) | Q(u) • T (u) ≥ √ 2 1 - Q (u) Q (u) 2 + 2E(u) , d dt J(u) ≤   Q (u) Q (u) 2 + 2E(u) -1   T (u) .
As long as E(u) > 0 our analysis has shown that a precise eigenfunction cannot be obtained. However, one can set a small threshold on E(u), for example, if E(u) ≤ δ we have:

M (u) th = √ 2 1 - Q (u) Q (u) 2 + 4δ d dt J(u) th =   Q (u) Q (u) 2 + 4δ -1   T (u) |cos(θ)| th = Q (u) Q (u) 2 + 4δ . ( 35
)
We note that the above discussion concerns the true evolution in the time continuous setting. However its implications can naturally affect the numerical algorithms.

  (a) J functional contour (b) The combined flow as finite state machine

Figure 1 :

 1 Figure 1: An illustration of the infinitesimally approximation of the combined flow.

Figure 2 :

 2 Figure2: A comparison of eigenpair generating algorithms for KdV equation with different noise levels. In Fig.2awe summarize the results of ANCGM and our flow when the initial function u 0 , red solid line, is an eigenfunction with additive white Gaussian noise (of mean of 0.1 and variance of 0.2). In Fig.2bwe repeat the experiment with higher level of noise (variance is 0.36).

  (a) N ∼ N (0.1, 0.2) (b) N ∼ N (0.1, 0.36)

Figure 3 :

 3 Figure 3: We show here that our algorithm indeed reaches an eigenpair. The results of T / T and sQ/ Q are shown when the algorithm converges for the KdV equation. One can observe they coincide very well. The black solid line stands for sQ(u)/ Q(u) and dashed yellow line is -∆u/ ∆u . In left top side the approximated λ * = -∆u * , u * / Q (u * ) , u * is plotted with dashed red line and the solid blue line is the ratio -∆u/Q(u).

Figure 4 :

 4 Figure 4: The value of J(u k ) and the mean of u k (right top axis) as a function of iteration number k when applying our flow on KdV for the experiments depicted in Fig. 2.

Figure 5 :

 5 Figure5: Experiment when the premise u * -u 0 << 1 does not hold. Fig.5asummarizes the result of NCGM and our flow when the initial function is far from the eigenpair (high level of noise). Whereas NCGM collapses to the trivial solution our flow reaches a meaningful solutions. In Fig.5bit is demonstrated that a good eigenpair approximation is attained. In Fig.5cthe evolution of J and the mean of u as a function of number of iterations are presented.

  a) The initial function u 0 (b) Result of NCGM (c) Angle Vs. Iterations, NCGM (d) J(u) and mean of u (right top axis) (e) Our flow's result (f) Angle Vs. Iterations, ours

Figure 7 :

 7 Figure 7: Two dimension NLS.

  (a) Initial function u 0 (b) Result of NCGM (c) Angle Vs. Iterations, NCGM (d) J(u) and mean of u (right top axis) (e) Our flow's result (f) Angle Vs. Iterations, ours

Figure 8 :

 8 Figure 8: Two dimension NLS with noisy initial function.

Figure 9 :

 9 Figure 9: Function decomposition. On the left column, the dashed red line is the initial signal u 0 , which is a combination of two eigenfunctions (λ = 0.25, λ = 4). The blue solid lines are the reconstructions of the signal up to the kth eigenfunction, i.e. fk = k i=1 α k • φ k as described in Algo. 3. On the right column, we show the eigenfunctions found by Algo. 3, in the respective order. The first two functions are very similar to the original ones. The third one is a residual function of low amplitude (dominated by the linear part).

2

 2

4 : 2 5 :

 425 α k = u k-1 ,φ k φ k u k = u k-1 -α k • φ k 6: end while 7: Compute the approximated function fk = k i=1 α k • φ k

We initiate the algorithm with u 0 = u * + N where u * is the eigenfunction u(x) = √ 2 sech ( √ λx) and N ∼ N (0, 0.45), see Fig. 6. We similarly initiate the algorithms NCGM and ANCGM. The first diverges and the second one converges to the trivial solution, λ = 0. Both of these results are not in the same scale as the initial signal and therefore are not plotted.

Two Dimensions NLS Equation

The following experiment is based on Example 3.4 from [START_REF] Yang | Newton-conjugate-gradient methods for solitary wave computations[END_REF] and demonstrates the performance of our algorithm in 2D. We apply our algorithm on NLS equation with periodic potentials which models nonlinear light propagation as well as Bose-Einstein condensate's dynamics in optical lattices. It is formulated as:

This equation is discussed in details in [START_REF] Dalfovo | Theory of bose-einstein condensation in trapped gases[END_REF][START_REF] Morsch | Dynamics of bose-einstein condensates in optical lattices[END_REF][START_REF] Skorobogatiy | Fundamentals of photonic crystal guiding[END_REF] and admits a rich variety of solitary waves of the form of U (x, y, t) = u(x, y)e -iµt , where u(x, y) satisfies the equation:

In our example µ = 4.11, V 0 = 6 , σ = 1 and Dirichlet boundary conditions are used. There is no known analytic solution for these equations. The initial


Ph.D Denis Gilbert
email: denis.gilbert@irstea.fr

Ph.D Edo Abraham
email: e.abraham@tudelft.nl

Idel Montalvo
email: imontalvo@ingeniousware.net

Ph.D Olivier Piller
email: olivier.piller@irstea.fr

Iterative Multi-Stage Method for a Large Water Network Sectorization into DMAs under Multiple Design Objectives

Keywords: Water Distribution Network, Graph Partitioning, DMA, Sectorization, Rehabilitation Design, Agent Swarm Optimization, Pressure Management, Water Quality, Engineering Judgment

This paper considers the sectorization of a large water distribution network into DMAs and simultaneously optimizes rehabilitation of the network with new pipes, control valves, and storage tanks. Since the available water resources are much smaller in the dry season, both the design and operational settings are optimized to satisfy water demand, water quality and pressure constraints, and efficiency indices under stringent conditions. Because of the heterogeneity of the multiple decision variables and the complicated way they interact through the multiple objectives (some complimentary and some conflicting), it is not possible to fully automate the simultaneous sectorization, rehabilitation and operational optimization. Therefore, we employ a multi-stage approach where engineering judgment and network graph simplification and visualization tools are employed to find a good feasible solution that is used as a first guess for further optimization of sectors and operational settings, to achieve feasible solutions with better cost of implementation, demand similarity among DMAs and better pressure uniformity in operations. A multi-objective Agent Swarm Optimization (ASO) framework is used to iteratively change the sectors at the boundaries. For the final configuration, sequential linear programming is used to find optimal valve and pump settings.

INTRODUCTION

This manuscript presents our work on the optimal sectorization of the E-town network, which was part of the Battle of Water Networks DMAs (BWNDMA 2016), a special contest with 14 participants that was presented at WDSA 2016 in Cartagena, Colombia, the sixth one in a series of contests since 1985 [START_REF] Giustolisi | Operational and tactical management of water and energy resources in pressurized systems: Competition at WDSA 2014[END_REF]. Firstly, we describe the network and the problem objectives and constraints.

Then, the methodology employed is explained, which consists of some engineering judgment and network analysis tools to simplify the problem, and the use of mathematical optimization and agent swarm optimization for design and operational decisions. We subsequently present our results and discuss lessons of general validity that have been learned from this exercise. Finally, we conclude with a summary.

The schematic shown in Figure 1 represents a large water network of the E-Town city in Colombia. As a result of a swelling tourist industry and overall economic growth, the city is no longer able to meet its increasing water demand. Having modelled the network and made an inventory of the forecasted demands, demand patterns, existing pump and tank characteristics, and the actual controls of valves, the municipality has determined that the current DMA configuration and operational and tactical management of the network do not allow it to meet demands in dry and rainy seasons (BWNDMA 2016). For example, as depicted in Figure 2, some of the tanks are empty and out of use because there are considerable differences in the pressure conditions of the city, and because the water use is not efficient. Therefore, the objective of this work was to sectorize the network into a new DMA configuration that allows the water utility to satisfy its customer needs while keeping a minimal number of DMAs (each with a similar number of users or demand), guarantying pressure uniformity across the municipality, meeting regulatory water quality specifications, and ensuring an efficient system operation during both the dry and wet seasons of the year. (𝑉 $%,$ -𝑉 :;<,$ -𝑉 => )

#2345 $@1

, where 𝑉 $%,$ , 𝑉 :;<,$ are the volume of water entering and leaving the ith DMA ands 𝑉 => volume consumed per DMA, averaged over all DMAs..

iii)

To minimize solution implementation cost 𝐶𝐶 %'< . This consists of the cost of pipe replacements or installations (as a function of length and discrete set of available diameters), the cost of control valve installations (as a function of diameters), and the cost of new tanks installed (as a function of volume).

iv)

To minimize the pressure Uniformity index (𝑃𝑈 %'< ) during the rainy and dry seasons as two separate objectives. These objectives specify that all demand nodes in the network should have similar pressure and as close as possible to the minimum set pressure constraint. Please see This design problem is a very difficult since there are multiple decision variables that are not easy to incorporate in one optimization problem and because the multiple objectives interact in complicated ways. For example, one of the design stage decision variables is the isolation valves to be closed (i.e.

sectorization is done by closing some of the pipes on the boundary of DMAs). Since each DMA is specified to have only a maximum of two pressure controlled inlets, this will implicitly determine the number of PRVs and so the objective 𝐶𝐶 %'< , in addition to 𝐷𝑀𝐴 $%&'( and 𝐷𝑆. The new tanks to be installed and their volume is not easy to incorporate in this same sectorization problem as it is not clear a priori, where additional storage is needed or what the volumes should be. Similarly, we have the discrete decisions of introducing additional new pipes. These design decisions also affect the optimal level of water 𝑃𝑈 %'< , 𝑊𝐴 %'< , 𝑂𝑝𝐶𝐻 %'< and to a lesser extent 𝑂𝑝𝐶𝐻 %'< , that can be achieved through operational optimization of PRVs, FCVs and pump settings for the two seasons. Therefore, because of the magnitude of the search space and the large number of decision variables, the multiple interactions between the objective functions (some complimentary and some conflicting) and heterogeneity of the decision variables, it is not possible to derive a design process that is fully automatic.

In the literature, several methodologies have been proposed to tackle the problem of water network sectorization, often with respect to one decision variable and fewer objectives. For example, motivated to maximize security against contamination events, (Di Nardo et al. 2013) employ graph theory principles and a heuristic optimization to form isolated district meter areas, each of which is supplied by its own source (or sources). Although they use graph theory tools for sectorization, in [START_REF] Scarpa | Elementary DMA Design of Looped Water Distribution Networks with Multiple Sources[END_REF] the two objectives are enhancing the quality of delivered water and reducing the risk of contaminant spread. Most other work also uses graph theory and energy criteria (Di Nardo & Di Natale 2011; [START_REF] Giustolisi | Operational and tactical management of water and energy resources in pressurized systems: Competition at WDSA 2014[END_REF] and some have applied multi-agent systems for finding the best combination of sectors (Herrera et al. 20112, Montalvo et al. 2014) for other multi-objective problems.

There are also approaches that try to tackle the problem directly with heuristic optimization algorithms (Di Nardo et al. 2014) and recently with concepts derived from social network theory [START_REF] Diao | Automated creation of district metered area boundaries in water distribution systems[END_REF].

A reference to a plethora of other related work on sectorization approaches and specific objectives for sectorization can be found in these references and literature cited therein.

The approach followed in our work is not exclusively based on any one method presented in the literature but a combination of several approaches together with engineering judgement to solve the overall problem in multiple stages. As previously demonstrated in [START_REF] Khedr | Comparing optimization techniques with an engineering judgment approach to WDN design[END_REF] for a water distribution system rehabilitation problem, the use of engineering judgement to simplify a rather complex design and operational optimization problem can be very effective. We use graph theory tools to better understand the network topology, visualize the elevation map and identify redundant and isolated network elements, which were corrected before automating some processes. We then optimize iteratively the DMAs using agent swarm optimization. The operational settings of valves and pumps are determined using convex optimization tools. In the following section, we describe the methodologies in detail.

METHODOLOGY Engineering principles and Porteau

The Porteau software solution (2017) is a hydraulic toolkit for Water Distribution Analysis that is developed by IRSTEA research institute in France. It was designed with the help of a Computer-Aided Software Engineering tool and object oriented programing. This software is an alternative to the well-known EPANET solution [START_REF] Rossman | EPANET 2 Users manual[END_REF] and it provides different additional tools for analyzing complex WDNs such as simplification of the network, detection of isolated parts and duplicated pipes, and it owns a stochastic module for assessing the hydraulic state at peak period (Piller and Bremond, 2002). The technologic choice and the software possibilities are described in [START_REF] Piller | Porteau: An Object-Oriented Programming Hydraulic Toolkit for Water Distribution System Analysis[END_REF].

Porteau was first used in the model validation step to identify parts that are isolated from water resources and to check errors in the data file itself. In the original network file provided for the competition, a set of pipes was discovered to be disconnected from the rest of the network. This may cause problems for a hydraulic solver with the presence of customer nodes with no reference of head.

One pipe, originally closed (i.e.: Pipe '8872' in the INP file provided), was opened to avoid problems in hydraulic analysis and force the hydraulic solver to compute a correct rather than arbitrary piezometric head. Another issue we identified was regarding the two tanks 'Tank_3' and 'Tank_4' connected by a short pipe. This caused mass oscillation between the two tanks at each iteration and extended period simulation (EPS) hydraulic equations like in Porteau and EPANET were not able to find a correct physical solution since they ignore pipe inertia; such problems would not be an issue if a rigid water-column model was used [START_REF] Piller | Slow transient pressure driven modeling in water distribution networks[END_REF][START_REF] Nault | Improved Rigid Water Column Formulation for Simulating Slow Transients and Controlled Operations[END_REF]. We also found that this modeling flaw increased considerably the executable time of the EPS. Since the supply flow direction is from 'Tank_3' to 'Tank_4', both in rainy and dry seasons, we placed check valves on the outgoing pipes from both tanks. This prevention of two-way flow in the numerical model prevented the mass oscillation, both decreasing the computational time and giving correct and physically meaningful results. In valve and pump setting optimization problems (see last subsection), the same set of issues were resolved by putting a unidirectional constraint on the flow of the links originating from 'Tank_3' and 'Tank_4'.

Several parallel pipes (with the same positions) were found in the original network data. An attentive check was made to ensure it was not an error coming from the GIS. In Porteau, a connectivity analysis highlights parallel pipes and visualizes them. Then they are manually checked to see if they are different pipes. Since we found no duplicated pipes, the parallel pipes were kept in the model. Five out of service tanks were transformed into junction nodes. Practically, this allows water to cross the valve chamber and bypass the tank. Another abandoned tank ('Tank_7') was isolated because its elevation was too high to be supplied with sufficient pressure coming from another tank or through nodes around it.

The second step was to analyze the supply routes (Figure 2) and elevation map (Figure 3) to understand the supply network (see Figure 1). In the rainy season, the water has three origins. The source in the south, Cuza WTP, constitutes more than half of the supply and provides water to a large part of the city. The source at the central east of the network, Bachue WTP, is at a low elevation and can only supply the low elevation parts close to it, including Tank_14. Bochica WTP, although it has half the capacity of the south supply, it has the highest elevation and so supplies areas with a large variation in elevation across the north of the network. for isolation from the rest of the network, and one main entry point was chosen for each DMA and tagged "meters/compteurs". A graph theoretic tool was used to check that each DMA was connected only through the tagged inlets: using the network incidence matrix of the connected components and the PRIM algorithm (a greedy depth-first-search algorithm) the inlets to each DMA were confirmed [START_REF] Bartnik | Graphes, algorithmes, logiciels[END_REF]. For double-checking the results, a modified version of the Lee algorithm [START_REF] Lee | An Algorithm for Path Connections and Its Applications[END_REF] was used for identifying existing sectors in the network model. The algorithm uses graph decomposition and the information of the network elements for traversing in depth search until the border of a potential sector is reached.

The third step is iterative and includes the use of two different steady state hydraulic models to check the pressure at demand nodes during the peak period. The conventional deterministic models predict too large flow rates in main pipes, which are improbable to occur simultaneously. On the contrary, flow rates in branched pipes supplying few consumers are underestimated by the latter. More representative values can be derived by considering stochastic demands and upper limits of different confidence intervals; we use the stochastic model Opointe in Porteau that exploits peak demand diversity/simultaneity curves [START_REF] Piller | A Stochastic Model for Peak Period Analysis of Pipe Networks[END_REF] to more accurately approximate the flow rate at pipes as a function (𝑎𝑛 + 𝑏 𝑛 ) of the number of domestic consumers, n, served. This way a less pessimistic pressure is calculated in the core network that corresponds to a satisfaction risk.

Additionally, EPANET was used for simulation from within a c-code with conventional analysis of low and peak periods. If the lower constraint limit of 15m for pressure was not satisfied it was necessary to adjust the DMA boundary and specify a new DMA with a more homogenous pressure. Following this, another steady state for low demand consumption was run to check for the pressure upper limit of 60m.

Among the main results, it was necessary to reinforce the network with two new pipes of length approximately 800m and 1.2km (connecting node '5675' to nodes '4466' and '5374', respectively) in order to limit the head loss in supplying the DMA below 'Tank_15'. This high elevation area (e.g., node

'4252' has an elevation of 47m) couldn't be supplied with a pressure above 15m without violating the 60m pressure constraints upstream of it -feasible solutions are obtained only with this reinforcement with new parallel pipes.

The fourth step was to tune the setting of the flow control valves (FCVs) and of the pressure regulating devices (PRVs) that are added at entry points of DMAs. The main role of the FCV is to control the distribution of water and it was possible to set each FCV to the supply of one or several DMAs. This has simplified the problem of finding feasible initial settings for the FCVs. An automatic procedure that requires solving a convex optimization problem was also devised for that task and linked to the pressure uniformity criterion (see the next subsection). This has led to a feasible solution for the rainy season.

The average demand of each DMA and the mass balance of each tank were used to determine the initial setting of the FCV.

The solution of the rainy season with pumps working and reduced flows from the water treatment plants was adapted for the dry season. It was necessary to make two kinds of change: closing and opening valves to introduce new water routes (shortcuts) and to supply some tanks with insufficient water inflow with alternative sources. We added a pressure-sustaining valve to raise the pressure in order to limit the water coming from the south source 'Cuza WTP' and force it to go to 'Tank_16' and decreased the flow going to 'Tank_1'. We also reinforced with set pipes of diameter 762 mm and total length 721 m, starting from 'Fagua_Pump_Station' to transport the supply to 'Tank_12'.

Sequential Convex Programming for Operational Optimization

In the problem description and rules for the competition (BWNDMA, 2016), the objectives 𝐷𝑀𝐴 $%&'( , 𝐷𝑆, 𝐶𝐶 %'< , 𝑂𝑝𝐶𝐻 %'< , (and implicitly 𝑃𝑈 %'< ) were taken into account only while choosing the DMA configurations in the design iterations as they are not a function of operational setting changes.

The operational decisions affect the two objective functions relating to pressure uniformity (𝑃𝑈 %'< ) and water age (𝑊𝐴 %'< ). By posing the water quality objective function 𝑊𝐴 %'< as a constraint, i.e. all designs would need to satisfy the stipulated water age limit, this objective function could be removed from the multi-objective design problem. For the dry season, in the operational optimization, we consider the simultaneous tuning of the settings for the FCVs, the PRVs and the flows from the fixed-speed pumps with the two objective functions 𝑃𝑈 %'< and 𝑊𝐴 %'< such that the following (linear constraints) are satisfied:

• pressure constraints of 15m-60m for the specified nodes at all time steps,

• flow constraints on all water sources between 0 and the maximum specified, and additional constraint for the two Mohan sources to sum to 206 L/s,

• tank levels between 10-90 % full, at all time steps, and

• mass balance at all nodes.

In the rainy season, all the above except for the flow from the pumps are considered.

This optimization problem is a difficult nonlinear programming problem (NLP), with non-convex constraints. An optimization method based on the strictly feasible sequential convex programming (SCP) described in [START_REF] Wright | Control of water distribution networks with dynamic DMA topology using strictly feasible sequential convex programming[END_REF] was used. This SCP method solves the non-convex NLP problem by sequentially solving convex approximations (linearized sub-problems). The convexity of the approximations means each sub-problem can be solved accurately and efficiently. In the optimization problem considered here, we start from feasible PRV and FCV and pump flow settings that were determined in the DMA design stage. Let 𝑥 ∶= 𝑞 O ℎ O 𝜂 O O be the vector of the flows in all links (𝑞), the heads at all unknown head nodes (ℎ) and the headloss across the PRVs (𝜂), respectively.

The NLP for calculating η and flows through FCVs at each time instant is as follows:

(1)

where the vector equation 𝑔 𝑥 contains the head loss across each link and mass balance equations at all nodes, respectively. Where the link indexed in 𝑔 𝑥 is a pump, the head gain across the pump is represented by the pump curve equations [START_REF] Rossman | EPANET 2 Users manual[END_REF]. The upper and lower bounds on the flows 𝑞 $ represent constraints on flows through FCV from different sources, and the minimum and maximum pressure levels at demand and non-demand nodes and tank level constraints are imposed using the constraints on the heads ℎ S . The respective 𝜂 $ values are constrained as non-negative for PRVs and positive for pressure sustaining valves (PSVs), where the direction of flow is already determined. For entrance to DMAs with only flow meters, the flow direction constraints can be omitted. At each iteration of the SCP, the nonlinear objective function 𝑓(⋅) and the nonlinear constraints 𝑔(𝑥) are approximated by their linearization and the following problem is solved:

(2)

where 𝑥 V is the variable at the k-th iteration of the SCP, 𝛻𝑓 𝑥 V and 𝛻𝑔 𝑥 V are the gradients (i.e. partial derivatives) of the objective function and constraints, respectively, computed at current iterate 𝑥 V and 𝑥 is the solution we seek, at which we shall compute the linearized gradients at the next iterate 𝑘 + 1. In addition to these, we impose flow direction constraints for the supply pipes (i.e., the FCVs from the sources), and inlet PRVs by setting the corresponding lower bound 𝑞 Y to zero. min ⌘,q,h f (h; ⌘, q)

subject to: g(x) = 0

q i < q i  qi , 8i 2 N F CV h j < h j  hj , 8j 2 N N ⌘ i 0, 8i 2 N P RV min x rf (x k ) T x subject to: rg(x k ) T (x x k ) + g(x k ) = 0, q i < q i  qi , 8i 2 N F CV h j < h j  hj , 8j 2 N N ⌘ i 0, 8i 2 N P RV
In the optimization problem (1), the flows from the pumps at Fagua and Mohan water treatment plants and the head difference across the pumps are decision variables, constrained by the corresponding pump curves. In the SCP subproblems, the heads are approximated using the linearization of the pump curve at a given flow. In Figure 4, we show different approaches that can be used to approximate the pump curves. A simple approximation of the pump-curve (left side of Figure 4a and Figure 4b) can be derived by drawing a line between the head and flow axis intercepts; this underestimates the head of the pump for a given flow and sometimes can result in violation of the maximum head constraints. Alternatively, (right side of Figure 4a and Figure 4b), a least-squares fit of the pump curve can overestimate or underestimate the head depending on the flow rate. A first-order Taylor linearization of the pump-curve around flow 𝑞 V (at the k-th iteration of the SCP) gives an overestimator for the head that is tight near 𝑞 V . From simulations, we found that the last linearization had the least hydraulic infeasibility, a faster convergence for the SCP and therefore it was adopted; the details of this pump approximation approach and how the tank levels are discretized using a first-order finite-difference model is described in [START_REF] Menke | Demonstrating demand response from water distribution system through pump scheduling[END_REF]. Three more aspects of the SCP application in our implementation that are important to mention are:

• The objective function 𝑊𝐴 %'< was ignored for the SCP because it was mostly zero when perturbing the feasible initial solution from the design subproblem. Moreover, ignoring this objective function means that the water quality simulations are not needed at each iteration of the SCP but only an extended period hydraulic simulation to check feasibility in taking a step in decision variables. This results in faster CPU time of the SCP. 𝑊𝐴 %'< was checked once the SCP had converged.

• The PRV technology used by the utility are specified to have fixed-head settings that change only once from dry to rainy season and vice versa. Therefore, in the operational optimization we enforce a constraint on the pressure at the downstream nodes of the regulating PRVs to be fixed in time. Most times, this caused infeasibility in the linearized subproblems of the SCP, even when starting from a feasible point. The linearized feasibility spaces at different times (and therefore at different states for flows and pressures) seem to no longer have an intersection with a fixed head constraint for the PRVs. This was relaxed by putting a penalty term for temporal changes in the downstream heads of PRVs instead. However, it is not clear how to find appropriate penalty weights to guarantee that we don't introduce too much suboptimality in the objective 𝑃𝑈 %'< . As a result, the SCP's results do not always give significant improvement. It is also not clear how optimal our initial feasible guess was since we are using a local optimization solver. Therefore, this is left as an interesting study for future work.

• The nonlinear objective 𝑃𝑈 %'< was studied by decomposing it into the two parts:

(3)

Where P min was fixed to be 15 mH2O, P i,j is the pressure at node i and time j; and P Avj is the average pressure at time j for all the nodes.

Using the semidefiniteness of the numerator in the second term of equation ( 3) and the positivity of its denominator, perhaps it is possible to show convexity of 𝑃𝑈 %'< , which simulations seem to imply. Nonetheless, we use a linearization of 𝑃𝑈 %'< in the SCP.

Network design optimization by an ASO method

Creating sectors in a water network is a challenging problem from the engineering point of view. Even more so when it is approached using optimization methods trying to find optimal DMA configurations:

• Classic optimization based on derivatives would be very limited for solving this problem as they need rigorous formulation and the decision variable should be continuous.

• Checking all alternatives by enumerating solutions is infeasible in terms of calculation time as the problem scales combinatorically (Pecci et al., 2016) • The use of "pure" evolutionary algorithms would have difficulties to find good solutions (in a reasonable amount of time/resources) for a "large size" problem like this one, because for example GA is sensitive to the initial population used and Genetic algorithms do not scale well with complexity.

• Hybrid evolutionary approaches that incorporate engineering experience and other more deterministic methods for solving sub-problems are considered as an option.

Agent Swarm Optimization (ASO) (Montalvo et al., 2010) introduced a new agent that reproduces the behavior of the Louvain method for community detection (Blondel et al., 2008). The Louvain method maximizes a modularity index resulting in the best possible grouping of nodes:

Q = 1 [\ 𝐴 $S - V ] V [\ 𝛿 𝑐 $ , 𝑐 S
, where (4)

• 𝐴 $S represents the edge weight between nodes i and j,

• 𝑘 $ and 𝑘 S are the sum of the weights of the edges attached to nodes i and j, respectively,

• m is half the sum of all edge weights in the graph,

• 𝑐 $ and 𝑐 S are the communities of the nodes i and j respectively,

• 𝛿 is a Dirac delta function that gives 1 if ci == cj, else zero.

Classic applications related to social network can use, for example, the number of messages between community members for weighting edges between nodes. In the case of water network sectors, this research has been using the water demand at nodes as a first weighting factor. For joining smaller communities when the modularity index was not able to be increased anymore the Z coordinate (or elevation) was used as weighting factor. These decisions had the intention to find sectors with similar demand and with nodes located at similar altitude. The action of the agent reproducing the Louvain method was combined with other type of agent reproducing the behavior of a modified Lee algorithm [START_REF] Lee | An Algorithm for Path Connections and Its Applications[END_REF] for finding the limits of potential sectors and deciding where to put PRV valves. This behavior was also used in this research for checking solutions created just based on engineering experience. ASO calculation and network visualization were done using Water-Ing, software from Ingeniousware GmbH. The software is open to extensions and available for free in a community edition (https://ingeniousware.com/software/about?name=watering).

The settings of the PRVs placed at DMA boundaries for separating sectors are not included directly in the process of creating sectors. Nevertheless, the results of the objective function after the PRV settings are defined is sent back to the ASO as reference in order to improve the sectors created.

During the time of this research, ASO was partially used to improve the sectorization solution as because of limited time, it was not possible to automate the whole process of creating sectors taking into account all the details of the objective functions. Nevertheless, the agents included in ASO helped on the evaluation/verification of the ideas emerging from engineering experience that were used to solve the problem.

SUMMARY and DISCUSSION OF RESULTS

In Figure 5, we show DMA configurations of the main North part of the network where the colors were chosen from a library by the Porteau software to distinguish each DMA. It can be observed that there is a large variation of elevation within small distances in the North. As a result, it was necessary to sectorize the areas into smaller DMAs so that the pressure can be controlled within the given limit of (15m-60m pressure) using PRVs, each DMA acting like an "irrigation terrace" used in agriculture as we go uphill. Although one of the objectives by the battle organizers was to have the smallest number of DMAs possible (but above 15 in number), this also explains the relatively large number of DMAs in this area, which has resulted in some DMAs having relatively smaller average volume demand compared to bigger DMAs with higher volume of consumption. Therefore, the demand similarity objective DS was increased to maximum to make sure the conflicting pressure similarity objective is reduced, and pressure constraints were satisfied.

For the south part (see Figure 6) and in our sectorization solution, the water is coming from Tank 15 (depicted with the dark-blue square in the center bottom). It is important to highlight the necessary reinforcement/reorganization is the most important intervention we have done in term of implementation costs. It consists of pipe reinforcement with roughly 2km of pipes (top thick green line) to limit the head loss at the lower elevation and satisfy pressure constraints at the high altitude in the ends of DMA 4. Because of the large variation in elevation in in this area, we created three DMAs to sectorize this area into three areas of more homogeneous elevation. from Bachue (the very right of Figure 7). In the dry season, this supply is halved in capacity and so was complemented by Mohan pumping station. In this competition, the pumping stations are specified to contain only fixed-speed pumps. To increase redundancy and flexibility in supply, we recommend pumps (to have higher head specifications) and introduce variable speed pumps, which can be operated more efficiently with resultant savings of energy [START_REF] Wu | Incorporation of variable-speed pumping in multiobjective genetic algorithm optimization of the design of water transmission systems[END_REF]). This can be exploited even further if the water resources (aquifers) at Mohan have capacity to supply more water. In summary, the following were also achieved in the final configuration and optimized operational settings:

• In Table 1, we show the decomposition of 𝑃𝑈 %'< into the first and second elements of the sum in Equation ( 3), labelled PUnet_Pmin PUnet_Pav, respectively. We note that some DMAs have high deviation from the minimum pressure because they are either in a sector of low elevation near the sea ('Secteur_005', 'Secteur_010' in NW) or possess high variation in elevation and are close to a source supplying other DMAs ('Secteur_012' and 'Secteur_020' at the entrance of the Northern part, see also Figure 8). Note also that 'Secteur_012' is close to 'Tank_16' with nodal pressures close to the 60m upper limit at times). • The water age objective (𝑊𝐴 %'< ) is close to zero (0.005 hours) in both as there are only 20 violations of the 60-hour maximum limit.

• By making use of a graph simplification tool and engineering judgement, it was possible to identify the main flows and bottlenecks. This has allowed changing the network with minimal pipe interventions (11 in total but with only 5 with significant lengths). Similarly, the operational changes from one season to another have been kept

to only 14 (with five main pipes status changes, eight FCV setting changes and three PRV setting changes). The details of all these are summarised in the Excel files submitted as supplementary material together with the optimized network models in EPANET INP format.

• Similarly, we have tried to minimize the capital cost of intervention with valves (see Table 2).

SUMMARY AND PERSPECTIVES

The objective of this paper is to share our experience with the civil engineering and academic community in facing the challenging Battle of Water Networks District Meter Areas problem.

We have adopted three main strategies:

The first one consists of using our experience and skills in network modeling as much as possible. For that, it was very important to understand how the supply of tanks and demand nodes may work in rainy and dry seasons. The steps taken included:

• model validation,

• supply route and altimetry analysis,

• use of a deterministic and stochastic hydraulic model to test the pressure constraints for candidate DMA configuration,

• tuning the setting of the FCVs and of the PRVs.

This strategy has allowed the derivation of a feasible solution with regards to the pressure and tank level-range constraints. It was the most time consuming and has required inventiveness. Our experience was also used to simplify the problem in terms of the number of variables and objectives and to help the ASO with appropriate rules.

The second strategy was to use a metaheuristic agent-swarm method to explore other solutions that are not intuitive and so were not exploited in engineering judgment. The method moves some elementary pieces of network and analysis the different combinations. A single hydraulic steady state was run to check the pressure constraint feasibility and water age constraints. No substantial improvement was achieved due to the lack of automation of the overall sectorization process with ASO.

The third and last strategy was to solve the operational optimization problem for the setting of PRVs and FCVs, and for the pump working point by sequential convex programming method. Further studies are needed to guarantee significant decreases in the objective or use relaxation approaches for global optimality studies [START_REF] Pecci | Outer approximation methods for the solution of co-design optimisation problems in water distribution networks[END_REF].

The final solution consists of a reinforcement with 3km of pipes at a cost of 141 k$. In sectorizing the network, we placed three kinds of control and metering devices: PRVs that are regulating the DMAs There is no tank cost since we do not recommend new tanks or increasing the capacity of existing ones; the total cost is only approximately 349 k$. The water age at every node is lower than 60 hours so that our solution is optimal for this objective, with only insignificant violations of 0.05 hours aggregated over all nodes over seven Satisfying the objectives to the average pressure and increasing the pressure uniformity in each and across the DMAs has led to the creation of additional smaller sized DMAs compared to the solution derived at the first stage. Finally, we have 31 DMAs in total.

Future research will study further the use of ASO and deterministic optimization tools, as well as their efficient coupling to further explore the space of possibilities. We have also provided as supplementary material to this paper, the INP files of our results and a summary of the results in Excel files containing the scores of the objective functions set, the definition of the sectors, nodal pressures, the tank level at each time-step, and all these for the two seasons. Visual Basic and C codes for running EPANET and generating the analysis on the networks are available on request.

It should be noted that, although the algorithms and software tools used have were vital in solving a large-scale problem that is too complex to tackle manually, they were complimentary to sound engineering judgement rather than completely substituting it. An essential core of the solutions was also the result of applying experience and engineering knowledge.

FIGURE CAPTIONS

TABLES

Figure 1 :

 1 Figure 1: E-Town network Map with resources and used tanks.

  vi)Minimize the total number of operational changes (𝑂𝑝𝐶𝐻 %'< ) in the valve settings (for pressure reducing valves (PRVs) and flow control valves (FCVs)) and opening and closing of isolation valves for boundary pipes from dry season to rainy season.

Figure 2 :

 2 Figure 2: The original supply routes for the E-Town network (BWNDMA 2016)

Figure 3 :

 3 Figure 3:. Altimetry for Downtown E-Town network (network pipes in white; sectors delimited by

Figure 4 :

 4 Figure 4: (a) Different linear approximations for the Pump Curve of Fagua compared to the given

Figure 5 :

 5 Figure 5: Final DMA splitting result for the main North part of E-Town network

Figure 7 :

 7 Figure 7: Final DMA splitting for the central part of E-Town network

Figure 8 :

 8 Figure 8: Mean pressure deviation from PUnet_Pmin for North part of E-Town network

  Although we use 52 valves in the network, 13 of the PRVs located at entry point of the DMAs are completely open and not doing any regulation. They only represent water meters at their respective DMA entry point, where we are required to measure the volume of water going through. EPANET does not offer the possibility to include water meters and that's why they are represented in our solution as totally open PRV without any regulation.

  are placed at the sector entrances, metering devices that we represent by open PRVs, and open PRVs at the outlets for metering what is leaving the DMAs. The valves cost 208 k$ in total.

Figure 1 :

 1 Figure 1: E-Town network Map with resources and used tanks.

Figure 2 :

 2 Figure 2: The original supply routes for the E-Town network (BWNDMA 2016)

Figure 4 :

 4 Figure 4: (a) Different linear approximations for the Pump Curve of Fagua compared to the given

Figure 5 :Figure 8 :Figure 9 :

 589 Figure 5: Final DMA splitting result for the main North part of E-Town network

Table 1 .

 1 Comparison of 𝑃𝑈𝑛𝑒𝑡 objective function elements for all the 31 DMAs.

		demand
	DMA	m3/7 days	nbNode	nbDemandNode	PUnet_Pmin	PUnet_Pav
	Secteur_028	1549.31	25	25	0.6987	2.2601
	Secteur_032	1759.63	29	29	0.0691	5.4905
	Secteur_021	2006.85	31	26	0.8269	2.5155
	Secteur_029	2301.72	26	26	0.0235	5.4707
	Secteur_011	4908.42	44	37	0.7508	2.4374
	Secteur_022	5191.81	48	36	1.3208	2.5766
	Secteur_004	5383.49	118	62	1.4741	4.0094
	Secteur_027	7897.8	88	74	2.3604	4.0048
	Secteur_007	9387.84	193	162	6.1162	5.6021
	Secteur_017	9578.36	100	94	3.2133	3.5324
	Secteur_019	10279.01	50	39	0.7096	3.1027
	Secteur_034	10739.73	14	12	0.1432	2.3412
	Secteur_030	11089.48	246	224	8.4334	7.414
	Secteur_018	12335.58	145	115	3.7259	5.3124
	Secteur_035	14782.33	34	30	0.4497	3.1806
	Secteur_025	15004.06	103	97	0.4734	9.1526
	Secteur_033	15567.98	38	34	1.2717	2.1086
	Secteur_026	23114.2	448	377	10.6525	4.4648
	Secteur_023	54159.39	514	440	10.3694	5.6524
	Secteur_031	54601.91	498	410	12.8299	8.7903
	Secteur_002	57795.21	557	409	12.5499	8.492
	Secteur_005	71648.25	572	507	9.5529	10.0936
	Secteur_012	73703.9	934	773	37.0166	19.2645
	Secteur_010	77886.19	852	713	25.5174	7.4012
	Secteur_001	81519.5	249	186	3.0262	7.093
	Secteur_020	83853.16	1055	898	20.8845	9.9377
	Secteur_009	84905.37	606	501	16.4216	2.6549
	Secteur_008	85158.52	250	198	4.678	3.6356
	Secteur_006	87365.95	338	217	6.5152	2.3366
	Secteur_003	88830.7	708	567	17.5792	3.0187
	Secteur_013	90513.91	1045	891	27.9382	4.5598
	noDma	9112.37	1164	41	3.0443	18.2869

Table 2 .

 2 Intervention types and associated cost.

	Intervention type	Implementation cost in $
	Pipes	141188
	Tanks	0
	Valves	208142
	Total	349330