Shilpa Khatri

Arnold D Kim

Ricardo Cortez

Camille Carvalho

Close evaluation of layer potentials in three dimensions

Keywords: Nearly singular integrals, close evaluation problem, potential theory, boundary integral equations, numerical quadrature

come

Introduction

The close evaluation problem arises in boundary integral methods for numerically solving boundary value problems for linear, elliptic partial differential equations. In boundary integral equation methods, the solution of the boundary value problem is given in terms of double-and single-layer potentials, which are integrals of a kernel multiplied by a density over the boundary of the domain. The kernel for the single-layer potential is the fundamental solution of the elliptic, partial differential equation and the kernel for the double-layer potential is the normal derivative of that fundamental solution. As a result, each of these kernels has an isolated singularity at a known point on the boundary. When evaluating double-and single-layer potentials at points close to the boundary, the associated kernel is regular, but is sharply peaked because the evaluation point is close to its singular point on the boundary. For this reason, we say that double-and single-layer potentials evaluated close to the boundary are nearly singular integrals.

Nearly singular integrals are more challenging to compute numerically than weakly singular ones. When computing weakly singular integrals, one must explicitly address the singularity in the kernel. For example, weakly singular integrals can be computed accurately using high-order product Nyström methods [START_REF] Atkinson | The Numerical Solution of Integral Equations of the Second Kind[END_REF][START_REF] Delves | Computational Methods for Integral Equations[END_REF]. These product Nyström methods are often used to solve the boundary integral equations for the density. However, for a nearly singular integral, there is no singularity to address. Nonetheless, without effectively addressing the peak in a nearly singular integral, one must increase n(y  ) B averaging operation the order of the quadrature rule to obtain accuracy commensurate with evaluation points that are far from the boundary.

The close evaluation problem has been studied extensively. Schwab and Wendland [START_REF] Schwab | On the extraction technique in boundary integral equations[END_REF] have developed a boundary extraction method based on a Taylor series expansion of the layer potentials. Beale and Lai [START_REF] Beale | A method for computing nearly singular integrals[END_REF] have developed a method that first regularizes the nearly singular kernel and then adds corrections for both the discretization and the regularization. The result of this approach is a uniform error in space. Helsing and Ojala [START_REF] Helsing | On the evaluation of layer potentials close to their sources[END_REF] developed a method that combines a globally compensated quadrature rule and interpolation to achieve very accurate results over all regions of the domain. Barnett [START_REF] Barnett | Evaluation of layer potentials close to the boundary for Laplace and Helmholtz problems on analytic planar domains[END_REF] has used surrogate local expansions with centers placed near, but not on, the boundary. This work led to the work by Klöckner et al. [START_REF] Klöckner | Quadrature by expansion: A new method for the evaluation of layer potentials[END_REF] that introduces Quadrature By Expansion (QBX). QBX uses expansions about accurate evaluation points far away from the boundary to compute accurate evaluations close to it. The convergence of QBX has been studied in [START_REF] Epstein | On the convergence of local expansions of layer potentials[END_REF]. Moreover, fast implementations of QBX have since been developed [START_REF] Klinteberg | A fast integral equation method for solid particles in viscous flow using quadrature by expansion[END_REF][START_REF] Rachh | Fast algorithms for quadrature by expansion i: Globally valid expansions[END_REF][START_REF] Wala | A fast algorithm for Quadrature by Expansion in three dimensions[END_REF], and rigorous error estimates have been derived for the method [START_REF] Klinteberg | Error estimation for quadrature by expansion in layer potential evaluation[END_REF]. Recently, the authors developed a method that involved matched asymptotic expansions for the kernel [START_REF] Carvalho | Asymptotic analysis for close evaluation of layer potentials[END_REF]. In that method, the asymptotic expansion that captures the peaked behavior of the kernel can be integrated exactly and the relatively smooth remainder is integrated numerically, resulting in a highly accurate method.

Most of these aforementioned studies are for two-dimensional problems. There are fewer results for three-dimensional problems. Beale et al. [START_REF] Beale | A simple method for computing singular or nearly singular integrals on closed surfaces[END_REF] have extended the regularization method to threedimensional problems. Additionally, QBX has been used for three-dimensional problems [START_REF] Klinteberg | A fast integral equation method for solid particles in viscous flow using quadrature by expansion[END_REF][START_REF] Wala | A fast algorithm for Quadrature by Expansion in three dimensions[END_REF]. In principle, the matched asymptotic expansion method by the authors can be extended to three dimensional problems. However, we do not pursue that approach because the method we present here is more direct and simpler to implement.

We study the close evaluation of double-and single-layer potentials in the context of the interior Dirichlet and exterior Neumann problems for Laplace's equation in three dimensions, respectively. Compared to the close evaluation of layer potentials in two dimensions, we find that layer potentials in three dimensions are actually better behaved because the double integral over the boundary naturally includes an averaging operation about y , the point on the boundary at which the kernels have a singularity (see Fig. 1). This averaging operation is a crucial aspect of the layer potential's representation of a harmonic function that smoothly approaches its boundary data. To find this inherent averaging operation, we perform a local analysis of the contribution of the layer potentials in the vicinity of y . This local analysis provides valuable insight into the close evaluation problem and suggests a natural method for numerically computing nearly singular layer potentials in three dimensions. We also find that this method can be used to compute weakly singular integrals, in particular for numerically solving the boundary integral equation for densities.

The remainder of this paper is as follows. We precisely define the close evaluation problem for the double-and single-layer potentials in Section 2. In Section 3, we compute the asymptotic behavior of the contributions to these layer potentials made by a small region containing the point where the kernels are peaked. By doing so, we find a natural numerical method to evaluate these nearly singular integrals. This numerical method is detailed in Section 4. We discuss in Section 5 how this method can also be used to compute matrix entries in the Galerkin method to solve the boundary integral equations for the densities. We give several examples demonstrating the accuracy of this numerical method in Section 6. Section 7 gives our conclusions. A gives details of how we rotate spherical integrals in the results presented in this paper, and B gives a useful derivation of the spherical Laplacian.

Close evaluation of double-and single-layer potentials

as the double-layer potential,

u(x) = 1 4π B n(y) • (x -y) |x -y| 3 µ(y)dσ y , x ∈ D, (2.2)
with n(y) denoting the unit, outward normal at y ∈ B, and dσ y denoting the surface element. Here, µ is the density for the double-layer potential satisfying

-∂ n u = g on B, (2.4)
as the single-layer potential,

v(x) = 1 4π B ρ(y) |x -y| dσ y , x ∈ E. (2.5)
Here, ρ(y) is the density for the single-layer potential satisfying

1 2 ρ(y ) - 1 4π B n(y ) • (y -y) |y -y| 3 ρ(y)dσ y = g(y ), y ∈ B, (2.6)
Boundary integral equation methods involve solving (2.3) and (2.6) for µ and ρ, respectively. These results are then used to compute u using (2.2) for any desired point in D and v using (2.5) for any desired point in E. Boundary integral equations (2.3) and (2.6) can be solved to very high accuracy using Nyström methods [START_REF] Atkinson | The Numerical Solution of Integral Equations of the Second Kind[END_REF][START_REF] Delves | Computational Methods for Integral Equations[END_REF]. Thus, the focus turns to accurately computing the double-and singlelayer potentials, (2.2) and (2.5), respectively. The high accuracy of Nyström methods transfers over to the evaluation of these layer potentials except for close evaluation points.

To define a close evaluation point precisely, let 0 < ε 1 denote a small, dimensionless parameter, and consider

x = y ∓ ε n , (2.7)
for x ∈ D and x ∈ E, respectively, with y ∈ B denoting a point on the boundary, n = n(y ) denoting the unit, outward normal from D at y , and denoting a characteristic length of B, e.g. its radius of curvature at y . Before substituting (2.7) into (2.2), we first make use of Gauss' law [START_REF] Guenther | Partial Differential Equations of Mathematical Physics and Integral Equations[END_REF],

1 4π ∂D n(y) • (x -y) |x -y| 3 dσ y =      -1 x ∈ D, -1 2 x ∈ B, 0 x ∈ E, (2.8)
to write the double-layer potential evaluated at the close evaluation point given in (2.7) as

u(y -ε n ) = -µ(y ) + 1 4π B n(y) • (y -y -ε n ) |y -y -ε n | 3 [µ(y) -µ(y )] dσ y .
(2.9)

We call (2.9) the modified double-layer potential resulting from applying a subtraction method [START_REF] Delves | Computational Methods for Integral Equations[END_REF].

For the close evaluation of the single-layer potential, we substitute (2.7) into (2.5) and obtain

v(y + ε n ) = 1 4π B ρ(y) |y -y + ε n | dσ y . (2.10)
Our objective here is the accurate and efficient numerical evaluation of (2.9) and (2.10). The challenge in this problem is that in the asymptotic limit as ε → 0, those integrals are nearly singular. Setting ε = 0 yields weakly singular integrals which can be numerically evaluated accurately. In fact, this is what is done in the Nyström method to solve boundary integral equations (2.3) and (2.6). However, when ε is small, but finite, the functions to be integrated are not singular, but they are sharply peaked about y = y which makes computing (2.9) and (2.10) challenging.

3 Motivation: local analysis of the double-and single-layer potentials at a close evaluation point

The nearly singular behavior of (2.9) and (2.10) about y = y is the major cause for error in their evaluation. For this reason, we analyze the contribution made by a small region about y = y . By doing so, we motivate the method that we use to evaluate (2.9) and (2.10). In general, it is more challenging to evaluate the double-layer potential than the single-layer potential because the kernel in (2.9) is more peaked than that in (2.10). However, we show below that the modified double-layer potential given in (2.9) will have a more accurate evaluation than (2.10) due to the subtraction method that has been used. For simplicity, suppose that B is an analytic, closed, and oriented surface. This assumption can be relaxed to just local smoothness about y with appropriate technical modifications. We introduce the parameters s ∈ [0, π] and t ∈ [-π, π] for a rotated spherical coordinate system that parametrizes B such that y /|y | corresponds to the north pole, i.e. y = y(s, t) and y = y(0, •). Details about this rotated coordinate system can be found in A. In terms of this parameterization, (2.9) is given by

u(y -ε n ) = -μ(0, •) + 1 4π π -π π 0 K(s, t; ε) [μ(s, t) -μ(0, •)] sin(s) dsdt, (3.1)
with μ(s, t) = µ(y(s, t)), μ(0, •) = µ(y(0, •)), and

K(s, t; ε) = ñ(s, t) • (y d (s, t) -ε n ) |y d (s, t) -ε n | 3 J(s, t), (3.2)
with ñ(s, t) = n(y(s, t)), y d (s, t) = y(0, •) -y(s, t), and J(s, t) = |y s (s, t) × y t (s, t)|/ sin(s). In the limit as ε → 0, this kernel becomes sharply peaked for

y d = O(ε) leading to K = O(ε -2
). This behavior characterizes this nearly singular integral.

