Carlos Canudas De Wit
email: carlos.canudas-de-wit@gipsa-lab.fr

Federica Garin
email: federica.garin@inria.fr

Ruggero Fabbiano
email: ruggero.fabbiano@inria.fr

Pierre Rouchon
email: pierre.rouchon@mines-paristech.fr

Antoine Rousseau
email: antoine.rousseau@inria.fr

Antoine Rousseau Source

Poisson Integrals Necsys

Source localization using Poisson integrals

Keywords: Source seeking, sensor networks, multi-agent systems

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

This paper deals with the problem of source localization via several sensor devices providing pointwise concentration measures; sensors are assumed to be arranged in circular arrays, they can be fixed along the array or they can turn along a circular path defined by the array.

Source localization is relevant to many applications of vapor emitting sources [START_REF] Porat | Localizing vaporemitting sources by moving sensors[END_REF] such as explosive detection, drug detection, sensing leakage or hazardous chemicals, pollution sensing and environmental studies; sound source localization [START_REF] Zhang | Maximum likelihood sound source localization and beamforming for directional microphone arrays in distributed meetings[END_REF] is pertinent for intelligent conference calls systems that identify the speakers to improve sound and video quality; other applications also include heat source localization and vent sources in underwater field.

There exist in the literature a variety of methods to treat the problem of source localization and related issues. Many techniques deal with formulations associated with isotropic diffusion processes described by diffusion equations for which a closed-form solution is known; as the explicit solution depends on the source location (among other parameters), several identification methods have been de-⋆ The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007(FP7/ -2013) ) under grant agreements n. 223866 FeedNetBack and n. 257462 HYCON2 network of excellence.

vised to estimate the source position: [START_REF] Matthes | Source localization based on pointwise concentration measurements[END_REF] proposed a two-step identification procedure, dealing with the inhomogeneous case and a fixed sensor array; [START_REF] Porat | Localizing vaporemitting sources by moving sensors[END_REF] formulated a similar problem but with moving sensors, using a maximum-likelihood approach to estimate the source position, and considering moving sensors which update their position so as to approximately minimize the estimation error, by following the gradient of the Cramér-Rao bound to error variance. More fundamental problems such as source identifiability and optimal sensor placement are discussed in depth by [START_REF] Khapalov | Source localization and sensor placement in environmental monitoring[END_REF] using concepts and ideas of control system theory. The above mentioned approaches can be viewed as inverse problems for partial differential equations, with the goal of finding the initial conditions or a forcing term. Because of their nature, all such methods share the common drawbacks of heavy computations, and of high sensitivity to the explicit knowledge of the closed-form solution of the pde describing the diffusion process.

Another related but quite different approach for source localization is based on extremum seeking techniques [START_REF] Ariyur | Real-Time Optimization by Extremun-Seeking Control[END_REF]. In contrast to the methods mentioned previously, this approach is not based on any particular structure or knowledge of the diffusion solution; the method only applies for moving sensors, as it relies on the idea of collecting rich enough information to approximate the gradient through the use of a periodic probing signal. Adaptations of this idea to the problem of source localization with a non-holonomic unicycle have been previously reported in [Cochran and Krstić, 2007] for the stability study, and by Cochran et al. (2007) for application of this idea to several scenarios.

The method presented here follows a distinct direction than the ones previously cited: on the one hand, it does not require specific knowledge of the solution of the diffusion process, but only exploits properties of the pde, and it can compute the gradient direction from the pointwise concentration samples with a small computation load; on the other hand, it does not make use of a probing signal and thus avoids the oscillations required by extremum seeking techniques. The main idea consists in using Poisson integrals to compute the gradient necessary to perform source search, and also other higher-order derivatives (in particular, the Hessian), which can be useful to implement different control laws; the main assumption is that the diffusion is described by the Laplacian pde, i.e., that the Laplacian of the measured value is null on the studied domain: such an assumption is justified for isotropic diffusive sources in steady-state. Other potential benefits are: extensions to 3d source localization (will be reported elsewhere), and intrinsic high-frequency filtering (derivatives are computed using integrals) that makes the methods low sensitive to measurement noise.

The paper formalizes and extends previous ideas from [START_REF] Moore | Source seeking via collaborative measurements by a circular formation of agents[END_REF][START_REF] Moore | Source seeking via collaborative measurements by a circular formation of agents[END_REF][START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF], where the gradient has been approximated by the sum of pointwise measurements around a circle weighted by the position vector of each sampler to its center of rotation.

PROBLEM FORMULATION