We consider the contribution to u made by a small portion of B about y = y , which we denote by U . To do so, we introduce the parameter δ > 0 and define

U (y ; ε, δ) = 1 4π π -π δ 0 K(s, t; ε) [μ(s, t) -μ(0, •)] sin(s) dsdt. (3.3)
Next, we substitute s = εS in (3.3), and find

U (y ; ε, δ) = 1 4π π -π δ/ε 0 K(εS, t; ε) [μ(εS, t) -μ(0, •)] sin(εS)ε dSdt. (3.4)
Recognizing that ñ(εS, t) = n + O(ε), and y d (εS, t) = -εSy s (0, •) + O(ε 2 ) with y s (0, •) = 0 lying on the plane tangent to B at y , we find by expanding K(εS, t; ε) about ε = 0 that

K(εS, t; ε) = - J(0, •) ε 2 1 (S 2 |y s (0, •)| 2 + 2 ) 3/2 + O(ε -1 ). (3.5)
Substituting this leading-order behavior into (3.4) along with sin(εS) = εS + O(ε 3 ), we obtain

U (y ; ε, δ) = - J(0, •) 4π π -π δ/ε 0 S (S 2 |y s (0, •)| 2 + 2 ) 3/2 + O(ε) [μ(εS, t) -μ(0, •)] dSdt. (3.6)
Since the leading-order asymptotic behavior given in (3.5) is independent of t, we use the 2π-periodicity of μ with respect to t to write

U (y ; ε, δ) = - J(0, •) 4π δ/ε 0 S (S 2 |y s (0, •)| 2 + 2 ) 3/2 + O(ε) π 0 μ(εS, t) + μ(εS, t -π) -2μ(0, •)dt dS.
(3.7) Next, we use the regularity of μ over the north pole to substitute μ(εS, t -π) = μ(-εS, t), and find that

U (y ; ε, δ) = - J(0, •) 4π δ/ε 0 S (S 2 |y s (0, •)| 2 + 2 ) 3/2 + O(ε) π 0 μ(εS, t) + μ(-εS, t) -2μ(0, •)dt dS, = - J(0, •) 4π δ/ε 0 S (S 2 |y s (0, •)| 2 + 2 ) 3/2 + O(ε) π 0 ε 2 S 2 μss (0, •) + O(ε 4 )dt dS, = - J(0, •) 4 δ/ε 0 ε 2 S 3 (S 2 |y s (0, •)| 2 + 2 ) 3/2 1 π π 0 μss (0, •)dt dS + O(ε 2 ). (3.8)
We find that in the asymptotic limit corresponding to 0 < ε δ 1 that

δ/ε 0 ε 2 S 3 (S 2 |y s (0, •)| 2 + 2 ) 3/2 dS = δ ε |y s (0, •)| 3 + O(ε 2 ), (3.9)
(recall that |y s (0, •)| = 0), and in B, we show that

1 π π 0 μss (0, •)dt = 1 2 ∆ S 2 µ(y ), (3.10)
with ∆ S 2 µ(y ) denoting the spherical Laplacian of µ evaluated at y . Thus, we have found that .11) This result gives the leading-order asymptotic behavior for U . The results (3.9) and (3.11) provide the following insights about the close evaluation problem.

U (y ; ε, δ) = -δ ε J(0, •) 8|y s (0, •)| 3 ∆ S 2 µ(y ) + O(ε 2 ). ( 3
• The function that is integrated in (3.9) has a smooth limit to a finite value at the north pole.

In other words, despite the sharply peaked behavior of the kernel, the leading-order asymptotic behavior of the integrand is smooth.

• Working in a rotated spherical coordinate system in which y /|y | serves as the north pole effectively removes the dependence with respect to the azimuthal angle, denoted here by t, through an averaging about the north pole.

• The above procedure provides a way to rewrite the layer potential that avoids the close-evaluation problem.

Similar conclusions can be drawn for the single-layer potential. Using the same parameterization as above, we write (2.10) as

v(y + ε n ) = 1 4π π -π π 0 G(s, t; ε)ρ(s, t) sin(s)dsdt, (3.12)
with ρ(s, t) = ρ(y(s, t)), and

G(s, t; ε) = J(s, t) |y d (s, t) + ε n | . (3.13)
Just as we have done for u, we consider the contribution to v made by a small portion of B about y = y , which we denote by V , defined with the parameter δ as

V (y ; ε, δ) = 1 4π π -π δ 0 G(s, t; ε)ρ(s, t) sin(s)dsdt. (3.14)
Next, we substitute s = εS and obtain

V (y ; ε, δ) = 1 4π π -π δ/ε 0 G(εS, t; ε)ρ(εS, t) sin(εS)εdSdt. (3.15)
By expanding G about ε = 0, we find that that

G(εS, t; ε) = 1 ε J(0, •) (S 2 |y s (0, •)| 2 + 2 ) 1/2 + O(1). (3.16)
Substituting this leading-order behavior into (3.15) along with sin(εS) = εS + O(ε 3 ), and ρ(εS, t) = ρ(0, •) + O(ε), and then integrating that result with respect to t, we obtain

V (y ; ε, δ) = J(0, •)ρ(0, •) 2 δ/ε 0 εS (S 2 |y s (0, •)| 2 + 2 ) 1/2 + O(ε 2 ) dS.
(3.17)

In the asymptotic limit in which 0 < ε δ 1, we find that

V (y ; ε, δ) = δ J(0, •) 2|y s (0, •)| ρ(0, •) + O(ε). (3.18)
This result gives the leading-order asymptotic behavior for V .

In light of the asymptotic analysis given above, we have found a natural procedure to compute the close evaluation of layer potentials that avoids the challenges associated with evaluating nearly singular integrals. There are three, crucial steps in this procedure.

1. Rotate the spherical coordinate system so that y /|y | serves as the north pole where y denotes the boundary point at which the kernel for the layer potential is peaked.

2. Integrate first with respect to the azimuthal angle to average the function about the north pole.

3. Integrate the remaining function with respect to the polar angle.

It is important in the third step to integrate with respect to s rather than z = cos(s), which is often done, because the asymptotic analysis above is explicitly aided by sin(s) term inherent in the spherical integral. These results suggest that following these three steps leads to a natural numerical method to evaluate double-and single-layer potentials close to the boundary. The leading-order asymptotic behaviors given by (3.11) and (3.18) show that this method yields an O(ε 2 ) error for the modified double-layer potential and an O(ε) error for the single-layer potential. The smaller error with the modified double-layer potential is due to the subtraction of µ(y ) in conjunction with the azimuthal averaging step. Additionally (3.11) and (3.18) provide explicit expressions for the limiting behavior of the double-and single-layer potentials, respectively, as s → 0. These results may be used in the computation of the integral with respect to the polar angle in the third step. However, we avoid needing these results in the specific method we give below by choosing an open quadrature rule for s that does not require evaluation at the end points of the interval.

Numerical method for the close evaluation of double-and single-layer potentials

We now give a numerical method for computing the close evaluation of double-and single-layer potentials by following the three-step procedure given above. Suppose we have parameterized B by y = y(θ, ϕ) with θ ∈ [0, π] and ϕ ∈ [-π, π] where y = y(θ , ϕ ). We seek to compute

I 3D (y ) = 1 4π π -π π 0 F (θ, ϕ; θ , ϕ ) sin(θ) dθdϕ, (4.1)
where F is nearly singular at (θ, ϕ) = (θ , ϕ ). For the modified double-layer potential resulting from applying a subtraction method, we have

F (θ, ϕ; θ , ϕ ) = ñ(θ, ϕ) • [y(θ , ϕ ) -y(θ, ϕ) -ε ñ(θ , ϕ )] |y(θ , ϕ ) -y(θ, ϕ) -ε ñ(θ , ϕ )| 3 J(θ, ϕ) [μ(θ, ϕ) -μ(θ , ϕ )] , (4.2)
and for the single-layer potential, we have

F (θ, ϕ; θ , ϕ ) = J(θ, ϕ) |y(θ , ϕ ) -y(θ, ϕ) + ε ñ(θ , ϕ )| ρ(θ, ϕ), (4.3)
where

J(θ, ϕ) = |y θ (θ, ϕ) × y ϕ (θ, ϕ)|/ sin θ.
Motivated by the local analysis leading to (3.11) and (3.18) above, we perform the following threestep procedure to numerically evaluate I 3D (y ).

1. Rotate this integral to another spherical coordinate system whose north pole is aligned with y /|y |. The details of this rotation are given in A, and lead to θ = θ(s, t) and ϕ = ϕ(s, t) with s ∈ [0, π] and t ∈ [-π, π] where θ(0, •) = θ and ϕ(0, •) = ϕ . Upon applying this rotation, we obtain the spherical integral

I 3D (y ) = 1 4π π -π π 0 F (s, t) sin(s) dsdt, (4.4)
with F (s, t) = F (θ(s, t), ϕ(s, t); θ , ϕ ).

2. Use the 2N -point Periodic Trapezoid Rule to compute

1 2π π -π F (s, t)dt = F (s) ≈ F N (s) = 1 2N 2N j=1 F (s, t j ), (4.5)
with t j = -π + π(j -1)/N for j = 1, . . . , 2N .

3. Use the N -point Gauss-Legendre quadrature rule to compute the integral in s. Let z i ∈ (-1, 1) for i = 1, • • • , N denote the Gauss-Legendre quadrature points with corresponding weights w i . We map the quadrature points according to s i = π(z i + 1)/2 to the interval (0, π) and compute the numerical approximation,

I 3D ≈ I N 3D = π 4 N i=1 F N (s i ) sin(s i )w i . (4.6)
In the first step of the procedure, the coordinate system is rotated so that (θ , ϕ ) corresponds to the north pole of the spherical coordinate system. This rotation sets up the second step, which computes the azimuthal average of the function with respect to t about the north pole. The factor of 1/(2π) is incorporated into this step to explicitly obtain an averaging operation yielding F N (s), which is a smooth function in the polar angle, s. For the third step, we have chosen to use the Gauss-Legendre quadrature rule for which the quadrature points satisfy 0 < s 1 < s 2 < • • • < s N < π. Consequently, we do not need to explicitly evaluate F N (0). Alternatively, we could have chosen a closed rule that includes the end points. For that case, we would use (3.11) and (3.18) to explicitly prescribe the needed values in the limit as s → 0.

This numerical method is based on the product Gaussian quadrature rule for spherical integrals [START_REF] Atkinson | Numerical integration on the sphere[END_REF]. That quadrature rule uses a 2N -point Periodic Trapezoid Rule in the azimuthal angle and N -point Gauss-Legendre quadrature rule in the cosine of the polar angle, z = cos(s) ∈ [-1, 1]. However, the results of the local analysis above show that integrating with respect to the polar angle, s, is better since the factor of sin(s) plays an important role in the asymptotic behavior. Thus, rather than transforming (4.1) by substituting z = cos(s), we map the Gauss-Legendre quadrature rule from [-1, 1] to [0, π] using s i = π(z i + 1)/2 for i = 1, • • • , N . Consequently, the mapped quadrature weights are πw j /2, which is accounted for in (4.6) by the π/4 factor. We show in the numerical results presented in Section 6 that this subtle and seemingly small change to the product Gaussian quadrature rule is significant and yields substantially better results.

Provided that the densities, µ and ρ, for the double-and single-layer potentials, respectively, are known to very high accuracy, the main cause of numerical error for the close evaluation problem is due to the sharp peak of the kernels about y = y . For that case, the leading-order asymptotic expansions given by (3.11) and (3.18) indicate that this numerical method yields an O(ε 2 ) error for the modified double-layer potential and an O(ε) error for the single-layer potential, respectively. These asymptotic error estimates are valid only when the value of N used in this numerical method is sufficiently large to resolve the density and the boundary. Otherwise, errors introduced by inaccuracy of the density and/or under-resolution of the boundary will dominate over these asymptotic error estimates.

Numerical solution of the boundary integral equations

The numerical method described above can also be used to solve the boundary integral equations for the densities, (2.3) and (2.6). We describe its use in the context of the Galerkin method [START_REF] Atkinson | The numerical solution Laplace's equation in three dimensions[END_REF][START_REF] Atkinson | Algorithm 629: An integral equation program for Laplace's equation in three dimensions[END_REF][START_REF] Atkinson | A survey of boundary integral equation methods for the numerical solution of Laplace's equation in three dimensions[END_REF][START_REF] Atkinson | The Numerical Solution of Integral Equations of the Second Kind[END_REF]]. As we have done above, we suppose that we have parameterized B by y = y(θ, ϕ) with θ ∈ [0, π] and ϕ ∈ [-π, π].

To compute the density, µ, for the double-layer potential, we write (2.3) compactly as

- 1 2 µ + K D [µ] = f. (5.1)
We introduce the approximation for µ,

µ(y(θ, ϕ)) = μ(θ, ϕ) ≈ N -1 n=0 n m=-n Y nm (θ, ϕ)μ nm , (5.2)
with {Y nm (θ, ϕ)} denoting the orthonormal set of spherical harmonics. Note that N in (5.2) corresponds also to the order of the quadrature rule detailed in Section 4 so that the order of this approximation is commensurate with the order of the quadrature rule to be used. Substituting (5.2) into (5.1) and taking the inner product with Y n m (θ , ϕ ), we obtain the Galerkin equations,

- 1 2 μn m + N -1 n=0 n m=-n Y n m , K D [Y nm ] μnm = Y n m , f , (5.3)
with

K D [Y nm ](θ , ϕ ) = 1 4π π -π π 0 ñ(θ, ϕ) • (y(θ , ϕ ) -y(θ, ϕ)) |y(θ , ϕ ) -y(θ, ϕ)| 3 J(θ, ϕ) sin(θ) Y nm (θ, ϕ)dθdϕ. (5.4)
We construct the N 2 × N 2 linear system for the unknown coefficients, μn m resulting from (5.3) evaluated for n = 0, • • • , N -1 with corresponding values of m . To evaluate (5.4), we use (2.8) to rewrite it as

K D [Y nm ](θ , ϕ ) = - 1 2 Y nm (θ , ϕ ) + 1 4π π -π π 0 ñ(θ, ϕ) • (y(θ , ϕ ) -y(θ, ϕ)) |y(θ , ϕ ) -y(θ, ϕ)| 3 J(θ, ϕ) sin(θ) [Y nm (θ, ϕ) -Y nm (θ , ϕ )] dθdϕ, (5.5)
and evaluate the integral in this expression using the numerical method described above. To compute the inner products, Y n m , K D [Y nm ] and Y n m , f , we use the product Gaussian quadrature rule for spherical integrals [START_REF] Atkinson | Numerical integration on the sphere[END_REF] with N Gauss-Legendre quadrature points and 2N Periodic Trapezoid Rule points.

To compute the density ρ for the single-layer potential, we write (2.6) compactly as

1 2 ρ -K S [ρ] = g. (5.6)
We introduce the approximation for ρ,

ρ(y(θ, ϕ)) = ρ(θ, ϕ) ≈ N -1 n=0 n m=-n Y nm (θ, ϕ)ρ nm .
(5.7)

Performing the same calculations as above, we find that ρnm satisfies the Galerkin equations,

1 2 ρn m - N -1 n=0 n m=-n K † S [Y n m ], Y nm ρnm = Y n m , g , (5.8)
where K † S denotes the adjoint of K S defined according to

K † S [Y n m ](θ, ϕ) = 1 4π π -π π 0 ñ(θ , ϕ ) • (y(θ , ϕ ) -y(θ, ϕ)) |y(θ , ϕ ) -y(θ, ϕ)| 3 J(θ , ϕ ) sin(θ ) Y n m (θ , ϕ )dθ dϕ .
(5.9)

Recall that ñ(θ , ϕ ) is the unit normal pointing out of D. For this integral, we use (2.8) to write

K † S [Y n m ](θ, ϕ) = - 1 2 Y n m (θ, ϕ) + 1 4π π -π π 0 ñ(θ , ϕ ) • (y(θ , ϕ ) -y(θ, ϕ)) |y(θ , ϕ ) -y(θ, ϕ)| 3 J(θ , ϕ ) sin(θ ) [Y n m (θ , ϕ ) -Y n m (θ, ϕ)] dθ dϕ , (5.10)
which we then compute using the numerical method described above. We use the product Gaussian quadrature rule with N Gauss-Legendre quadrature points and 2N Periodic Trapezoid Rule points to compute the remaining inner products, K † S Y n m , Y nm and Y n m , g , needed to form the N 2 × N 2 linear system for the unknown coefficients, ρn m for n = 0, • • • , N -1 and corresponding values of m .

Numerical results

We present numerical results for computing the close evaluation of double-and single-layer potentials using the method detailed in Section 4. For the double-layer potential examples discussed in Section 6.1 below, we have computed the density, μ, by solving boundary integral equation (5.1) using the Galerkin method discussed in Section 5. Then, we compute the double-layer potential using that density for two ellipsoids and study the error. In Section 6.2 below, we compute the representation formula for a harmonic function given as the difference between the double-and single-layer potentials using densities that are known exactly. We study the error in computing this representation formula in two different domains and evaluate those results.

Double-layer potential

We present results in computing the solution of the interior Dirichlet problem (2.1) for two different ellipsoids. For all of these examples, we set the harmonic function,

u ex (x) = |x -x 0 | -1 , x ∈ D, (6.1)
with x 0 ∈ E to be the exact solution. Therefore, we prescribe

f (y) = |y -x 0 | -1 , y ∈ B, (6.2)
as the Dirichlet boundary data.

With the boundary data given in (6.2), we first use the Galerkin method described in Section 5 to solve for μnm , the expansion coefficients for the density introduced in (5.2). We then use (5.2) to evaluate the density, μ, that is needed for computing the double-layer potential. To test our method for close evaluation points, we evaluate the double-layer potential at the points x = y -εn (we have set = 1) for various points y on the boundary of the domain and for various ε. The error is determined by computing the absolute value of the difference between the double-layer potential result and the exact harmonic function, u ex (y -εn ) = |y -εn -x 0 | -1 . We show error results for computing u(x) over a range of ε values spanning several orders of magnitude using four different methods listed below using the kernel, K defined in (3.2). For each of these methods, we have applied the rotation of the coordinate system described as the first step in the procedure detailed in Section 4. Methods 1 and 2 compute the double-layer potential given by (2.2) directly and integrate with respect to z = cos(s) and s, respectively. Methods 3 and 4 use the subtraction method, so they compute the modified double-layer potential given in (2.9). Method 3 integrates with respect to z = cos(s) and Method 4 integrates with respect to s. Method 4 is the one we have proposed here to compute the close-evaluation of the double-layer potential.

In the examples that follow, we consider two different ellipsoids used by Atkinson [START_REF] Atkinson | The numerical solution Laplace's equation in three dimensions[END_REF]. Let (x 1 , x 2 , x 3 ) denote an ordered triple in a Cartesian coordinate system. The boundaries for these ellipsoids are given by y(θ, ϕ) = (α sin θ cos ϕ, β sin θ sin ϕ, γ cos θ), θ ∈ [0, π], ϕ ∈ [-π, π]. (6.3)

For the first example, we set α = 1, β = 1.5, and γ = 2, and for the second example, we set α = 1, β = 2, and γ = 5. These two ellipsoids, which we call Ellipsoids 1 and 2, are shown as the left and right plots in Fig. 2, respectively. For both of these examples, we have set x 0 = (5, 4, 3) ∈ E in (6.1) and (6.2). In Fig. 2 we also identify two boundary points, labeled A and B, lying on the x 1 x 3 -plane and the x 1 x 2 -plane, respectively, which we will refer to below when we discuss the error results. The coordinates for these points are A = (-1.16, -0.64, 0) and B = (-0.72, 0, 1.75) for Ellipsoid 1, and A = (-0.74, -0.93, 0) and B = (-1.81, 0, 2.13) for Ellipsoid 2.

We have used the Galerkin method with N = 48 to compute the density, μ. For Ellipsoid 1 shown in the left plot of Fig. 2, we find that this order resolves μ with an estimated truncation error on the order of 10 -15 . Ellipsoid 2 shown in the right plot of Fig. 2 is more elongated along the x 3 axis, so it deviates farther from a sphere than Ellipsoid 1. As a result, we find that the Galerkin method with N = 48 yields a larger estimated truncation error for μ on the order of 10 -12 . We have subsequently used N = 48 in the quadrature rules for the four methods given above.

Results for Ellipsoid 1

In the top row of Fig. 3, we plot log 10 of the error made using Method 1. The left two plots show the error over the vertical x 1 x 3 -plane and a close-up of the small region near the boundary point, A = (-1.16, -0.64, 0). The right two plots show the error over the horizontal x 1 x 2 -plane and a close-up of the small region near the boundary point, B = (-0.72, 0, 1.75). The bottom row of Fig. 3 shows the same plots for Method 4.

These results show that using Method 1 leads to an O(1) error close to the boundary which is symptomatic of the close evaluation problem. In contrast, the results using Method 4 show an error that is several orders of magnitude smaller. To examine this error more closely, we show in Fig. 4 errors using Methods 1, 2, 3, and 4 evaluated along the normal at point A = (-1.16, -0.64, 0) on the boundary and the x 1 x 3 -plane, and at point B = (-0.72, 0, 1.75) on the boundary and the x 1 x 2 -plane. In this way, we clearly observe that the subtraction method used in Method 3 and Method 4 leads to substantially smaller errors compared with those made using Method 1 and Method 2. The error for Methods 1, 2, and 3 all appear to asymptote to a fixed value as ε → 0. In contrast, we see that the error for Method 4 decreases as ε → 0. To study the asymptotic behavior of the error for Method 4, we compute the best linear fit to the error (see dashed black lines in Fig. 4) for the four smallest ε values. These linear fits have slopes of 1.73 for point A and 1.83 for point B consistent with the O(ε 2 ) error estimate from the asymptotic analysis in (3.11).

Results for Ellipsoid 2

The error results for Ellipsoid 2 plotted in Figs. 5 and6 are of the same format as those for Ellipsoid 1 in Figs. 3 and4, respectively. The close-ups in Fig. 5 show small regions near the boundary points A = (-0.74, -0.93, 0) and B = (-1.81, 0, 2.13).

Figure 5 shows that the error for Method 4 is again several orders of magnitude smaller than that for Method 1. However, the magnitude of the errors appears to be larger for both methods than the previous example. This is due to the larger truncation error associated with the Galerkin method. In Fig. 6, we see that Method 4 clearly has the smallest error. The dashed black lines show the best linear fit to the errors for the four smallest ε values. These linear fits have slopes of 1.83 for point A and 1.92 for point B consistent with O(ε 2 ) error estimate in (3.11).

Summary of double-layer potential results

For both Ellipsoids 1 and 2, we have observed that Method 4 provides the most accurate evaluation of the double-layer potential for close evaluation points. Method 3 is the next most accurate. However, the errors for Methods 1, 2, and 3 all appear to asymptote to a finite value rather than decay to zero as the error for Method 4 does. For this reason, we conclude that Method 4 is an effective method for the close evaluation of the double-layer potential.

The errors for Method 4 plotted in Figs. 4 and6 show a characteristic behavior. For very small values of ε corresponding to points very close to the boundary, Method 4 is highly accurate with an error approaching machine precision. For intermediate values of ε, the error obtained with Method 4 increases as the value of ε starts to coincide with the spacing of the Gauss-Legendre quadrature points in the polar angle s. For bigger values of ε, Methods 1, 2, 3, and 4 have errors of a similar order of magnitude because the evaluation point is far enough from the boundary that the double-layer potential is no longer nearly singular.

The representation formula

When the truncation error for the Galerkin method is not sufficiently small, it dominates over the asymptotic error for the close evaluation of the double-layer potential given by (3.11) and the singlelayer potential given by (3.18). For the exterior Neumann problem, (2.4), boundary data are in terms of the gradient of a harmonic function. Consequently, those boundary data tend to be less smooth than those for the interior Dirichlet problem. Because of this, the Galerkin method to solve boundary integral equation (5.6) for a boundary that deviates significantly from a sphere requires N in (5.7) to be sufficiently large that the truncation error does not interfere with the O(ε) asymptotic error estimate for the single-layer potential. It is in this way that the Galerkin method is a limitation and another high-order method can be used to compute the densities for the double-and single-layer potentials for more general domains [START_REF] Bruno | A fast, high-order algorithm for the solution of surface scattering problems: basic implementation, tests, and applications[END_REF].

The focus of this paper is in the close evaluation problem. Thus, we consider an alternate example to test the asymptotic error estimates for the double-and single-layer potentials rather than consider another high-order method to solve weakly singular boundary integral equations for the densities. Consider the representation formula for a harmonic function in domain D with boundary B,

u(x) = 1 4π B 1 |x -y| n(y) • ∇u(y) - n(y) • (x -y) |x -y| 3 u(y) dσ y , x ∈ D. (6.4)
Given a prescribed harmonic function u and its normal derivative on B, (6.4) gives a formula for computing u(x) for any x ∈ D in terms of the difference between the single-and double-layer potentials, where the prescribed boundary data can be interpreted as the corresponding densities, ρ and µ. For this case there is no boundary integral equation to solve since the densities are known exactly.

For the following examples, we have set the harmonic function

u ex (x 1 , x 2 , x 3 ) = e x3 (sin x 1 + sin x 2 ), (6.5)
following Beale et al. [START_REF] Beale | A simple method for computing singular or nearly singular integrals on closed surfaces[END_REF], to be the exact solution. We set (6.5) on the boundary to be µ and its normal derivative on the boundary to be ρ, which are used in the double-and single-layer potentials, respectively, to evaluate (6.4). For the following examples, we use the four methods described below. Here, K is the kernel defined in (3.2) and G is the kernel defined in (3.13).

Method 1: Compute u(y -εn ) = 1 4π π -π 1 -1 G(cos -1 (z), t; ε)ρ(cos -1 (z), t) -K(cos -1 (z), t; ε)μ(cos -1 (z), t)
Method 3: Compute u(y -εn ) = μ(0, •)+ 1 4π π -π 1 -1 G(s, t; ε)ρ(cos -1 (z), t) -K(cos -1 (z), t; ε) μ(cos -1 (z), t) -μ(0, •) dzdt using an N -point Gauss-Legendre quadrature rule for z ∈ [-1, 1], and a 2N Periodic Trapezoid Rule for t ∈ [-π, π]. Method 4: Compute u(y -εn ) = μ(0, •) + 1 4π π -π π 0 [G(s, t; ε)ρ(s, t) -K(s, t; ε) [μ(s, t) -μ(0, •)]] sin(s)dsdt
Those four methods are a combination of the four methods used in Section 6.1 for the double-layer potential, and the corresponding evaluation of the single-layer potential. Methods 1 and 2 do not use the subtraction method for the double-layer potential, and Methods 3 and 4 use the subtraction method. Methods 1 and 3 integrate with respect to z = cos(s), and Methods 2 and 4 integrate with respect to s. According to the asymptotic error estimates given by (3.11) and (3.18), we expect Method 4 to exhibit an O(ε) error in evaluating (6.4). In the examples that follow, we consider two domains used by Atkinson [START_REF] Atkinson | The numerical solution Laplace's equation in three dimensions[END_REF][START_REF] Atkinson | Algorithm 629: An integral equation program for Laplace's equation in three dimensions[END_REF]. The boundaries for these two domains are given by y(θ, ϕ) = R(θ)(sin θ cos ϕ, 2 sin θ sin ϕ, cos θ), θ ∈ [0, π], ϕ ∈ [-π, π].

(6.6)

For the first example, which resembles a peanut, we have set R(θ) = cos 2θ + 1.1 -sin 2 2θ. (6.7)

For the second example, which resembles an inverted mushroom cap, we have set

R(θ) = 2 - 1 1 + 100(1 -cos θ) . (6.8)
The peanut and inverted mushroom cap domains are shown as the left and right plots in Fig. 7, respectively. In Fig. 7, we identify two boundary points, labeled A and B, lying on the x 1 x 3 -plane and the x 1 x 2 -plane, respectively. For the peanut, A = (-0.16, -0.21, 0) and B = (-1.04, 0, 0.94). For the inverted mushroom cap, A = (-3.07, -1.27, 0) and B = (-1.55, 0, 1.41). We refer to these points below when we discuss the error results. Since we use exact densities, we are more easily able to use highly resolved computations for these two examples. For the peanut domain, we have set N = 120 and for the inverted mushroom cap, we have set N = 80.

Results for peanut-shaped domain

The error results plotted in Figs. 8 and9 are the same format as those in Figs. 3 and4, respectively. Just like those for the double-layer potential, the results in Fig. 8 clearly show that the error for Method 4 is several orders of magnitude smaller than that for Method 1. Moreover, Fig. 9 shows that Methods 3 and 4 have significantly smaller errors than Methods 1 and 2 demonstrating that using the modified double-layer potential resulting from applying the subtraction method is effective. These errors show that Methods 1, 2, and 3 asymptote to a fixed value as ε → 0. In contrast, the error for Method 4 decreases as ε → 0. To study the asymptotic behavior of the error for Method 4, we compute the best linear fit to the error (see dashed black lines in Fig. 9) for the four smallest ε values along the normal for two chosen points, A and B. These linear fits have slopes of 1.00 for point A and 1.01 for point B consistent with the expected O(ε) error estimate.

Results for inverted mushroom cap domain

The error results plotted in Figs. 10 and 11 are the same format as those in Figs. 8 and9, respectively. Fig. 10 also shows that Method 4 clearly outperforms Method 1. Fig. 11 shows that the subtraction method used for Methods 3 and 4 yields much better results compared with those for Methods 1 and 2. Again, Methods 1, 2, and 3 all appear to asymptote to a fixed value as ε → 0, where Method 4 decreases. The dashed black lines in Fig. 11 show the best linear fits to the error for Method 4 for the four smallest ε values. The slopes of these linear fits are 1.01 for point A and 1.00 for point B consistent with the expected O(ε) error estimate given by (3.18).

Summary of representation formula results

For both the peanut and inverted mushroom cap domains, we have observed that Method 4 provides the most accurate evaluation of the representation formula for close evaluation points. Method 3 is the next most accurate due to the use of the subtraction method. The errors of Methods 1, 2, and 3 all asymptote to a finite value as ε → 0, whereas the error for Method 4 decays to zero as O(ε). For this reason, we conclude that Method 4 is effective for the close evaluation of the representation formula. The error for Method 4 demonstrates the same characteristic behavior discussed for the double-layer

Conclusions

We have presented a simple and effective numerical method for the close evaluation of double-and single-layer potentials in three dimensions. The close evaluation of these layer potentials are challenging to numerically compute because they are nearly singular integrals. Based on a local analysis of these layer potentials about the point at which their kernel is sharply peaked, we have identified a natural way to compute those nearly singular integrals. This method consists of three steps:

1. Rotate the coordinate system so that the point at which the integrand is sharply peaked is aligned with the north pole.

2. Integrate first with respect to the azimuthal angle using the Periodic Trapezoid Rule.

3. Integrate the remaining function with respect to the polar angle using the Gauss-Legendre Quadrature Rule mapped from [-1, 1] to [0, π].

A key point to Step 3 is that we do not integrate with respect to the cosine of the polar angle because the factor of sine of the polar angle (the natural Jacobian for the spherical coordinate system) is important in effectively computing these integrals. Numerical results presented here show that the last step is essential for accurately computing these integrals. Moreover, the numerical results confirm the asymptotic analysis results: the error for the modified double-layer potential due to the subtraction method decreases quadratically with respect to the distance from the boundary and the error for the single-layer potential decreases linearly with respect to the distance from the boundary. This method is effective because integration with respect to azimuth in the rotated coordinate system is a natural averaging operation leading to a smooth function in the polar angle. This averaging operation is a feature of the three-dimensional problem. It is not available for two-dimensional problems. The results from the local analysis given by (3.11) and (3.18) show that this azimuthal averaging yields a function in the polar angle that has a finite limit at the north pole of the rotated coordinate system. Rather than use these results in a quadrature rule for the polar angle, we have chosen to use the Gauss-Legendre Quadrature Rule, which is an open rule, so that the quadrature points never include the north and south poles. Therefore, these limiting values are not explicitly needed.

The local analysis used here can be extended to compute the full asymptotic expansions for the close evaluation of double-and single-layer potentials. These asymptotic expansions can be used to find approximations to a desired order of accuracy. Details of these asymptotic expansions can be found in [START_REF] Carvalho | Asymptotic approxiamtion for the close evaluation of doublelayer potentials[END_REF]. These asymptotic expansions may be useful if one requires more accuracy for the close evaluation of layer potentials. Nonetheless, we find that O(ε 2 ) asymptotic error estimate for the modified double-layer potential and the O(ε) asymptotic error estimate for the single-layer potential to be quite accurate. An asymptotic expansion for the single-layer potential may be useful to obtain an asymptotic error that matches that of the modified double-layer potential, for example.

We have assumed that the boundary is closed, oriented, and analytic in this paper, but as long as the surface near the point on the boundary at which the integrand is singular is relatively smooth so that there exists a natural averaging operation, one can apply this approach. Doing so may require consideration of different quadrature rules for that specific context.

Finally, the analysis shown here can be extended to other problems. In particular, this method is broadly applicable to weakly singular and nearly singular integrals over two-dimensional surfaces. These weakly singular and nearly singular integrals may correspond to solutions of other elliptic partial differential equations including Helmholtz's equation and Stokes' equations.

A Rotations on the sphere

We give the explicit rotation formulas over the sphere used throughout this paper. Consider y, y ∈ S 2 , the unit sphere. We introduce the parameters θ ∈ [0, π] and ϕ ∈ [-π, π] and write y = y(θ, ϕ) = sin θ cos ϕî + sin θ sin ϕ + cos θ k.

(A.1)

The specific parameters θ and ϕ are set so that y = y(θ , ϕ ). We would like to work in the rotated, uvw-coordinate system in which Let us rewrite (A.4) compactly as y(θ, ϕ) = R(θ , ϕ )y(s, t) with R(θ , ϕ ) denoting the 3 × 3 matrix given above. It is a rotation matrix. Hence, it is orthogonal. We now seek to write θ = θ(s, t) and ϕ = ϕ(s, t). To do so, we introduce ξ(s, t; θ , ϕ ) = cos θ cos ϕ sin s cos t -sin ϕ sin s sin t + sin θ cos ϕ cos s, (A.5) η(s, t; θ , ϕ ) = cos θ sin ϕ sin s cos t + cos ϕ sin s sin t + sin θ sin ϕ cos s, (A.6)

ζ(s, t; θ , ϕ ) = -sin θ sin s cos t + cos θ cos s.

B Spherical Laplacian

In this Appendix, we establish the result given in (3.10). We first seek an expression for ∂

Figure 1 :

 1 Figure 1: The point y on the boundary surface denoted by B, indicated by the outward, unit normal, n(y ). Even though the fundamental solution for Laplace's equation and its normal derivative are singular at y , integration over B includes a natural averaging operation about y indicated by the dashed closed curve. This averaging operation is important for computing the close evaluation of double-and single-layer potentials in three dimensions.

  Consider a simply connected, open set denoted by D ⊂ R 3 with boundary, B, and let D = D ∪ B. For some given smooth data f , we write the function u ∈ C 2 (D) ∩ C 1 ( D) satisfying the interior Dirichlet problem, ∆u = 0 in D, u = f on B,

  ) • (y -y) |y -y| 3 µ(y)dσ y = f (y ), y ∈ B. (2.3) Let E = R 3 \ D, and Ē = E ∪ B. Let n still denote the unit normal pointing outward from D. For some given smooth data g, we write the function v ∈ C 2 (E) ∩ C 1 ( Ē) satisfying the exterior Neumann problem, ∆v = 0 in E,