We consider diffusion processes where the source is isotropic and the diffusion is homogeneous, so that the diffusion is described by the well-known diffusion equation (1). For the sake of simplicity, we focus on the 2dimensional case, where the diffusion happens in the (x, y)plane; we will use the notation z = (x, y) to denote a point z ∈ R 2 , and, whenever convenient, we will also use the notation z = ρe iθ to denote the point (x, y) = (ρ cos θ , ρ sin θ), using the usual bijection between points in R 2 and elements of C. However, all operations will be intended as operations on reals and on vectors in R 2 ; in particular, all derivatives are intended as (partial) derivatives of real-valued functions, and integrals are intended as entry-wise integrals of the real-valued entries of the vectors. For an open set Ω ⊂ R 2 , we will denote its border with ∂Ω and its closure with Ω = Ω ∪ ∂Ω.

The 2-dimensional isotropic and homogeneous diffusion is described by the following linear parabolic pde with constant coefficients, known as "isotropic diffusion equation" or as "heat equation"

∂f (z, t) ∂t -κ∆f (z, t) = 0 , ∀z ∈ Ω , t ≥ 0 (1)
where Ω is an open subset of R 2 , ∆ is the Laplacian operator defined by ∆f (z) = ∂ 2 f ∂x 2 + ∂ 2 f ∂y 2 , κ is the isotropic diffusion coefficient. Such an equation can describe various diffusion phenomena: for instance, f can represent the temperature or the concentration of a chemical (e.g., a pollutant, or salinity).

In this work, we assume that the diffusion process is fast enough, so that our main interest is in studying only the steady-state behavior resulting from equation (1); therefore, we limit our attention to solutions of the following equation, known as Laplace equation ∆f (z) = 0 , ∀z ∈ Ω .

Our interest is in the case where the steady-state has been reached but a source is still emitting somewhere at a constant rate, and our goal is to find the source location. As a model for such a source, we assume that the source occupies a portion of space not belonging to Ω, and it affects the values of f in Ω by imposing a boundary condition. More precisely, we consider an open domain Ω = Ω \Ω s , where Ω ⊂ R 2 is a connected bounded set representing the region we are interested in studying, and Ω s is a small connected subset of Ω which represents the area occupied by the source. Hence, the boundary of Ω is formed by two parts: an external one, equal to the boundary of Ω and denoted by ∂Ω ext , and an inner one, equal to the boundary of the source region Ω s and denoted by ∂Ω in .

As an illustrative family of examples, we consider a source of heat in a room, where the source is modeled as a small circular region Ω s imposing a constant given value of f at the inner boundary ∂Ω in ; within this family, we consider at first a simple example, where f can be written in closed form, and then a richer example, where f can be computed numerically with the finite elements method.

Ω s = B ρ (c),
where the center c is the same for both balls, but the radius ρ is much smaller than r. In this example, we assume that the boundary conditions are a constant temperature T s around the source and a smaller constant temperature T ext on the wall, so that we are looking for a solution of the following Dirichlet problem:

   ∆f (z) = 0 in Ω = B r (c) \ B ρ (c) f (z) = T s on ∂Ω in = ∂B ρ (c) f (z) = T ext on ∂Ω ext = ∂B r (c) .
(2)

It is easy to verify that all functions of the form f (z) = α log z -c + β satisfy ∆f (z) = 0 for all z = c, so in particular for all z ∈ Ω; moreover, it is easy to see that such functions are constant along circles centered in c. Then, by imposing the given boundary conditions, one can find the correct values of α and β, and can find the following solution f to the problem (2):

f (z) = T ext log( z -c /ρ) -T s log( z -c /r) log(r/ρ) .
Thanks to the properties of the Laplacian operator, and to choice of the boundary conditions, this is actually the only solution of the problem (2). a rectangular room: similarly to Example 1, a circular heater Ω s lies in the middle of the room, and imposes a constant boundary condition f (z) = T s for all z ∈ ∂Ω in .

On the walls, the boundary condition is imposed not on f , but on the derivative of f in the direction of the outward unit normal n: we assume that ∇f (z) • n = 0 along the walls, which models isolating walls; in the middle of the back wall, a large window is open, and imposes a constant boundary condition f (z) = T ext due to the outer colder temperature. The temperature profile at the steady state can be obtained by solving the Laplacian equation with the above-described boundary conditions by the use of a finite-element method. Figure 1 shows the temperature obtained by solving the above-described problem via the FreeFem++ software (see [START_REF] Hecht | FreeFem++ manual[END_REF]), for a room 10 m wide and 6 m long; the window is 6 m large and the heater occupies a circle with a radius of 50 cm; the heater's temperature is T s = 45 ℃, and the temperature outside the window is T ext = 10 ℃.

The solutions of the Laplace equation ∆f = 0 are called harmonic, and have particular properties; in the next section we will recall some of these properties, which we will then exploit in order to design algorithms for source seeking.

Because in our setup the inner boundary ∂Ω in represents a source, we assume that values of f on ∂Ω in are higher than values of f on ∂Ω ext ; under this assumption, the maximum principle (see Proposition 1) ensures that the maximum value of f on Ω is attained on ∂Ω in . Hence, the problem of finding the source can be described as the problem of finding the maximum value of f on Ω; more precisely, having assumed a constant value for f along ∂Ω in , we will consider that the source seeking problem is solved if any point along ∂Ω in has been reached.

GRADIENT COMPUTATION USING POISSON INTEGRAL FORMULA

In this section we present the preliminaries, and the method based on Poisson integrals that we will use to compute the function f (z) and its n th order derivatives at any point inside a circular region by using only measurements along a circular path (the border of the circle where the sensors are placed); in particular, we will provide specific formulas to estimate the gradient ∇f (z 0 ) at the center of the sensor array by using only informations from the circular path.

Harmonic functions

We start by recalling some important properties of harmonic functions, which can be found for instance in textbooks such as [START_REF] Folland | Introduction to Partial Differential Equations[END_REF] and [START_REF] Axler | Harmonic Function Theory[END_REF].

Definition 1. (Harmonic function) Let Ω ⊆ R 2 be an open set. A function f : Ω → R is harmonic in Ω if f ∈ C 2 (Ω) and ∆f (z) = 0 for all z ∈ Ω.
The above definition of a harmonic function requires the function to be twice-differentiable. Then, it is possible to prove that harmonic functions have a much higher regularity: if f is harmonic on Ω, then f is analytic (Folland, 1995, Coroll. 2.11), and in particular f ∈ C ∞ (Ω). Moreover, the assumption that f ∈ C 2 (Ω) in the definition of a harmonic function is convenient for having welldefined second-order partial derivatives, but it is not an essential assumption: any distribution f ∈ D(Ω) which satisfies ∆f = 0 in D(Ω) is indeed a harmonic function (Folland, 1995, Coroll. 2.20), and hence belongs to C ∞ (Ω).

Harmonic functions satisfy the so-called maximum principle (and an analogous property for the minimum), which imposes strong limitations on the location of extrema of such functions, and thus helps us to ensure that our optimization method does not get trapped into local maxima. Two useful versions of the maximum principle are the following (see [Axler et al., 2001, Coroll. 1.9] or [Folland, 1995, Coroll. 2.14] for the former and [Axler et al., 2001, Claim 1.29] for the latter). (3)

Under some regularity assumptions on the border ∂Ω of the domain and on the function g describing the imposed border conditions, there exists a unique solution of (3) which is continuous on Ω. For some particular domains, such a solution can be characterized in the form of an integral, involving the values of g on ∂Ω and a function (called "Poisson kernel") depending on the shape of the domain. When the domain is the unit ball centered at the origin, i.e., Ω = B 1 (0) := { z : z < 1 }, the Poisson kernel is the following (Axler et al., 2001, eq. 1.15)

P B1(0) (z, ζ) := 1 -z 2 z -ζ 2 , z ∈ B 1 (0) , ζ ∈ ∂B 1 (0) .
This kernel allows to write the solution of the Dirichlet problem on the unit ball, as follows. Proposition 3. (Axler et al., 2001, Thm. 1.17) Given a continuous function g : ∂B 1 (0) → R, define the function

f : B 1 (0) → R as follows f (z) :=      1 2π 2π 0 P B1(0) (z, e iθ )f (e iθ ) dθ z ∈ B 1 (0) g(z) z ∈ ∂B 1 (0) .
Then, f is continuous on B 1 (0) and is harmonic on B 1 (0).

Moreover, the solution of the Dirichlet problem is unique, and hence, for any function f which is harmonic on the unit ball and continuous on its closure, f (z) at points inside the ball can be computed with a formula involving only the values of the restriction of f to the border ∂B 1 (0). Proposition 4. (Axler et al., 2001, Thm. 1.21) Let f : B 1 (0) → R be continuous on B 1 (0) and harmonic on B 1 (0). Then, for all z ∈ B 1 (0),