Method 1 :4π π -π 1 - 1 K

 111 Computeu(y -εn ) = 1 (cos -1 (z), t; ε)μ(cos -1 (z), t)dzdt using an N -point Gauss-Legendre quadrature rule for z ∈ [-1, 1] and a 2N Periodic Trapezoid Rule for t ∈ [-π, π].

Method 3 : 1 K

 31 t; ε)μ(s, t) sin(s)dsdt using an N -point Gauss-Legendre quadrature rule for s ∈ [0, π] and a 2N Periodic Trapezoid Rule for t ∈[-π, π]. Compute u(y -εn ) = -μ(0, •) (cos -1 (z), t; ε) μ(cos -1 (z), t) -μ(0, •) dzdtusing an N -point Gauss-Legendre quadrature rule for z ∈ [-1, 1] and a 2N Periodic Trapezoid Rule for t ∈ [-π, π]. Method 4: Compute u(y -εn ) = -μ(0, •) , t; ε) [μ(s, t) -μ(0, •)] sin(s)dsdt using an N -point Gauss-Legendre quadrature rule for s ∈ [0, π] and a 2N Periodic Trapezoid Rule for t ∈ [-π, π].

Figure 2 :

 2 Figure 2: (Left) Boundary surface of the ellipsoid given by (6.3) with α = 1, β = 1.5, and γ = 2 which we call Ellipsoid 1. (Right) Boundary surface of the ellipsoid given by (6.3) with α = 1, β = 2, and γ = 5, which we call Ellipsoid 2. In each figure, the two black curves indicate where the boundary of the ellipsoid intersects with the vertical x 1 x 3 -plane, and the horizontal x 1 x 2 -plane, and the black arrows give the unit outward normal vectors on those curves. The two boundary points labeled A and B lying on the x 1 x 3 -plane and x 1 x 2 -plane, respectively, give the locations where we investigate the error made by four different methods. The coordinates for these points are A = (-1.16, -0.64, 0) and B = (-0.72, 0, 1.75) in the left plot, and A = (-0.74, -0.93, 0) and B = (-1.81, 0, 2.13) in the right plot.