f (z) = 1 2π 2π 0 P B1(0) (z, e iθ )f (e iθ ) dθ .
By a simple dilation and translation of coordinates (mapping z to (z -z 0 )/R) it is possible to obtain an analogous formula for the ball B R (z 0 ) := { z : z -z 0 < R } (see e.g. [Folland, 1995, Chapt. 2, Exercise 1]). Proposition 5.

(Poisson integral formula for B R (z 0 )) Let f : B R (z 0 ) → R be continuous on B R (z 0 ) and harmonic on B R (z 0 ). Then, for all z ∈ B R (z 0 ), f (z) = 1 2π 2π 0 P B R (z0) (z, z 0 + Re iθ )f (z 0 + Re iθ ) dθ ,
where P B R (z0) (z, ζ) is the Poisson kernel for B R (z 0 ), defined as follows:

P B R (z0) (z, ζ) := R 2 -z -z 0 2 z -ζ 2 z ∈ B R (z 0 ) , ζ ∈ ∂B R (z 0 ) .

Gradient computation using Poisson integral formula

Consider an open set Ω ⊆ R 2 and a function f harmonic on Ω. Poisson integral formula given in Proposition 5 can be applied to any ball B R (z 0 ) such that its closure B R (z 0 ) is contained in Ω, because this ensures that f is harmonic and continuous in B R (z 0 ); this allows to compute the value of f (z) at points z inside the ball by using measurements of f along the circle ∂B R (z 0 ). The Poisson integral formula also gives a technique to compute derivatives (gradient, Hessian etc.) of f at any point inside a ball B R (z 0 ), with an integral involving only the values of f along the circle ∂B R (z 0 ), as follows. Proposition 6.

Let Ω ⊆ R 2 be an open set, and f : Ω → R be harmonic on Ω. For any z 0 ∈ Ω, for any R > 0 such that B R (z 0 ) ⊆ Ω, for any z = (x, y) ∈ B R (z 0 ), and for any non-negative integers m, n,

∂ m+n f (z) ∂x m ∂y n = 1 2π 2π 0 ∂ m+n ∂x m ∂y n P B R (z0) (z, z 0 +Re iθ )f (z 0 +Re iθ ) dθ .
Proof : The assumptions ensure that f is harmonic and continuous on B R (z 0 ), so that the Poisson integral formula from Proposition 5 holds true. Then, it is immediate to notice that in the Poisson integral formula the only dependence on z is in the Poisson kernel, so that one can exchange integration (which is with respect to θ) and derivation (which is with respect to x and/or y).

In particular, the gradient can be computed as follows: for all z ∈ B R (z 0 ),

∇f (z) = 1 2π 2π 0 ∇ z P B R (z0) (z, z 0 + Re iθ )f (z 0 + Re iθ ) dθ,
where the gradient of the Poisson kernel with respect to z = (x, y) is the following

∇ z P B R (z0) (z, ζ) = 2 z -ζ 2 (z 0 -z)+ 2P B R (z0) (z, ζ) z -ζ 2 (ζ -z).
Remark 1. The integral formulas given in Propositions 5 and 6 become significantly simpler when z = z 0 , the center of the circle B R (z 0 ) on which the integral is computed; for instance, one can obtain the following formula for f (z 0 ) (also known as the Mean Value Theorem)

f (z 0 ) = 1 2π 2π 0 f (z 0 + Re iθ ) dθ
and the following expressions for the gradient and for the Hessian matrix, respectively

∇f (z 0 ) = 1 πR 2π 0 f (z 0 + Re iθ ) cos θ sin θ dθ , (4)
H(z 0 ) = 2 πR 2 2π 0 f (z 0 + Re iθ ) cos 2θ sin 2θ sin 2θ -cos 2θ dθ .
(5)

Remark 2. A discrete approximation of Eq. ( 4) had already been used by [START_REF] Moore | Source seeking via collaborative measurements by a circular formation of agents[END_REF][START_REF] Moore | Source seeking via collaborative measurements by a circular formation of agents[END_REF][START_REF] Briñón Arranz | Collaborative estimation of gradient direction by a formation of AUVs under communication constraints[END_REF], for source seeking in the context of underwater vehicles (where each vehicle is assumed to be rotating over a circle); in both works only very mild assumptions where done on the function f , and the use of the integral formula (4) was justified by the fact that, for every

f which is C 1 in a neighborhood of z 0 , lim R→0 1 πR 2π 0 f (z 0 + Re iθ ) cos θ sin θ dθ = ∇f (z 0 )
and, for every integer

N ≥ 3, lim R→0 2 RN N j=1 f (z 0 + Re i 2jπ N ) cos 2jπ N sin 2jπ N = ∇f (z 0 ) ,
as it was proved in Briñón Arranz et al., 2011, Lemma 1 using the Taylor expansion of f and trigonometric properties. Our paper provides a framework in which the integral formula for the gradient computation is indeed exact for any given R > 0 (provided that the ball B R (z 0 ) is contained in Ω), and not only in the limit for R → 0; a non-vanishing radius is necessary in order to attenuate the effect of quantization and measurement noise. Other advantages of our approach are that the Poisson integral formula allows to obtain higher order derivatives, and allows a generalization to higher dimension (3d applications will be the object of our future work).

SOURCE-SEEKING STRATEGY

In this section we describe a method for solving the sourceseeking problem described in Section 2. The main idea is to use gradient, and possible high-order derivatives, information obtained from the values sensed along a circle, thanks to the Poisson formulas presented in previous sections.

Mobile sensing devices

We consider several mobile robots, each endowed with one or more sensors providing pointwise concentration measures; we consider the problem where each robot is required to perform the source localization task. The presence of multiple robots can be useful to ensure redundancy, so as to protect against failures, and also to better describe the source boundary, in a scenario where each robot is able to find only one point on ∂Ω s .

The robots' dynamics can be modeled in various ways, depending on the application at hand: as an example, we consider robots modeled as a nonholonomic unicycle

ż(t) = ve iθ(t) θ(t) = u(t) (6)
controlled by the torsional torque u(t). z(t) describes the rotational point of the robot in the plane, and θ(t) is the heading angle; the heading velocity v is assumed here to be constant, but different strategies can be alternatively devised. Here, the central point is how to design the control u(t) using information from the sensors devices to reach the source origin. Ideally, we would like each robot to compute the gradient at its center position z(t) by using the Poisson formula (4), with z 0 = z(t); to this aim, it needs to collect measurements on a circle ∂B R (z(t)) and to compute the Poisson integral. Two practical ways to obtain good approximations are the following.

A first effective sensing device (depicted on the left in Figure 2) consists in N sensors, arranged along a circular array of radius R, centered at the robot's central position z(t) with equispaced angles, namely, the j th sensor is at position z j (t) = z(t) + Re i 2jπ N ; then, the integral is approximated by the Riemann sum using the N measured values. Notice that, thanks to the rotational invariance of the Laplacian, it is not crucial to specify whether the array of sensors maintains an orientation fixed with respect to the absolute coordinate system, or it is solidly connected to the robot and it rotates with the robot heading angle (in the latter case, the robot will compute the gradient in the local coordinate system). We denote the measurement of the j th sensor at time t by f j (z(t)) = f (z j (t)) + w j (t) , where w j (t) is the measurement noise; with this notation, the approximated version of the Poisson formula (4) is

∇f (z(t)) = 2 N R N j=1 f j (z(t)) cos 2jπ N sin 2jπ N . (7)
Similarly, the approximated Hessian can be obtained from Eq. ( 5) as follows:

Hf (z(t)) = 4 N R 2 N j=1 f j (z(t)) cos 4jπ N sin 4jπ N sin 4jπ N -cos 4jπ N . (8)
A second setup (depicted on the right in Figure 2) is obtained by considering only one sensor instead of N , but allowing for a rotation of the sensor around the center of the robot; in this case, either it is supposed that the robot stays still during such a rotation, so that the integrals in Eqs. ( 4) and ( 5) are perfecly computed (apart from the measurement noise), or the robot moves during the rotation but with a speed sufficiently slow with respect to the rotation of the sensor, so that only a small error is introduced due to the deviation from the perfect circle.

In this paper we focus on the first setup.

Source-localization feedback design

The main idea is to perform a gradient ascent, with the gradient being computed by Eq. ( 7). In general, a harmonic function f might not be convex; however, the local maximum principle (see Proposition 2) ensures that f does not have any local maximum inside Ω: hence, search is ensured not to get trapped in any local maximum, except possibly on the outer boundary ∂Ω ext . The termination of the search on a local maximum on ∂Ω ext can be avoided by introducing some simple rule that allows to distinguish the external boundary from the internal one; for instance, one might have a knowledge of a rough lower bound on the value of f at the source, which is also an upper bound for values on the external boundary.