Figure 3 :

 3 Figure 3: The errors made by Method 1 (top row) and Method 4 (bottom row) for Ellipsoid 1. The first and second columns show log 10 of the errors made by these two methods in the x 1 x 3 -plane, with the plots in the second column showing a close-up in the vicinity of point A = (-1.16, -0.64, 0). The third and fourth columns show log 10 of the errors made by these two methods in the x 1 x 2 -plane, with the plots in the fourth column showing a close-up in the vicinity of point B = (-0.72, 0, 1.75).

Figure 4 :

 4 Figure 4: Comparison of errors made by Methods 1 -4 for Ellipsoid 1 for evaluation points along the normal at distance ε from A = (-1.16, -0.64, 0) on the boundary and in the x 1 x 3 -plane (left plot) and for evaluation points along the normal at distance ε from point B = (-0.72, 0, 1.75) on the boundary and in the x 1 x 2 -plane (right plot). The dashed black lines show the best linear fit through the Method 4 errors for the four smallest ε values. Both fits are consistent with the O(ε 2 ) asymptotic error estimate given by (3.11).

Figure 5 :

 5 Figure 5: The errors made by Method 1 (top row) and Method 4 (bottom row) for Ellipsoid 2. The first and second columns show log 10 of the errors made by these two methods in the x 1 x 3 -plane with the plots in the second column showing a close-up in the vicinity of point A = (-0.74, -0.93, 0). The third and fourth columns show log 10 of the errors made by these two methods in the x 1 x 2 -plane with the plots in the fourth column showing a close-up in the vicinity of point B = (-1.81, 0, 2.13).

Figure 6 :

 6 Figure 6: Comparison of errors made by Methods 1 -4 for Ellipsoid 2 for evaluation points along the normal at distance ε from point A = (-0.74, -0.93, 0) on the boundary and in the x 1 x 3 -plane (left plot) and for evaluation points along the normal at distance ε from point B = (-1.81, 0, 2.13) on the boundary and in the x 1 x 2 -plane (right plot). The dashed, black lines show the best linear fit through the Method 4 errors for the four smallest ε values. Both fits are consistent with the O(ε 2 ) asymptotic error estimate given by (3.11).

  dzdt using an N -point Gauss-Legendre quadrature rule for z ∈ [-1, 1], and a 2N Periodic Trapezoid Rule for t ∈ [-π, π].

  s, t; ε)ρ(s, t) -K(s, t; ε)μ(s, t)] sin(s)dsdt using an N -point Gauss-Legendre quadrature rule for s ∈ [0, π], and a 2N Periodic Trapezoid Rule for t ∈ [-π, π].

  using an N -point Gauss-Legendre quadrature rule for s ∈ [0, π], and a 2N Periodic Trapezoid Rule for t ∈ [-π, π].

Figure 7 :

 7 Figure7: (Left) Boundary surface of the peanut-shaped domain given by (6.6) with R(θ) given by (6.7). (Right) Boundary surface of the inverted mushroom cap domain given by (6.6) with R(θ) given by (6.8). In both domains, the two black curves indicate where the boundary intersects with the vertical x 1 x 3 -plane, and the horizontal x 1 x 2 -plane, and the black arrows give the unit outward normal vectors on those curves. The two boundary points labeled A and B lying on the x 1 x 3 -plane and x 1 x 2plane, respectively, give the locations where we investigate the error made by four different methods. For the peanut-shaped domain, A = (-0.16, -0.21, 0) and B = (-1.04, 0, 0.94), and for the inverted mushroom cap, A = (-3.07, -1.27, 0) and B = (-1.55, 0, 1.41).

Figure 8 :

 8 Figure 8: The errors made by Method 1 (top row) and Method 4 (bottom row) for the peanut-shaped domain shown on the left plot of Fig. 8. The first and second columns show log 10 of the errors made by these two methods in the x 1 x 3 -plane with the plots in the second column showing a close-up in the vicinity of point A = (-0.16, -0.21, 0). The third and fourth columns show log 10 of the errors made by these two methods in the x 1 x 2 -plane with the plots in the fourth column showing a close-up in the vicinity of point B = (-1.04, 0, 0.94).

Figure 9 :Figure 10 :

 910 Figure 9: Comparison of errors made by Methods 1 -4 for the peanut-shaped domain for evaluation points along the normal at distance ε from point A = (-0.16, -0.21, 0) on the boundary and in the x 1 x 3 -plane (left plot) and for evaluation points along the normal at distance ε from point B = (-1.04, 0, 0.94).on the boundary and in the x 1 x 2 -plane (right plot). The dashed, black lines show the best linear fit through the Method 4 errors for the four smallest ε values. For both fits, the slopes are consistent with the O(ε) asymptotic error estimate.