The gradient ascent strategy can be implemented by defining a reference heading θ r (z) in the direction of the gradient of the source at the point z, using formula (7). As the system is second-order, it requires a damping term involving time derivatives of both the heading angle and its reference; time derivatives of θ r can be computed from the approximated gradient (7) and Hessian (8).

The proposed feedback is

       u(t) = k P [θ r (z) -θ(t)] + k D [ θr (z) -θ(t)] , θ r (z) = arg ∇f (z) , θr (z) = v cos θ(t), sin θ(t) Hf (z) 0 -1 1 0 ∇f (z) ∇f (z) 2 , ( 9
)
where ∇f (z(t)) and Hf (z(t)) are defined by Eqs. ( 7) and (8).

Clearly, the reference heading angle θ r (z) is an approximation of the gradient's angle θ g (z) = arg ∇f (z). Moreover, the expression for θr (z) is an approximation of θg (z) = Here, ( ẋ, ẏ) = v cos θ(t) , sin θ(t) by ( 6), while the spatial derivatives of θ g (z) are computed from the expression θ g = arctan ∂f ∂y / ∂f ∂x , as follows

∂θg (z) ∂x ∂θg (z) ∂y = 1 ∇f (z) 2 ∂ 2 f (z) ∂y∂x ∂f (z) ∂x - ∂f (z) ∂y ∂ 2 f (z) ∂ 2 x ∂ 2 f (z) ∂ 2 y ∂f (z) ∂x - ∂f (z) ∂y ∂ 2 f (z) ∂x∂y = Hf (z) 0 -1 1 0 ∇f (z) ∇f (z) 2 .

Simulation results

We consider the problem of locating a heat source in a room, described in Example 2. The agents are robots such as those depicted on the left in Fig. 2, each having N = 12 sensors arranged on a circle with radius R = 20 cm, at equally spaced angles of 30°; we consider the motion and control laws given by Eqs. ( 6) and ( 9), where the constant velocity is chosen as v = 0.1 m/s. We suppose that the measurements are affected by white Gaussian noise of standard deviation σ = 0.5.

Figure 3 shows the trajectories of a set of robots starting from different locations: all agents reach the source, and their trajectories are perpendicular to the contour lines of the temperature; the small dithering in the robot trajectories is due to the noise in the measurements. The initial orientation is chosen as θ(0) = arctan(y 0 /x 0 ), and the initial angular velocity is θ(0) = 0; the values for the control constants are chosen as k P = 100 and k D = 20.

CONCLUSION

In this paper, we have addressed the problem of localization of a source, with mobile robots endowed with sensors providing pointwise concentration measures. We have focused our attention on sources that can be modeled by the Laplacian pde, which describes the steady state of homogeneous isotropic diffusion or heat; we have exploited the properties of the solutions of such equation (without making use of any explicit expression for the solution itself) in order to find formulas to compute the gradient of the measured quantity, as well as higherorder derivatives, using pointwise measurements along a circle. This allows the robots to implement a control law which drives them towards the source, either by a simple gradient-ascent strategy, or by some higher-order control law involving other derivatives. As an illustrative example, we have considered the search of a heat source in a room, with a robot modeled as a unicycle whose torque angle is controlled with a law involving the gradient and Hessian of the temperature.

  Example 1. (Heater in a circular room with constant Dirichlet boundary conditions) Consider a circular room with a circular heater at its center c. Denoting by B r (c) the open ball B r (c) := { z : z -c < r }, we have that the room is Ω = B r (c) and the heater location is

Example 2 .

 2 Fig. 1. Steady-state temperature in the room from Example 2

  Proposition 1. (Maximum principle) Let Ω ⊂ R 2 be a bounded open set, and let f : Ω → R be continuous in Ω and harmonic in Ω. Then, the maximum value of f on Ω is achieved on ∂Ω. Proposition 2. (Local maximum principle) Let Ω ⊆ R 2 be a connected open set, and let f : Ω → R be harmonic on Ω. If f has a local maximum in Ω, then f is constant. Consider now the Dirichlet problem with homogeneous boundary condition ∆f (z) = 0 in Ω f (z) = g(z) on ∂Ω .

Fig. 2 .

 2 Fig. 2. Possible sensor deployment over a robot. Left figure: N sensors are fixed along a circular path at the top of the robot, with equispaced angles; right figure: a single sensor is rotating fast on the top of the robot.

Fig. 3 .

 3 Fig. 3. Trajectories of the heat-seeking agents described in sect. 4.3.


Romain Benassi

Julien Bect

Emmanuel Vazquez

Bayesian optimization using sequential Monte Carlo