Figure 11 :

 11 Figure 11: Comparison of errors made by Methods 1 -4 for the inverted mushroom cap for evaluation points along the normal at distance ε from point A = (-3.07, -1.27, 0) on the boundary and in the x 1 x 3 -plane (left plot) and for evaluation points along the normal at distance ε from point B = (-1.55, 0, 1.41) on the boundary and in the x 1 x 2 -plane (right plot). The dashed, black lines show the best linear fit through the Method 4 errors for the four smallest ε values. For both fits, the slopes are consistent with the O(ε) asymptotic error estimate given by (3.18).

  û = cos θ cos ϕ î + cos θ sin ϕ  -sin θ k, v = -sin ϕ î + cos ϕ , ŵ = sin θ cos ϕ î + sin θ sin ϕ  + cos θ k. (A.2)In this rotated system we have ŵ = y . For this rotated coordinate system, we introduce the parameters s ∈ [0, π] and t ∈ [-π, π] such that y = y(s, t) = sin s cos t û + sin s sin t v + cos s ŵ.(A.3)It follows that ŵ = y = y(0, •). This corresponds to setting y to be the north pole. By equating (A.1) and (A.3) and substituting (A.2) into that result, we obtain   sin θ cos ϕ sin θ sin ϕ cos θ ϕ -sin ϕ sin θ cos ϕ cos θ sin ϕ cos ϕ sin θ sin ϕ -

  2 s [•]| s=0 in terms of θ and ϕ. By the chain rule, we find that Note that at s = 0, we have θ = θ. Substituting (B.2) -(B.5) into (B.1) and replacing θ by θ, we obtain

	∂ 2 ∂s 2 [•]	s=0	= cos 2 t	∂ 2 ∂θ 2 + sin 2 t	1 sin 2 θ	∂ 2 ∂ϕ 2 + 2 cos t sin t	1 sin θ	∂ 2 ∂θ∂ϕ	+ sin 2 t	cos θ sin θ	∂ ∂θ	-sin 2t	cos θ sin 2 θ	∂ ∂ϕ	,
																			(B.6)
	from which it follows that
				1 π	0	π	∂ 2 ∂s 2 [•]	s=0	dt =	1 2	∂ 2 ∂θ 2 +	cos θ sin θ	∂ ∂θ	+	1 sin 2 θ	∂ 2 ∂ϕ 2 =	1 2	∆ S 2 ,	(B.7)
	which establishes the result.
	∂ 2 ∂s 2 [•]	s=0	=		∂θ ∂s	2 ∂ 2 ∂θ 2 +	∂ϕ ∂s	2 ∂ 2 ∂ϕ 2 + 2	∂θ ∂s	∂ϕ ∂s	∂ 2 ∂θ∂ϕ	+	∂ 2 θ ∂s 2	∂ ∂θ	+	∂ 2 ϕ ∂s 2	∂ ∂ϕ s=0	.	(B.1)
	Using θ defined in (A.8) and ϕ defined in (A.9), we find that
										∂θ(s, t) ∂s s=0	= cos t,	(B.2)
										∂ 2 θ(s, t) ∂s 2	s=0	=	cos θ sin θ	sin 2 t,	(B.3)
										∂ϕ(s, t) ∂s	s=0	=	sin t sin θ	,	(B.4)
										∂ 2 ϕ(s, t) ∂s 2	s=0	= -	cos θ sin 2 θ	sin 2t.	(B.5)
															21

Acknowledgments

This research was supported by National Science Foundation Grant: DMS-1819052. A. D. Kim is also supported by the Air Force Office of Scientific Research Grant: FA9550-17-1-0238.


Andrea Pimpinella

Alessandro E C Redondi

Matteo Cesana

Load Balancing and Performance Optimization in wM-Bus Smart Meter Networks

Smart meters networks are rapidly becoming a reality in many developed countries. In this paper, we focus on the optimization of a network of smart meters operated by the wM-Bus protocol, which is the de-facto metering standard in Europe. In such a scenario, data concentrators receive data from smart meters using the wM-Bus protocol and relay it to a central server using a legacy mobile cellular backhauling technology such as GSM/GPRS. Due to the massive amount of data produced by meters installed in urban scenarios and the association-less nature of the wM-Bus protocol, data concentrators may be overloaded with many duplicate measurement packets, causing capacity problems on the backhauling links and computational overload at the central server. To solve these issues, we propose a datadriven optimization framework to populate forwarding whitelists at each data concentrator so that (i) load is balanced among the different concentrators and (ii) the overall performance of the network in terms of packet reception rate and received signal strength are maximized. We also propose a heuristic algorithm to generate near optimal forwarding whitelists in acceptable computing time. Extensive experiments are performed on a real scenario consisting of a city wide gas meter network deployed in northern Italy. Results show that the proposed heuristic is able to produce whitelists that reduce the average backhauling traffic as much as 80%, with a corresponding network quality within 4% of the one computed by the optimal solution.

I. INTRODUCTION

The last few years have witnessed a rapid growth of systems and technologies related to the Internet of Things, in particular for what concerns applications related to Smart Cities. In this context, the field of Smart Metering is emerging as one of the most concrete examples of how connected things may bring improvements to the everyday life of both citizens and municipality administrators. As a matter of fact, smart meters roll-outs plans are already active in several countries. To mention a few significative examples, Germany's electricity grid operator will install about 7.5 million smart energy meters across Germany between 2018 and 2020 as part of the transition from centralized energy generation to distributed energy resources [START_REF]Gesetz zur digitalisierung der energiewende[END_REF]. Italy is also following this trend: the Italian Authority for Energy, Gas and Water has recently opted for a mandatory roll-out according to which services operators will provide nearly 40 million energy customers and 22 million gas and water customers with a smart meter network over the next few years.

Smart meters have unique requirements in terms of bandwidth, energy consumption and costs. On the one hand, they do not require high data rate as the number of measurements has very coarse time granularity (e.g. 2/4 measurements per day). On the other hand, smart meters are generally battery operated and required by national regulations to have a minimum of 10-15 years of autonomy. At the same time, their cost should be limited in order to make it convenient for services operators to deploy them in place of the traditional meters. These requirements clearly have an impact on the choice of the communication technology to be used for smart meters. Solutions based on GSM/GPRS, where the meters are directly able to communicate their measurements to a central server have the drawback of being costly and power eager. A solution which respects the aforementioned requirements and has become the de-facto standard for many smart meter deployments in Europe is the wireless M-Bus (wM-Bus) protocol.

In the reference architecture of a city-wide wM-Bus network, depicted in Figure 1, hundreds of thousands of meters communicate their measurements to a small number (tens or hundreds) of data concentrators (DC) which are deployed in order to provide network coverage. Due to the harsh propagation conditions of urban environments, data concentrators are generally deployed in elevated positions (e.g. light poles or other structures) and use GSM/GPRS as backhauling technology for relaying the meter data to a central server, as other wireless or wired backhaul technologies with higher throughputs are either too expensive or impractical to set up. Due to the association-less nature of the wM-Bus protocol, data transmitted by one meter can be received at multiple data concentrators. This redundancy has a twofold harmful effect on the entire system: (i) on the uplink, it reduces the capacity of the already limited GSM/GPRS backhauling links and it causes computational overload at the DCs and at the central server; (ii) on the downlink, it creates an ambiguity about which DC should deliver system commands to each meter (e.g. increase/decrease measurement rate, close the main gas/water valve for insolvent customers).

In principle, a careful network planning could avoid the redundancy caused by duplicate packets: this could be obtained by minimizing the overlap between the coverage areas of different DCs or with frequency reuse techniques. However such planning strategies are seldom used, because both DCs and meters are deployed incrementally where needed rather than in a global single deployment process. For this reason, we focus in this paper on the creation of whitelists to be used Fig. 1: Reference architecture: smart meters of different utilities (gas/water) communicate their measurements to multiple data concentrators using the wM-Bus protocol. Data concentrators relay such information to a central data management service through GPRS communication for automatic billings and management procedures on each DC for deciding which meter should be managed (by forwarding its data or sending commands to it) or not. We propose an optimization framework to create whitelists with the following features:

1) they balance as much as possible the amount of meters that each DC manages; 2) they provide the best associations between meters and DCs, in terms of link quality; 3) they can be easily updated in order to cope with the installation of new meters or with changes in the RF propagation environment; The problem is tackled using a data-driven approach to obtain information on the current propagation conditions between meters and DCs. Such data is then fed into an optimization problem that produces balanced whitelists, also maximizing the link quality between meters and DCs. A heuristic is also proposed to provide quick updates even on big instances of the problem, such as the ones arising from real deployments in urban scenarios.

The remainder of this paper is organized as follows: Section II briefly summarize the wM-Bus protocol; Section III introduces the optimization framework, gives the ILP formulation and focuses on the proposed heuristic algorithm. Results obtained on a real smart meter deployment are shown in Section IV. Finally, Section V reports on the related literature and Section VI concludes the paper.

II. WIRELESS M-BUS PROTOCOL

The wM-Bus protocol was introduced in 2005 as a European standard (EN 13757-4:2005), specifying the communication primitives between utility meters and data loggers at the physical and data link layer. The presence of a European standard to regulate how such communications should be performed allows to create an ecosystem where utility providers, meter manufacturers and data management companies can coexist to create better services for the end users as well as business opportunities. The current version of the standard (EN 13757-4:2013) specifies several modes of operation for smart meters. Each mode is characterized by a specific carrier frequency (around 434MHz, 868MHz or 169MHz), data rate, duty cycle and whether it supports unidirectional or bidirectional communication 1 . Table I summarizes the available communication modes. The choice of the particular mode to use is left to each country regulation: Italy and France, two of the main smart meters markets in Europe, use narrowband mode N at 169 MHz. The advantages of using a narrowband communication consist mainly in the limited noise resulting at the receiver, which greatly increases the available link budget and allows to decode frames even at a sensitivity as low as -115 dBm. The reduction in data rate compared to modes operating at higher carrier frequencies is not critical since the amount of data to be transmitted by each meter is very low. At the same time, using a lower carrier frequency allows for a greater transmission range, which can reach 1Km even in urban scenarios [START_REF] Spinsante | Evaluation of the wireless m-bus standard for future smart water grids[END_REF]. At the data link layer, the standard specifies two different frame formats (A and B): for uplink communications from meters to DCs, both formats provide an 8 byte sender address field, but due to the association-less nature of the protocol, no destination address field is contemplated 2 . A frame transmitted by a meter can be therefore accepted by more than one DC, and relayed on the backhauling infrastructure more than once. Therefore, without proper filtering mechanisms implemented at DCs, the efficiency of the entire network can be greatly affected. In the following sections we target precisely this issue.

III. OPTIMIZATION FRAMEWORK

A. Integer programming formulation

Keeping in mind the reference architecture of Figure 1, we consider a network of C data concentrators and M meters uniformly deployed over a urban area. Let C be the set of the concentrators while M the set of meters. We define the Range Set R c for the c-th concentrator as the set of in-range meters, containing all those meters from which concentrator c has received at least one frame in a specified time window. Similarly, let R m for the m-th meter be the set of in-range concentrators. Let x m,c be a binary association variable which is equal to 1 if the m-th meter is associated with the cth concentrator and 0 otherwise. The whitelist of the c-th concentrator therefore contains all those meters m for which x m,c = 1. Note that, by definition, x m,c = 0 if m / ∈ R c , that is a meter cannot be in the whitelist of a concentrator if it is not in its Range Set. Finally, let q m,c a measure of the link quality existing between meter m and concentrator c. Note that q m,c = 0 if m / ∈ R c . Under these definitions, we propose a nonlinear integer programming (NIP) formulation for obtaining balanced whitelists while maximizing the overall quality of the wireless links:

max xm,c J = α • c∈C log(L c ) + β • m∈M c∈C x m,c • q m,c (1) s.t. ∀c∈C x m,c = 1 ∀m ∈ M (2)
L c = ∀m∈M x m,c ∀c ∈ C (3)
x m,c = 0 ∀c ∈ C, m / ∈ R c (4)
The objective function ( 1) is a weighted sum of two contributions. The first term aims at balancing the amount of meters L c in each concentrator's whitelist. Balancedness is attained through the use of the logarithmic function, which is known to provide fairness and thus it is a common choice in load balancing problems [START_REF] Bu | Generalized proportional fair scheduling in third generation wireless data networks[END_REF], [START_REF] Ye | User association for load balancing in heterogeneous cellular networks[END_REF]. The second term of the objective function maximizes the total quality of the wireless links resulting after the whitelist creations. Moreover, the objective function can be shaped according to a user's preference through setting proper values for parameters α and β. Constraint (2) ensures that one meter is associated to the whitelist of only one concentrator, constraint (3) defines the load of each concentrator as the sum of associated meters and finally constraint (4) enforces that a meter can be inserted in the whitelist of a concentrator only if it is in that concentrator Range Set.

B. Heuristic algorithm

The optimization problem ( 1)-( 4) can be considered a Generalized Proportional Fairness (GPF) problem, whose solution is known to be NP-Hard [START_REF] Bu | Generalized proportional fair scheduling in third generation wireless data networks[END_REF]- [START_REF] Rakotomanana | Fair load balancing in heterogeneous cellular networks[END_REF]. At the same time, meters and gateways are often deployed incrementally rather than in a single solution, thus requiring to often recompute the aforementioned whitelists. Therefore, it is important to provide a way for quickly obtaining a good solution to the optimization problem, even though non optimal. Here we introduce a heuristic algorithm able to do so. The basic idea of the heuristic is to separately optimize the two terms of the objective function, taking first the perspective of the DCs and secondly that of the meters, alternating between the two until a stop condition is reached.

The heuristic starts by creating an initial solution S 0 , obtained by sorting the DCs in ascending order of the cardinality of the corresponding R c and then assigning meters to DCs in a round robin fashion. At each turn of the round robin schedule, the c-th concentrator gets assigned the meter among the nonassigned ones for which q m,c is maximized. Note that S 0 is balanced by construction, as each concentrator is assigned the same number of meters (but one). The heuristic then runs for max iter iterations, each one consisting of two phases. In the first phase, meters are sorted in descending order of the associated Maximum Metric Gain G m , defined as it follows:

G m = | c∈C x m,c • q m,c -q m,cm | (5)
where

(6) c m = argmax c∈Rm q m,c (7)
As one can see, G m captures the maximum improvement of the link quality of a meter if it is removed from its current whitelist and inserted in the whitelist of the best concentrator quality-wise, c m . The algorithm then proceeds in performing such exchanges one at a time. However, while each exchange improves the second term of the objective function [START_REF]Gesetz zur digitalisierung der energiewende[END_REF], it does so without taking balancedness (i.e., the first term of eq. ( 1)) into account. Therefore, to avoid potential heavy load unbalances among DCs, after each exchange the algorithm evaluates the term c∈C log(L c ) and compares it with its initial value L 0 , according to S 0 . Meter exchanges are interrupted if the decrease in the sum of the logarithms is greater than a predefined threshold ξ max , that is when:

L 0 -c∈C log(L c ) L 0 > ξ max . (8)
When this happens, all meters that have been moved are inserted into a Blocking List: such meters will be prevented to be exchanged again in future iterations. In the second phase, the heuristic re-balances the loads through redistributing the non-blocked meters to concentrators with the same rationale used in producing the initial solution S 0 . In details, concentrators are sorted in ascending order of |R c | (this time without considering the meters in the Blocking List) and non-blocked meters are assigned one at a time in a round robin fashion, each time to the concentrator maximizing the correspondent q m,c . We proceed with this phase until all the non-blocked meters have been redistributed, producing a better balancedness of the loads. This concludes the i-th iteration of the algorithm. At each iteration, a partial solution S i and its corresponding objective function value J(S i ) is evaluated.

When max iter iterations are executed, the best solution S is returned. A formal description of the whole algorithm is given in Algorithm 1.

Algorithm 1 Heuristic algorithm

Input: max iter; ξ; q m,c , R c , R m ∀m ∈ M, ∀c ∈ C Output: Association Indicators x m,c (whitelists) {Step 0} x m,c = 0, ∀m ∈ M, ∀c ∈ C Sort the DC in ascending order on R c cardinality Set i = 1 while (∃m | c∈Rm x m,c = 0 ) do Pick from the list concentrator c i m = argmax m q m,ci x (m ,ci) = 1 i = i + 1 mod C end while Set S 0 = {x m,c ∀m ∈ M, ∀c ∈ C} Set S = S 0 Compute L 0 = c∈C log(L c ) Set Blocking List = ∅ for (i = 1 : 1 : max iter) do {Phase 1} Set k = 1 Compute G m ∀m ∈ M Sort meters in descending order of G m while ( L0-c∈C log(Lc) L0 < ξ max ) and k < M do Pick from the list meter m k if m k / ∈ Blocking List then x (m k ,c) = 0, x (m k ,cm k ) = 1 Insert m k in the Blocking List end if k = k + 1 end while {Phase 2}
x (m,c) = 0 ∀m / ∈ Blocking List, ∀c ∈ C Sort the DC in ascending order on R c cardinality, ignoring the meters in the Blocking List

Set k = 1 while (∃m | c∈Rm x m,c = 0) do Pick from the list concentrator c k m = argmax m q m,c k x (m ,c k ) = 1 k = k + 1 mod C end while Set S i = {x m,c ∀m ∈ M, ∀c ∈ C} if J(S i ) > J(S) then S = {S i } end if end for

C. Data collection

The NIP formulation and the proposed heuristic algorithm require the knowledge of the link quality q m,c between each meter and a DC. We observe that, beside the application data (i.e., the meter measurements themselves), data frames received at the different DCs provide a priceless source of information for estimating the link quality q m,c . We propose to perform local processing of the received frames at each DC with the objective of retaining important information from meters while at the same time reducing the amount of data to be transmitted at the data management server. In addition to the aforementioned whitelist, each DC maintains a local database for storing meters related information. For each meter m, a concentrator stores the number of received messages from the meter n m,c and the average Received Signal Strength (RSS) µ m,c . Upon reception of a frame from a meter m with RSS s, a data concentrator c performs the following operations:

1) Check if m is already present in the local database. If not, create a new entry for that meter and set n m,c = 0, µ m,c = 0. 2) Update the statistics relative to m in online fashion:

n m,c = n m,c + 1 (9) µ m,c = (n m,c -1) • µ m,c + s n m,c (10)
3) Check if m is present in the current whitelist. In case it is, forward the frame to the data management server. Otherwise, discard the frame. At the end of a predefined time window (e.g., 2 weeks), the data concentrator transmits a copy of the local database to the data management server and clears all information. Such an approach allows to limit both the memory requirements on the data concentrator and the amount of resources needed for transmitting such information. Note that since the transmission rate of meters is known a priori (e.g., 2/4 frames per day), observing the number of received messages n m,c over a predefined period of time can be considered as a good estimator of the packet error rate (PER) p m,c . Being f the transmission rate of a meter (in frames per day) and w the length in days of the time window, the PER can be estimated as:

p m,c = 1 - n m,c f • w (11)
Normally, the PER is limited between 0 and 1. However, if a node is subject to clock-related issues (e.g., transmits more messages than what set per day), the PER value may also be lower than 0. In this case, simply checking the PER value may identify malfunctioning meters that may deplete their batteries sooner than expected. Finally, the indicators p m,c and µ m,c can be conveniently fused to estimate the link quality q m,c :

q m,c = u • µ m,c -v • p m,c (12)
where u and v are tunable weighting coefficients which can be adjusted according to a user's preference. Also, considering that variation ranges of µ m,c and p m,c can be very different (µ m,c typically varies between -10 dbm and -120 dbm while p m,c is normally limited between 0 and 1), for numerical stability they are normalized to unity, such that both vary within 0 and 1.

A. General Description and Experiments Setup

We consider three instances of a real smart meter network deployed over a urban area in the north of Italy. The network is composed of about 10000 smart gas meters and 201 concentrators, operating the wM-Bus protocol using mode N. The backhaul link from each concentrator to the network server is based on GPRS technology. DCs collect channel information as detailed in Section III-C and transmit it every 14 days using the MQTT protocol to a central server. There, a centralized MySQL database is used to store all received information. From this scenario, we created three network instances by sampling spatially three windows of different sizes, containing M = [373, 1931, 9512] meters, always considering all available concentrators. Windows sizes have been empirically set with the rationale of considering problems with an increasing number of variables such to compare the quality of the obtained solutions and the corresponding execution times. We refer to these instances as Small (S), Medium (M) and Large (L) in Table II. The experiments are characterized by four degrees of freedom which regard on the one hand the channel quality metric definition (u, v in ( 12)) and on the other hand the weights of the utility objective function (α, β in ( 1)). Believing that neither optimal channel conditions with unbalanced loads among the DCs nor perfect load balancedness with poor channel conditions can be considered desired network states for the application under study, we set α and β such that the tradeoff between the two system aspects we consider is balanced. In practice, α is set to c∈C |R c |/10 4 while β is set to 1. For what concerns the channel quality part of the objective function, we investigate three cases, considering i) only µ m,c ((u, v) = (1, 0)), ii) only p m,c ((u, v) = (0, 1)) and iii) a convex combination of the two ((u, v) = (0.3, 0.7)). These three cases are referred to in Table II as: Load+RSS (S1,M1,L1), Load+PER (S2,M2,L2) and Load+RSS&PER (S3,M3,L3). A further investigated case is referred to in Table II as Blind Optimum, where we suppose that channel quality information is not available and the channel quality metric q m,c is replaced in the optimization problem ( 1)-( 4) by -d m,c , where d m,c is the physical distance between the m-th meter and the c-th concentrator. This case reflects the scenario of the very first network deployment, where whitelists are computed assuming that the shorter the meter-DCs distance, the higher the channel quality. For what concerns other settings, the NIP is formulated with AMPL and solved with the BONMIN 1.3 solver on a Intel Core 2 Duo processor with P9700 chipset, 2.8 GHz CPU and 4 GB RAM running on the platform Microsoft Windows 7. We choose to limit the maximum running time of the NIP to 8 hours; after this, the solver returns the best solution found so far, if the problem is feasible. The heuristic algorithm is implemented in Matlab, with max iter set to 1000 and ξ max to 5%.

B. NIP Output