come

Bayesian optimization using sequential Monte Carlo

Romain Benassi, Julien Bect, and Emmanuel Vazquez SUPELEC, Gif-sur-Yvette, France Abstract. We consider the problem of optimizing a real-valued continuous function f using a Bayesian approach, where the evaluations of f are chosen sequentially by combining prior information about f , which is described by a random process model, and past evaluation results. The main difficulty with this approach is to be able to compute the posterior distributions of quantities of interest which are used to choose evaluation points. In this article, we decide to use a Sequential Monte Carlo (SMC) approach.

Overview of the contribution proposed

We consider the problem of finding the global maxima of a function f : X → R, where X ⊂ R d is assumed bounded, using the expected improvement (EI) criterion [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. Many examples in the literature show that the EI algorithm is particularly interesting for dealing with the optimization of functions which are expensive to evaluate, as is often the case in design and analysis of computer experiments [START_REF] Santner | The design and analysis of computer experiments[END_REF]. However, going from the general framework expressed in [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF] to an actual computer implementation is a difficult issue.

The main idea of an EI-based algorithm is a Bayesian one: f is viewed as a sample path of a random process ξ defined on R d . For the sake of tractability, it is generally assumed that ξ has a Gaussian process distribution conditionally to a parameter θ ∈ Θ ⊆ R s , which tunes the mean and covariance functions of the process. Then, given a prior distribution π 0 on θ and some initial evaluation results ξ(X 1 ), . . . , ξ(X n0 ) at X 1 , . . . , X n0 , an (idealized) EI algorithm constructs a sequence of evaluations points X n0+1 , X n0+2 , . . . such that, for each n ≥ n 0 ,

X n+1 = argmax x∈X ρn := θ∈Θ ρ n (x; θ)dπ n (θ) , (1)
where π n stands for the posterior distribution of θ, conditional on the σ-algebra F n generated by X 1 , ξ(X 1 ), . . . , X n , ξ(X n ), and

ρ n (x; θ) := E n,θ ((ξ(X n+1 ) -M n ) + | X n+1 = x) is the EI at x given θ, with M n = ξ(X 0 ) ∨ • • • ∨ ξ(X n
) and E n,θ the conditional expectation given F n and θ. In practice, the computation of ρ n is easily carried out (see [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]) but the answers to the following two questions will probably have a direct impact on the performance and applicability of a particular implementation: a) How to deal with the integral in ρn ? b) How to deal with the maximization of ρn at each step?