We first analyse the results obtained by solving the NIP formulation to optimality. Figure 2a and 2b show the distributions of the RSS and PER resulting after the computation of the optimal whitelists, considering different quality metric parameters. When u = 1 and v = 0 the channel metric optimizes the RSS, while for v = 1 and u = 0 the problem focuses on the PER. As one can see, in the former case (red curves) the average RSS is 3dBm higher than in the latter case, whereas the average PER is almost 10% better when explicitly targeted (blue curves). Looking at the results reported in Table II for the small instance, we observe that targeting only the PER is computationally cheaper than targeting only the RSS (or both) by almost four times. For larger instances however, due to the combinatorial structure of the ILP, execution times are greater than 8 hours. This cannot be acceptable considering that the network is expanded incrementally according to specific needs and thus requiring quick and easy updates of the association patterns. Table II shows that both the M and the L instances have led the NIP to exit for time limit, investigating just a small subgroup of the solution space; only the S instance returns the optimum association in a reasonable, still too large, time frame. Moreover, as shown in Table II, the results obtained from the Blind case are not close to the optimal ones: the gaps are greater than 8%, approaching in some cases 14%. This is expected, considering the harsh propagation conditions usual of urban environments. In such a scenario, the meter-DC couple with the shorter distance is not necessarily the best association.

C. Heuristic vs Optimal Solution

The heuristic algorithm is meant to give acceptable suboptimal solutions in shorter execution times compared to those of the NIP. The algorithm tries to improve as much as possible the objective function value from the starting point of the S 0 solution. Figure 3 plots the evolution of the objective function value during the execution of the heuristic for one of the considered case study. As reported in Table II, the solutions obtained by the proposed heuristic are very close to the optimal ones, even after relatively shorter execution times: gaps are limited to less than 4% of the optimal solution. It is worth looking at the first term of the objective function and at the final quality metrics resulting from the computed whitelists, which contribute to define the second term of the objective function. Figure 4 shows that the load term of J in the heuristic case is close to the optimum one, being even greater in some cases: what happens is that when the heuristic load term dominates the optimal one then the opposite happens for the quality metric term and vice versa, the sum of the two being always greater in the optimum case. Figure 5a illustrates the average RSS resulting from the Optimum and Heuristic in the three scenarios, while Figure 5b does the same for what concerns the PER. As one can see, when the focus is on RSS (red bars), the average RSS is between 1 and 2 db higher

D. Reduction of the Backhauling Load

One of the main motivation for the adoption of whitelists, disregarding the strategy implemented to create them, is to limit the amount of traffic to be forwarded from data concentrators on their GPRS-based backhauling links to the central server, as well as the computational load the server is subject to. Figure 6 reports the average and standard deviation of the traffic reduction that can be obtained on each gateway after the activation of whitelists in the three considered network scenarios and according to the three different combinations of u and v parameters, with respect to a whitelists-less situation. As one can see, the average per-concentrator reduction is between 65% and 80% regardless of the considered scenario, with peaks of more than 100% (i.e., the backhauling traffic is halved). At the central server, such reduction is even more important, as it is multiplied by the number of concentrators.

V. RELATED WORK

Several works in the last few years have given attention to the wireless M-Bus protocol, either from a performance evaluation point of view [START_REF] Masek | Wireless m-bus: an attractive m2m technology for 5g-grade home automation[END_REF] or from a system design perspective [START_REF] Mihajlović | Implementation of wireless m-bus concentrator/gateway for remote reading of smart gas meters[END_REF]. In [START_REF] Spinsante | Evaluation of the wireless m-bus standard for future smart water grids[END_REF], [START_REF] Squartini | Wireless m-bus sensor nodes in smart water grids: the energy issue[END_REF] and later in [START_REF] Spinsante | Wireless m-bus sensor networks for smart water grids: analysis and results[END_REF], the application of wM-Bus protocol to smart water grids is analyzed. Authors focus in particular on energy requirements of wM-Bus compliant smart meters, and reason on the possibility of using water kinetic energy as an energy harvesting source, also taking into consideration the transmission range and the collision probability of each meter. Using both simulations and experimental tests with a system prototype, authors conclude that the wM-Bus protocol is an adequate tradeoff for the applicative context under study. In [START_REF] Carratù | Performance analysis of wm-bus networks for smart metering[END_REF], authors analyze statistical properties of frames transmitted by wM-Bus gas/water meters and received at different DCs during one year. The analysis reveals that regulating the transmission power solely according to the distance from the nearest DC is not enough to provide adequate transmission quality, and suggests to use the monthly average value of the RSS as a predictor to tune the transmission power for the next observation period. In this work, on the one hand we confirm that distance-based approaches do not capture channel quality in all its aspects and on the other hand we investigate a new channel quality metric, evaluating its efficiency to drive the network towards better system performances.

The problem of balancing the number of users connected to each gateway/base station is of great importance in the context of both mobile cellular and wireless sensor networks. In cellular networks, the problem of unbalanced loads among cells deteriorating system-wide performances is of great concern and is usually considered in the context of resources allocation to new users approaching the network and looking for a base station to associate. In [START_REF] Bu | Generalized proportional fair scheduling in third generation wireless data networks[END_REF] and later in [START_REF] Ye | User association for load balancing in heterogeneous cellular networks[END_REF], authors propose an ILP formulation which aims at finding the best associations between users and base stations in order to achieve load balancing among base stations and some level of fairness among the bandwidth allocated to users. Authors underline that while linear utility functions results in throughput-optimal non-balanced solutions, concave objective functions, such as the logarithm, naturally achieve fairness among users. An interesting option is given also in [START_REF] Wang | Dynamic load balancing and throughput optimization in 3gpp lte networks[END_REF], where a weigthed utility objective function is proposed, such that weigthing coefficients on a load balance index and network throughput are properly tuned according to the desired level of tradeoff between the two. In sensor networks, where equipments are usually static and requirements in terms of throughput are much coarser, the problem of balancing the load among a-priori chosen DCs is tackled by properly clustering the sensor nodes such that each sensor belongs to one and only one cluster and communicates with the central server only via the head of that cluster, as presented in [START_REF] Gupta | Load-balanced clustering of wireless sensor networks[END_REF]- [START_REF]Energy efficient load-balanced clustering algorithm for wireless sensor networks[END_REF]. Load balancing among cluster heads is achieved by either minimizing the variance of the number of sensors connected to each cluster head [START_REF] Gupta | Load-balanced clustering of wireless sensor networks[END_REF] or by taking a min-max approach [START_REF] Low | Load-balanced clustering algorithms for wireless sensor networks[END_REF]. In [START_REF]Energy efficient load-balanced clustering algorithm for wireless sensor networks[END_REF] a distance based, load aware clustering algorithm is presented where gateways are sorted at each iteration in ascending order on the number of so far assigned sensors (min-heap tree), which are then associated to the nearest gateways following the aforementioned list. A similar algorithm but working from sensors' side is shown in [START_REF] Kuila | Improved load balanced clustering algorithm for wireless sensor networks[END_REF], where sensors are sorted in ascending order on the number of in-range gateways and then assigned to the best gateway according to the chosen metric. We note that there is no evidence in literature about a joint use of the above perspectives and this work shows a possible approach in this context. Moreover, beside the focus on load balancing, our work focuses also on system performances, introducing new aspects different than the distance in order to achieve better results.

VI. CONCLUSION

In this paper we propose a data-driven optimization framework to jointly balance the loads and maximize the channel quality of a smart meter network. This is crucial for the reduction of the traffic on the limited-capacity backhauling links, the reduction of the computational load at the network gateways towards the central server and the establishment of a safer and less-ambiguous association-based communication paradigm between meters and gateways. We show that the information contained in the data frames sent by the sensors can be stored at DC side and exploited to lead the system to better performances in terms of Packet Error Rate, Received Signal Strength and load balancedness. Then, we propose a heuristic approach which is able to give good solutions in proportionately smaller execution times, which is crucial relatively to the need of smart metering networks. In fact, considering that usually such networks face more than one deployment phases, they strongly need to quickly update their network association states if the case requires it and so computational times cannot be too long. The gaps between the Optimum and the Heuristic objective function values show that heuristic solutions are close enough to optimal ones, being therefore a good approximation. Future works could regard, among the others, the definitions of finer channel quality metrics, the formulation of quicker and more-performing heuristics and the consideration of different utility functions.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Sensibility of the ILP to weights u and v

Fig. 4 :

 4 Fig. 4: First term of J in Optimum and Heuristic cases, light and dark bars respectively

TABLE I :

 I wM-Bus communication modes

	Mode	Frequency	Max. Duty Cycle	Description
	S (Stationary)	868 MHz	0.02 -1%	Meters transmit few samples per day
	T (Frequent Tx)	868 MHz	0.1 -1%	Meters transmit short frames every few seconds
	R (Frequent Rx)	868 MHz	1%	Meters listen every few seconds for reception of a wake up message and starts transmission. Used for readout of several meters at once.
	C (Compact)	868 MHz	0.1%	Optimized version of mode T
	N (Narrowband)	169 MHz	10%	Long range, narrowband system
	F (Frequent Tx/Rx)	433 MHz	10%	Frequent bidirectional communication

TABLE II :

 II Numerical Results. S,M and L instances have respectively: 1436, 10828 and 82616 number of variables. OF stands for Objective Function, ET stands for Execution Time

			Optimum	Blind Optimum	Heuristic	Heur vs Opt
			OF	ET (min)	OF	Gap from Opt (%)	OF	ET (min)	Gap (%)
		1	229,94	256,2	211,37	8,1	223,91	0,18	2,6
	S	2 3	414,25 387,19	66 285,6	358,40 342,53	13,5 11,5	400,49 373,43	0,17 0,19	3,3 3,6
		1	1448,63	>480	1319,71	8,9	1400,46	1,46	3,3
	M	2 3	1855,52 1763,15	>480 >480	1593,26 1550,54	14,1 12,1	1794,36 1700,76	1,45 1,46	3,3 3,5
		1	9446,19	>480	8437,46	10,7	9173,90	36	2,9
	L	2 3	12006,6 11173.9	>480 >480	10470,82 9860,81	12,8 11,8	10770,54	33 35	3,1 3,6
			IV. EXPERIMENTS

Each mode is named with a letter and a number. The letter specifies the mode and the number specifies unidirectional (=1) or bidirectional (=2) communication

The EN 13757-4 standard actually specifies an optional manual association between a meter and a DC through a three-way handshake process. However, no specifications are provided on how the DC should be selected or what happens in case multiple DCs receive the association request.

Comparison of Opt/Heur performances on RSS

Large

Medium Small

Backhauling Traffic Reduction after Heuristic

Fig. 6: Mean per-GW Backhauling Traffic Reduction in Optimum than in Heuristic in M and L instances while the Heuristic gives slightly higher RSS in S instance. When PER is targeted (blue bars), in L instance the Heuristic gives almost 15% higher PER than the Optimum while the gap in M instance is reduced to 6%. The situation is reversed in S instance, where the Heuristic gives 9% better PER with respect to the Again, where Heuristic performance are better from a quality metric perspective, the Optimum does better on the load side. We anyway underline that the tradeoff between the two terms of the objective function depends on the choice regarding the setting of the parameters α and β, which can be tuned according to the user's preference.