We can safely say that most implementations-including the popular EGO algorithm [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]-deal with the first issue by using an empirical Bayes (or plug-in) approach, which consists in approximating π n by a Dirac mass at the maximum likelihood estimate of θ. A plug-in approach using maximum a posteriori estimation has been used in [START_REF] Lizotte | An experimental methodology for response surface optimization methods[END_REF]; fully Bayesian methods are more difficult to implement (see [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF] and references therein). Regarding the optimization of ρn at each step, several strategies have been proposed (see, e.g., [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Gramacy | Particle learning of Gaussian process models for sequential design and optimization[END_REF][START_REF] Bardenet | Surrogating the surrogate: accelerating Gaussianprocess-based global optimization with a mixture cross-entropy algorithm[END_REF][START_REF] Ginsbourger | DiceOptim: Kriging-based optimization for computer experiments[END_REF]).

This article addresses both questions simultaneously, using a sequential Monte Carlo (SMC) approach [START_REF] Chopin | A sequential particle filter method for static models[END_REF][START_REF] Del Moral | Sequential Monte Carlo samplers[END_REF] and taking particular care to control the numerical complexity of the algorithm. The main ideas are the following. First, as in [START_REF] Gramacy | Particle learning of Gaussian process models for sequential design and optimization[END_REF], a weighted sample

T n = {(θ n,i , w n,i ) ∈ Θ × R, 1 ≤ i ≤ I} from π n is used to approximate ρn ; that is, I i=1 w n,i ρ n (x; θ n,i ) → I ρn (x)
. Besides, at each step n, we attach to each θ n,i a (small) population of candidate evaluation points {x n,i,j , 1 ≤ j ≤ J} which is expected to cover promising regions for that particular value of θ and such that max i,j ρn (x n,i,j ) ≈ max x ρn (x).

Algorithm and results

At each step n ≥ n 0 of the algorithm, our objective is to construct a set of weighted particles

G n = γ n,i,j , w ′ n,i,j , γ n,i,j = (θ n,i , x n,i,j ) ∈ Θ × X, w ′ n,i,j ∈ R , 1 ≤ i ≤ I, 1 ≤ j ≤ J (2)
so that i,j w ′ n,i,j δ γn,i,j → I,J π ′ n , with

dπ ′ n (γ) = gn (x | θ) dλ(x) dπ n (θ) , x ∈ X , θ ∈ Θ , γ = (θ, x),
where λ denotes the Lebesgue measure, gn

(x | θ) = g n (x | θ)/c n (θ), g n (x | θ
) is a criterion that reflects the interest of evaluating at x (given θ and past evaluation results), and c n (θ) = X g n (x | θ)dx is a normalizing term. For instance, a relevant choice for g n is to consider the probability that ξ exceeds M n at x, at step n. (Note that we consider less θs than xs in G n to keep the numerical complexity of the algorithm low.)

To initialize the algorithm, generate a weighted sample T n0 = {(θ n0,i , w n0,i ), 1 ≤ i ≤ I} from the distribution π n0 , using for instance importance sampling with π 0 as the instrumental distribution, and pick a density q n0 over X (the uniform density, for example). Then, for each n ≥ n 0 :

Step 1: demarginalize -Using T n and q n , construct a weighted sample G n of the form (2), with x n,i,j

iid ∼ q n , w ′ n,i,j = w n,i gn(xn,i,j|θn,i) qn(xn,i,j )cn,i , and c n,i =

J j ′ =1 gn(x n,i,j ′ |θn,i) qn(x n,i,j ′ ) .
Step 2: evaluate -Evaluate ξ at X n+1 = argmax i,j

I i ′ =1 w n,i ′ ρ n (x n,i,j ; θ n,i ′ ).
Step 3: reweight/resample/move -Construct T n+1 from T n as in [START_REF] Chopin | A sequential particle filter method for static models[END_REF]: reweight the θ n,i s using w n+1,i ∝ πn+1(θn,i) πn(θn,i) w n,i , resample (e.g., by multinomial resampling), and move the θ n,i s to get θ n+1,i s using an independant Metropolis-Hastings kernel.

Step 4: forge q n+1 -Form an estimate q n+1 of the second marginal of π ′ n from the weighted sample X n = {(x n,i,j , w ′ n,i,j ), 1 ≤ i ≤ I, 1 ≤ j ≤ J}. Hopefully, such a choice of q n+1 will provide a good instrumental density for the next demarginalization step. Any (parametric or non-parametric) density estimator can be used, as long as it is easy to sample from; in this paper, a tree-based histogram estimator is used.

Nota bene: when possible, some components of θ are integrated out analytically in (1) instead of being sampled from; see [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF].

Experiments. Preliminary numerical results, showing the relevance of a fully Bayesian approach with respect to empirical Bayes approach, have been provided in [START_REF] Benassi | Robust Gaussian process-based global optimization using a fully Bayesian expected improvement criterion[END_REF]. The scope of these results, however, was limited by a rather simplistic implementation (involving a quadrature approximation for ρn and a non-adaptive grid-based optimization for the choice of X n+1 ). We present here some results that demonstrate the capability of our new SMC-based algorithm to overcome these limitations.

The experimental setup is as follows. We compare our SMC-based algorithm, with I = J = 100, to an EI algorithm in which: 1) we fix θ (at a "good" value obtained using maximum likelihood estimation on a large dataset); 2) X n+1 is obtained by exhaustive search on a fixed LHS of size I × J. In both cases, we consider a Gaussian process ξ with a constant but unknown mean function (with a uniform distribution on R) and an anisotropic Matérn covariance function with regularity parameter ν = 5/2. Moreover, for the SMC approach, the variance parameter of the Matérn covariance function is integrated out using a Jeffreys prior and the range parameters are endowed with independent lognormal priors.

Results. Figures 1(a) and 1(b) show the average error over 100 runs of both algorithms, for the Branin function (d = 2) and the log-transformed Hartmann 6 function (d = 6). For the Branin function, the reference algorithm performs better on the first iterations, probably thanks to the "hand-tuned" parameters, but soon stalls due to its non-adaptive search strategy. Our SMC-based algorithm, however, quickly catches up and eventually overtakes the reference algorithm. On the Hartmann 6 function, we observe that the reference algorithm always lags behind our new algorithm.

We have been able to find results of this kind for other test functions. These findings are promising and need to be further investigated in a more systematic large-scale benchmark study.

Fig. 1 .

 1 Fig. 1. A comparison of the average error to the maximum (100 runs)