Yassine Laguel

Jérôme Malick

Zaid Harchaoui

FIRST-ORDER OPTIMIZATION FOR SUPERQUANTILE-BASED SUPERVISED LEARNING

Keywords: supervised learning, risk measure, distributional robustness, nonsmooth optimization

Classical supervised learning via empirical risk (or negative log-likelihood) minimization hinges upon the assumption that the testing distribution coincides with the training distribution. This assumption can be challenged in modern applications of machine learning in which learning machines may operate at prediction time with testing data whose distribution departs from the one of the training data. We revisit the superquantile regression method by proposing a first-order optimization algorithm to minimize a superquantile-based learning objective. The proposed algorithm is based on smoothing the superquantile function by infimal convolution. Promising numerical results illustrate the interest of the approach towards safer supervised learning.

INTRODUCTION

Classical supervised learning assumes that, at training time, we have access to examples (x 1 , y 1 ), . . . , (x n , y n ) drawn i.i.d. from a distribution P, and that at testing time, we may face a new example, also drawn from P. The learned predictor or function can be used by humans or machines to make decisions, or used in as an intermediate component in a greater data processing and computing system. This common framework is currently challenged by important domain applications [START_REF] Recht | Do imagenet classifiers generalize to imagenet?[END_REF], in which several of the standard assumptions turn out to be unrealistic or simply incorrect. We may not face the same distribution at test time as we did at training time (train-test distribution shift). Recent failures of learning systems when operating in unknown environments [START_REF] Metz | Microsoft's neo-Nazi sexbot was a great lesson for makers of AI assistants[END_REF][START_REF] Knight | A self-driving Uber has killed a pedestrian in Arizona[END_REF] underscore the importance of reconsidering the learning objective used to train learning machines in order to ensure robust behavior in the face of unexpected distributions at prediction time.

The generalized regression framework presented in [START_REF] Rockafellar | Risk tuning with generalized linear regression[END_REF] provides an attractive ground to design learning machines dis-playing increased robustness in the face of unexpected testing distributions. The framework hinges upon the notion of superquantile, a statistical summary of a distribution tail [START_REF] Lee | Minimax statistical learning with Wasserstein distances[END_REF][START_REF] Duchi | Variance-based Regularization with Convex Objectives[END_REF][START_REF] Kuhn | Wasserstein distributionally robust optimization: Theory and applications in machine learning[END_REF]. This notion of robustness is aligned with the one in distributionally robust optimization [START_REF] Ben-Tal | Robust optimization[END_REF] and empirical likelihood estimation [START_REF] Owen | Empirical Likelihood, Chapman & Hall/CRC Monographs on Statistics & Applied Probability[END_REF]. It is, however, different, from notions of robustness commonly considered in robust statistics [START_REF] Ben-Tal | Robust optimization[END_REF]Sec. 12.6].

The superquantile is a risk measure, a family of statistical summaries of distribution tails, well studied in economics and finance [START_REF] Rockafellar | Optimization of Conditional Value-at-Risk[END_REF][START_REF] Ben-Tal | An old-new concept of convex risk measures: The optimized certainty equivalent[END_REF]. The quantity is, however, a nonsmooth function. We present here a simple approach, based on infimal convolution smoothing, which allows one to easily adapt state-of-the-art gradient-based optimization algorithms for classical supervised learning to the superquantilebased learning framework. Moreover, we provide a companion software package in Python available here https: //github.com/yassine-laguel/spqr .

Superquantile

Risk measures play a crucial role in optimization under uncertainty, involving problems with an aversion to worst-cases scenarios. Among popular convex risk measures, superquantile (also called Conditional Value at Risk) has received a special attention because of its nice convexity properties; see e.g. the textbook [START_REF] Shapiro | Lectures on stochastic programming: modeling and theory[END_REF]Chap. 6].

We use here the notation and terminology of Rockafellar and Royset [START_REF] Rockafellar | Superquantiles and their applications to risk, random variables, and regression[END_REF]. The p-quantile Q p (U ) of a random variable U is defined as the general inverse of the cumulative distribution of U . More precisely, for a random variable U (admitting a second order moment), the cumulative distribution function F U : R → [0, 1] is defined as F U (x) = P(U ≤ x). For any p ∈ [0, 1], the p-quantile Q p (U ) and the p-superquantile Qp (U ), are respectively defined by

Q p (U ) = min{x ∈ R, F U (x) ≥ p} Qp (U ) = 1 1 -p 1 p =p Q p (U )dp . ( 1
)
The superquantile is, therefore, a measure of the upper tail.

The parameter p allows one to control the sensitivity to risk.

The superquantile enjoys a dual representation [START_REF] Föllmer | Convex measures of risk and trading constraints[END_REF] Qp (U ) = max

0≤q(•)≤ 1 1-p Ω q d P(ν)=1 ν∈Ω U (ν)q(ν)dP(ν) (2)
Interestingly, the dual formulation uncovers another interpretation of the superquantile learning objective relating it to the re-weighting of the terms in the empirical risk. In practice, the ambiguity on the data distribution may be formalized before training, for instance by incorporating side information (geographical and/or temporal for instance) that drives the heterogeneity of the data. Superquantile learning is expected to produce models that perform better in case of distributional shifts between the training time and the testing time, compared to models trained using standard empirical risk minimization.

Superquantile-based learning

We are interested in a supervised machine learning setting with training data D = (x i , y i ) 1≤i≤n ∈ (R p × R q ) n , a prediction function ϕ : R d × R p → R q (such as a linear model or a neural network) and a loss function : R q × R q → R (such as the logistic loss or the least-squares loss). The classical empirical risk minimization writes

min w∈R d E (xi,yi)∼D ( (y i , ϕ(w, x i )) . (3)
A natural approach consists then in replacing the expectation in (3) by the superquantile (1) in the case of discrete distributions standing for the training data

min w∈R d [ Qp ] (xi,yi)∼D (y i , ϕ(w, x i ) (4)
Introducing L i (w) = (y i , ϕ(w, x i )) and L(w) = (L i (w)) i , we simply write the superquantile optimization problem as

min w∈R d f (w) = Qp (L(w)). (5)
Note that the objective function can specified by expressing the superquantile in its dual formulation (2) for the discrete distribution

f (w) = sup q∈Kp n i=1 q i L i (w) with K p = q ∈ R n , n i=1 q i = 1, q i ∈ 0, 1 n(1 -p) ∀i .
This representation is central to the implementation as we shall see in Sec. 3. Existing works on minimizing superquantiles considered linear programming or convex programming including interior point algorithms; see [START_REF] Rockafellar | Superquantile regression with applications to buffered reliability, uncertainty quantification, and conditional value-at-risk[END_REF]. Our approach considers first-order algorithms instead; although natural, this work seems to be the first one to do so.

SMOOTHING THE SUPERQUANTILE

In this section, we study the differentiability properties of the superquantile objective [START_REF] Lee | Minimax statistical learning with Wasserstein distances[END_REF]. We first derive the expression of the subdifferential, when the L i -s are convex1 . Then, when L i are smooth (and possibly nonconvex), we show how to smooth the superquantile by infimal convolution in order to apply gradient-based optimization algorithms [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF].

Subdifferential expression

The superquantile risk measure ( 5) is usually nonsmooth and computing its subdifferential (or even a single subgradient) is not straightforward. Using the dual reformulation (2), we get the expression of the entire subdifferential for the convex case. Note that gradients of superquantile-based functions for general distributions are obtained, with advanced tools, in [START_REF] Ruszczyński | Optimization of convex risk functions[END_REF]. Interestingly, the nonsmoothness of these functions arises only with discrete distributions.

Proposition 2.1. Assume the model ϕ and the loss are such that the L i are convex. For w ∈ R d , let I p (w) be the set of indices i ∈ {1, . . . , n} such that L i (w) = Q p (L(w)). Then the subdifferential reads as a Minkowski sum

∂f (w) = 1 1 -p i∈{1,...n} L i (w)>Qp(L(w)) ∂L i (w) n + 1 1 -p i∈Ip(w) α i ∂L i (w) n , α i ∈ [0, 1] ∀i ∈ I p (w) 1 n i∈Ip(w) α i = 1 n n i=1 1 L i (w)≤Qp(L(w)) -p
In particular, when L is differentiable at w, f is differentiable at w if and only if the set I p (w) is reduced to a singleton.

Proof. The proof consists in applying various convex calculus rules, taken from the textbook [18, Chap D]. First we apply Theorems 4.1.1 and 4.4.

2 to h i (w, η) = max(L i (w) -η) ∂h i (x, η) = {(∂L i (w), -1)(1 L i (w)>η +α1 L i (w)=η ), α ∈ [0, 1]} We apply Theorem 4.1.1 with h(w, η) = η+ 1 n(1-p) n i=1 h i (w, η) ∂h(w, η) = 1 1 -p n i=1 ∂L i (w) n (1 L(w)i>η + α i 1 L i (w)=η ), 1 - 1 1 -p n i=1 1 n (1 L i (w)>η + α i 1 L i (w)=η ) , α i ∈ [0, 1], ∀i ∈ {1, . . . , n} .
By [START_REF] Rockafellar | Optimization of Conditional Value-at-Risk[END_REF], f satisfies f (w) = min η∈R h(w, η), with Q p (L(w)) = arg min η∈R h(w, η). We can thus apply Corollary 4.5.3 to get

∂f (w) = 1 1 -p n i=1 ∂L i (w) n δ i (w, α) with α s.t. 0 = 1 - 1 1 -p n i=1 δ i (w, α) n and α i ∈ [0, 1], ∀i with δ i (w, α) = (1 L i (w)>Qp(L(w)) + α i 1 L i (w)=Qp(L(w)) ).
Observe finally that for any sequence

(α i ) 1≤i≤n 0 = 1- 1 1 -p n i=1 δ i (w, α) n ⇔ 1 n i∈Ip(z) α i = 1 -p - n i=1 1 n 1 L i (w)>Qp(L(w)) ⇔ 1 n i∈Ip(z) α i = 1 -p -(1 -P[L(w) ≤ Q p (L(w))]) ⇔ 1 n i∈Ip(z) α i = 1 n n i=1 1 L i (w)≤Qp(L(w)) -p
which yields the result.

Thus, the computation of a subgradient can be performed in linear time: the cost essentially stems from the computation of the quantile Q p (L(w)) and the sum of vectors in R d (assuming such sums can be computed in constant time).

Gradient of smoothed approximation

As shown in Proposition 2.1, the objective function is not differentiable in general (even when L is differentiable), and we propose to smooth it using infimal convolution as in [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. More precisely, we follow the methodology of [START_REF] Beck | Smoothing and first order methods: A unified framework[END_REF] and we propose to smooth only the superquantile Qp rather than the whole function f . Given formulation (2), we introduce

f µ (w) = max q∈Kp n i=1 q i L i (w) -µ d(q) for µ > 0 (6)
where d : R n → R is a fixed non-negative strongly convex function that satisfies min q∈K d(q) = 0. In this paper, we consider the euclidean distance to the uniform probability measure and the entropic penalty function

d(q) = 1 2 q - 1 n e 2 and d(q) = log(n)+ n i=1 q i log(q i )
where e = (1, . . . , 1) is the usual vectors of all ones. As a direct application of [16, Th. 1], we have the following proposition establishing that f µ is a smooth approximation of f . Algorithm 1: Fast subroutine for smoothed oracle

1 Initialization. u = L(w) + µ n e, = 1 n(1-p) , q µ = 0 ∈ R n 2 P = {u i , i ∈ {1, . . . , n}} ∪ {u i -µ , i ∈ {1, . . . ,
9 for 1 ≤ k ≤ n do 10 if λ < u k -µ then 11 [q µ ] k = ; 12 else if u k -µ ≤ λ < u k then 13 [q µ ] k = u k -λ µ ; 14 else 15 [q µ ] k = 0 16 end 17 end Output: q µ ∈ R n : solution of (6)
Proposition 2.2 (Gradient of smoothed approximation). Assume the model ϕ and the loss are such that the L i are smooth for any i. In the above setting, the convex function f µ provides a global approximation of f , i.e.

f µ (w) ≤ f (w) ≤ f µ (w) + µ 2 for any w ∈ R d . If L is differentiable, then f µ is differentiable as well, with ∇f µ (w) = JL(w) T q µ (w), (7)
where JL(w) is the Jacobian of L at w and q µ (w) is the optimal solution of (6), unique by strong convexity of d.

To be made practical, the previous result needs to be equipped with a fast and efficient procedure to solve [START_REF] Duchi | Variance-based Regularization with Convex Objectives[END_REF]. As stated in the next proposition, Algorithm 1 addresses this issue. The procedures follows closely the ones in [START_REF] Condat | Fast projection onto the simplex and the l1 ball[END_REF], where convex duality and one-dimensional search ideas are fruitfully combined. Proposition 2.3. Algorithm 1 computes the optimal solution of the problem (6) (with the euclidean or the entropic penalty) at a cost of O(n) operations.

Proof. We detail the proof for d(q) = 1 2 q -1/n e 2 ; the second case of the entropy follows the same lines. We dualize the constraint n i=1 q i -1 = 0 to get the Lagrangian:

L(q, λ) = n i=1 q i L i (w)- µ 2 n i=1 q i - 1 n 2 +λ 1 - n i=1 q i .
With the notation and u introduced in the algorithm, the dual function writes:

θ(λ) = max q∈R n 0≤qi≤l L(q, λ) = λ- µ 2n + n i=1 max 0≤qi≤l (u i -λ)q i - µ 2 q 2 i
For λ ∈ R and i ∈ {1, . . . , n} fixed, let us introduce the function h i (q i ) = (u i -λ)q i -µ 2 q 2 i . Then, we get arg max

0≤qi≤l h i (q i ) =    0 if λ ≥ u i ui-λ µ if u i ≥ λ ≥ u i -µ if λ ≤ u i -µ (8)
As a result, we get the explicit expression of θ(λ). Observing that it is differentiable, we get

θ (λ) = 1 - n i=1 u i -λ µ 1 ui≥λ≥ui-µ + 1 ui-µ >λ .
Observe now that lim λ→+∞ θ (λ) = 1 and since n =

1 1-p > 1, lim λ→-∞ θ (λ) < 0.
Therefore, θ is a nondecreasing and continuous (piecewise affine) function that takes negative and positive values: by the intermediate value theorem, there exists a solution λ ∈ R such that θ (λ ) = 0. By duality theory, the associated q (the optimal solution of ( 8) for λ = λ ) is the solution of the primal problem [START_REF] Duchi | Variance-based Regularization with Convex Objectives[END_REF]. Finally, we compute λ zeroing θ . Since θ is piecewise affine, we just need to evaluate θ at points belonging to the set P and at a and b as defined in Algorithm 1. One can then find λ by testing three simple cases (i

) if θ (a) = 0, take λ * = a, if θ (b) = 0, take λ * = b, else, take λ * = a -θ (a)(b-a) θ (b)-θ (a)
. Regarding computational costs, this algorithm boils down to the search of a and b, and the assignment of the coordinates of q µ . This also sums up to a O(n) cost.

Thus Algorithm 1 provides an efficient oracle for minimizing of f with first-order algorithms.

A PYTHON TOOLBOX FOR SUPERQUANTILE OPTIMIZATION

We provide a Python software package called SPQR to the =community for research in superquantile-based optimization and learning. The software package includes optimization and modeling tools to solve problems of the form ( 5) with just a few lines of code. The implementation builds off basic structures of scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF].The code is publicly available at https://github.com/ yassine-laguel/spqr. We describe here the optimization methods used in in the toolbox and how to call the basic functions. We refer to the online documentation for more details, custom options, and parameter settings.

First-order optimization algorithms

Although stochastic gradient algorithms are popular methods to solve empirical risk minimization problems at scale (3), replacing the expectation by the superquantile in (4) completely changes the situation making these algorithms not directly applicable. Indeed computing the function values and gradients requires sorting loss values on the whole data set, which is not directly amenable to classical stochastic gradient algorithms. This rehabilitates batch optimization algorithms in our context. We cover a variety of methods • when L is convex: subgradient method and dual averaging. We implement in particular the "weighted" version of dual averaging with a Euclidean proxfunction [START_REF] Nesterov | Primal-dual subgradient methods for convex problems[END_REF]. For an iterate x k and a gradient g k of f at x k , the update writes:

x k+1 = -α k s k+1 with s k+1 = k i=0 g k g k
where (α k ) k≥0 denotes the tuned step-size of the method. The tuning is carried through a line-search strategy performed at the first iteration. To use these algorithms, we provide a subgradient oracle (from Proposition 2.1) with the same complexity as computing a quantile (ie. O(n) with n the number of data points).

• when L is smooth, we can use the smoothed objective: gradient method, accelerated gradient method and quasi-Newton (BFGS). In particular the accelerated gradient method relies on the following scheme [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF]:

α 0 = 0, α s = 1 + 1 + 4α s-1 ) 2 and γ s = 1 -α s α s+1
x s+1 = y s -1 β ∇f (y s ), y s+1 = (1 -γ s )x s+1 + γ s x s with x 0 = y 0 = 0. To use these algorithms, we provide a gradient oracle using Algorithm 1, again with a O(n) complexity (Proposition 2.3).

Basic usage: input format and execution

The user provides a dataset (X, Y ) ∈ R p × R m and an oracle for the function L and its gradient. The dataset is stored into two python lists (or numpy arrays) X and Y; for instance, for realizations of random variables:

import numpy as np X = np.random.rand(100, 2) alpha = np.array([1., 2.]) Y = np.dot(X, alpha) + np.random.rand(100)

The two python functions L and L prime are assumed to be functions of the triplet (w,x,y) where w is the optimization variable and (x,y) a data point. For instance, one can perform risk-sensitive linear regression with:

# Define the loss and its derivative def L(w,x,y): return 0. RiskOptimizer inherits from scikit-learn's estimators: we use the fit method to run the optimization algorithm on the data, providing a solution of [START_REF] Lee | Minimax statistical learning with Wasserstein distances[END_REF].

# Running the algorithm optimizer.fit(X,Y) lst_iterates = optimizer.list_iterates sol = optimizer.solution

NUMERICAL ILLUSTRATIONS

We compare the proposed approach (4) with the common approach using empirical risk minimization on synthetic and real data. We solve the ordinary least squares problem

min w∈R d E (xi,yi)∼D (y i -w x i ) 2
using the corresponding function of scikit-learn (by calling LinearRegression.fit(X,Y) method). We solve its risk-sensitive counterpart min w∈R d

[ Qp ] (xi,yi)∼D (y i -w x i ) 2

using our toolbox with risk-sensitive linear regression, Euclidean smoothing (with µ = 1000), and L-BGFS as optimizer (see Sec. 3).

Synthetic Dataset

We consider a regression task on a synthetic training dataset of n = 10 4 points in R 40 × R. The design matrix X = (x i ) 1≤i≤n is generated with the make low rank matrix procedure of scikit learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] with a rank 30. For a given model parameter w ∈ R, we generate the data according to

y i = x i w + ε i .
The noise ε i is defined here as a mixture

ε i = βε N + (1 -β)ε L
where all random variables are independent, ε N follows a standard normal distribution, ε L follows a Laplace distribution with location µ = 10 and scale s = 1, and β follows a Bernoulli distribution with parameter p = 0.8. Define the squared residuals (or losses)

r 2 i = (y i -w x i ) 2 for i = 1, . . . , n
and the p-quantiles of the empirical distribution of (r 2 i ) i=1,...,n for p = 0.5 and p = 0.9.

Model

Mean p-quantile of the loss p = 0.5 p = 0.9 We report the p-quantiles and the distribution of losses obtained on the training dataset. and on a test dataset of 2000 data points independently generated with the same procedure; see Table 1 and Figure 1. As p grows, the superquantile-based or risk-sensitive model shifts the upper tail on errors to the left, which shows an improved performance on extreme cases. This comes with the price of lower performances on inputs well managed by the standard approach (see metrics for p = 0.5 in Table 1).

E

Real Dataset

We consider the superconductivity dataset [START_REF] Hamidieh | A data-driven statistical model for predicting the critical temperature of a superconductor[END_REF] which contains the information of 21, 263 superconductors. The learning task is to predict the critical temperature of a superconductor from the 10 most important features as selected by [START_REF] Hamidieh | A data-driven statistical model for predicting the critical temperature of a superconductor[END_REF]. We split the dataset into a training set and a testing set with a ratio 80%/20%.

We report in Figure 2 the comparison between the quantiles of the testing and training loss distribution respectively. In terms of the quantile at 90%, the proposed approach display better statistical behavior on the testing loss than the common approach based on empirical risk minimization. This is in line with the aim of the formulation considered, which seeks to gain a better control on the tails of the loss distribution.

Model Mean

p-quantile of the loss p = 0.9 p = 0.95 p = 0.99

E

CONCLUSION

Risk-sensitive optimization plays a major role in the design of safer models for decision-making and has recently gained interest in machine learning. We provide a toolbox to tackle superquantile-based learning problems using first-order optimization algorithms. Numerical illustrations on regression tasks show an improved statistical behavior in terms of higher quantiles of the testing loss.

n}} 3 4 b 6 λ := a; 7 else 8 λ

 3468 Find a := max {s ∈ P, θ (s) ≤ 0} := min {s ∈ P, θ (s) > 0}; 5 if θ (a) = 0 then := a -θ (a)(b-a) θ (b)-θ (a)

Fig. 1 :

 1 Fig.1: Quantiles of the empirical distribution of residuals on the test. The risk-sensitive model was trained with p = 0.9.

Fig. 2 :

 2 Fig.2: Distribution of the loss values r i on the train and test superconductivity dataset. The risk-sensitive model is trained with p = 0.9.5. CONCLUSION

  Before solving the problem (5), we have to instantiate the RiskOptimizer object of SPQR with the two oracles, following standard usage of scikit-learn. The basic instantiation is as follows.

	def L_prime(w,x,y):
	return -1.0 * (y -np.dot(x,w)) * x
	from SPQR import RiskOptimizer
	# Instantiate a risk optimizer object
	optimizer = RiskOptimizer(L, L_prime)
	5 * np.linalg.norm(y -
	np.dot(x,w)) ** 2

Table 1 :

 1 Quantiles of the empirical distribution of residuals on the test.

	16.45	5.55	60.2
	Qp -p = 0.5 18.75	13.9	41.2
	Qp -p = 0.7 22.3	20.7	36.6
	Qp -p = 0.9 23.7	22.5	37.7

Table 2 :

 2 Metrics of the distribution of the loss values ri on the test

					16.5		35.8		42.7		55.7
		Qp -p = 0.8	17.4		34.7		41.0		53.8
		Qp -p = 0.9	18.1		35.6		41.0		53.6
		Qp -p = 0.95 18.9		36.5		41.4		53.6
	superconductivity dataset
	0.05 0.06		Training Loss Distribution Ordinary Least Squares Risk-averse Regression	0.025 0.030		Test Loss Distribution Ordinary Least Squares Risk-averse Regression
	0.04						0.020
	0.03						0.015
	0.02						0.010
	0.01						0.005
	0.00	40	50	60	70	80	0.000	40	50	60	70	80

The authors gratefully acknowledge support from NSF CCF 1740551, DMS 1839371, and faculty research awards.

978-1-7281-6662-9/20/$31.00 c 2020 IEEE

Convexity of the L i -s is guaranteed when e.g. the model ϕ is linear and the loss is convex with respect to its second variable, as for the l

-squared loss and the the cross-entropy loss.


Alessandro Renzaglia

Jilles Dibangoye

Vincent Le Doze

Olivier Simonin

A Common Optimization Framework for Multi-Robot Exploration and Coverage in 3D Environments

Keywords: Multi-Robot Systems, Cooperative Exploration, Optimal Coverage

This paper studies the problems of static coverage and autonomous exploration of unknown three-dimensional environments with a team of cooperating aerial vehicles. Although these tasks are usually considered separately in the literature, we propose a common framework where both problems are formulated as the maximization of online acquired information via the definition of single-robot optimization functions, which differs only slightly in the two cases to take into account the static and dynamic nature of coverage and exploration respectively. A common derivative-free approach based on a stochastic approximation of these functions and their successive optimization is proposed, resulting in a fast and decentralized solution. The locality of this methodology limits however this solution to have local optimality guarantees and specific additional layers are proposed for the two problems to improve the final performance. Specifically, a Voronoi-based initialization step is added for the coverage problem and a combination with a frontier-based approach is proposed for the exploration case. The resulting algorithms are finally tested in simulations and compared with possible alternatives.

Introduction

Multi-robot teams, especially when involving aerial vehicles, are extremely efficient systems to help humans in acquiring information on large and complex areas [START_REF] Scaramuzza | Vision-controlled micro flying robots: from system design to autonomous navigation and mapping in GPS-denied environments[END_REF][START_REF] Erdelj | Help from the sky: Leveraging uavs for disaster management[END_REF]. In these scenarios, two fundamental tasks are static coverage and exploration. In the coverage problem, the robots have to find the static set of positions that optimizes a certain coverage criterion, e.g. the portion of the Fig. 1 Mission goal: a team of cooperating UAVs has to gather information on an unknown 3D environment (gray surface), guided by the observations collected online (blue region). Depending on the mission, the final goal is either to explore the environment or to statically cover it. monitored environment. In the exploration, the objective is to generate the paths that allow the robots to observe the entire environment in a minimum time. In both cases, the region of interest is assumed to be unknown and the robots have to retrieve the necessary information during the mission (Fig. 1). For their relevance, these problems have received wide attention in the robotics community and numerous solutions have been proposed in the last years. However, they have been usually treated separately, proposing different formulations and approaches, and not as particular cases of a more general information-based problem. In both cases, the robots are indeed called to navigate through an unknown environment and cooperate to maximize the observed area.

The main contributions of this paper are: to present a common framework to formulate visual coverage and exploration, and to design two variants of a decentralized optimization-based solution for the two considered problems. The principal advantage of formulating both problems in such a unified framework is the transfer of knowledge, e.g., theories, methods and insights can be more easily extended from one setting to the other. In this paper, we make use of an observation-based optimization problem solved employing local stochastic approximation of a suitable optimization function which was initially introduced for the coverage problem as a building block of our solution for the exploration task. More specifically, the two problems are presented as optimization problems where, based on collected data, the robots try to online maximize the portion of observed environment. To do so, we define singlerobot optimization functions whose maximization leads to a decentralized solution of the considered problems. A common optimization framework based on stochastic approximations, exploiting the data acquired during the task and exchanged among neighbors, is proposed to efficiently obtain a solution and respect the constraints of the problem. However, as any local approach, if no prior information on the environment is exploited, the solution suffers from local optima limitations. To overcome this problem, we propose two specific solutions for the coverage and exploration problems.

Local minima in a coverage task make the team converge to configurations that can be significantly far from the global optimum of the problem. These solutions are strongly dependant on the initial state of the team and a way to overcome this drawback can be found in initializing the optimization with a favorable configuration. To this end, the presence of a partial knowledge on the environment can be an important source of information to exploit. This kind of information is available in many real situations, where a coarse knowledge of the shape, or of few specific features, is available or can be gathered before the mission starts. Based on this information, we propose a Voronoi-based partition of an approximated version of the environment to cover to obtain an initial configuration from which the local online optimization can start. On the other hand, in the exploration case, the optimization can be locally stuck and the robots would need external global inputs to identify new directions to pursue and be sure not to leave any part of the environment unexplored. To this end, we present a solution that combines the stochastic optimization with a frontier-based approach able to provide to each robot this global information whenever required to ensure the completeness of the exploration. This solution has been recently presented in [START_REF] Renzaglia | Combining stochastic optimization and frontiers for aerial multi-robot exploration of 3d terrains[END_REF].

The rest of the paper is organized as follows. The next section reviews the related work for both the considered problems. Section 3 presents the common formulation and the definition of the two optimization functions. In Section 4, the stochastic approximation based optimization method adopted in this work is presented. The two specific solutions to overcome the limitations due to the locality of the proposed approach are then described in Section 5. The resulting algorithms are finally tested in simulation and the results for the coverage and exploration problems are provided in Sections 6 and 7 respectively.

Related Work

In this section, we provide an overview of the most relevant related literature on static coverage and autonomous exploration, with a particular focus on the works taking into account three-dimensional environments and involving aerial vehicles. Since these two problems have been usually studied separately in the literature, we here present two distinct discussions on the respective related works and their connections with our contribution.

Static Coverage

The optimal coverage of a given area with a team of mobile sensors is a widely studied problem in literature and several formulations have been considered over the years. [START_REF] Gage | Command control for many-robot systems[END_REF] presented a first classification of coverage problems, including three varieties: blanket, barrier and sweep coverage. The version considered in our work corresponds to the blanket coverage, defined as finding the static configuration which optimizes a given coverage criterion. A fundamental work on this problem was presented by Cortés et al. in [START_REF] Cortés | Coverage control for mobile sensing networks[END_REF], where an approach based on the centroidal Voronoi partition was presented to drive the agents, moving in two-dimensional convex areas, toward the final configuration. Successively, numerous papers proposed extended versions of this solution, taking into account more complex environments and sensing models, but still making the assumption of perfect initial knowledge on the environment. [START_REF] Ganguli | Visibility-based multi-agent deployment in orthogonal environments[END_REF] studied the case of a coverage criterion based on the visibility of the environment, assumed to be orthogonal, and using omnidirectional sensors with unlimited range. More recent works constantly increased the complexity of the sensing models and start studying the scenario of sensors moving in three dimensions. In [START_REF] Papatheodorou | Collaborative visual area coverage[END_REF], the authors presented a solution for the optimal visual coverage of a planar surface with a team of aerial vehicles equipped with down-facing cameras. The solution is based on a distributed control law maximizing a coverage-quality criterion, where the quality is considered decreasing with the altitude. A similar scenario is considered in [START_REF] Funada | Visual coverage control for teams of quadcopters via control barrier functions[END_REF], but with the objective of not leaving any area unobserved in between the fields of view of the visual sensors mounted on the UAVs, while penalizing overlaps between them. Still considering planar surfaces, [START_REF] Bousias | Collaborative visual area coverage using aerial agents equipped with ptz-cameras under localization uncertainty[END_REF] studies a more complex model for the UAVs sensors, which are here Pan-Tilt-Zoom (PTZ) cameras, leading to a more complex planning problem. Passing to 3D surfaces, the case of visual optimal coverage with a network of omnidirectional sensors deployed on the surface is tackled in [START_REF] Temel | On deployment of wireless sensors on 3d terrains to maximize sensing coverage by utilizing cat swarm optimization with wavelet transform[END_REF], where the authors present a solution based on Cat Swarm Optimization. A similar problem but using UAVs is instead studied in [START_REF] Thanou | Distributed visibility-based coverage using a swarm of UAVs in known 3D-terrains[END_REF], where the solution is again based on the concept of Voronoi tessellation of the environment. Contrary to our case, in this work the environment to cover is assumed to be initially known. [START_REF] Zhang | 3d model-based multi-camera deployment: a recursive convex optimization approach[END_REF] proposed an alternative approach, based on recursive convex optimization, for multi-camera deployment to visually cover a three-dimensional object, which is still considered perfectly known from the beginning.

Most of these solutions, exploiting geometric and gradient based approaches, often require strong assumptions on the environment or on the constraints of the problem. To deal with more challenging scenarios, derivative-free approaches can be suitable alternatives, as shown in [START_REF] Renzaglia | Multi-robot threedimensional coverage of unknown areas[END_REF] where a centralized, stochastic optimization algorithm was adopted to drive a team of UAVs to cover unknown 3D terrains. A modified version of this algorithm, proposing sparse regression techniques, was then presented in [START_REF] Tseng | Near-optimal probabilistic search via submodularity and sparse regression[END_REF] where a target detection problem is formalized as a coverage problem with the objective of finding minimum-time trajectories. The main limitation of this class of solutions remains the possibility of having a strong dependency on the initial conditions and their effect on the final solution. Our contribution on this problem aims to significantly reduce this problem and so overcome this important limit.

Autonomous Exploration

The autonomous exploration is a fundamental task in multi-robot applications and for this reason it has been largely studied and presents a very rich existing literature. [START_REF] Juliá | A comparison of path planning strategies for autonomous exploration and mapping of unknown environments[END_REF] provides an extensive overview and comparison of some of the main solutions. Another interesting survey, more focused on the communication aspects, can then be found in [START_REF] Amigoni | Multirobot exploration of communicationrestricted environments: A survey[END_REF].

Most of the existing strategies for autonomous exploration exploit the frontiers between the already explored area and the still unknown environment to guide the robots during the mission. The first frontier-based solution, for a single robot, was presented by Yamauchi in his seminal paper [START_REF] Yamauchi | A frontier-based approach for autonomous exploration[END_REF], and then extended to the multi-robot case in [START_REF] Yamauchi | Frontier-based exploration using multiple robots[END_REF]. In this work, each robot simply selected its closest frontier point at each iteration until the environment is completely explored. Following these first solutions, numerous variations, especially considering the exploration of two-dimensional areas, have been proposed in the following years. An expensive but efficient centralized approach to solve the problem of simultaneously considering traveling costs and gain utilities was presented [START_REF] Burgard | Coordinated multi-robot exploration[END_REF]. [START_REF] Bautin | Minpos: A novel frontier allocation algorithm for multi-robot exploration[END_REF] proposed a decentralized approach that allocates frontier points based on a rank, in terms of travel distance, to obtain a well balanced spatial distribution of robots in the environment. The problem of frontiers assignment in multi-robot exploration is also studied in [START_REF] Faigl | Goal assignment using distance cost in multi-robot exploration[END_REF], where the authors present a traveling salesman problem formulation considering all current possible goals for each robot. Even though this solution can provide shorter exploration times with respect to alternative solutions, it implies a very high computational cost. Besides the frontier allocation problem, another crucial aspect of multi-robot applications can be to comply with communication requirements. In [START_REF] Banfi | Strategies for coordinated multirobot exploration with recurrent connectivity constraints[END_REF], the authors consider this problem by presenting a strategy that allows agents to coordinate in order to form a network that satisfies recurrent connectivity constraints.

Considering 3D scenarios, a direct extension of a frontier-based method has been presented in [START_REF] Dornhege | A frontier-void-based approach for autonomous exploration in 3d[END_REF]. In three dimensions, even the simple process of identifying the frontiers during the exploration becomes more challenging and costly and has been a subject of research [START_REF] Zhu | A 3d frontier-based exploration tool for mavs[END_REF]. In [START_REF] Tang | An autonomous exploration algorithm using environment-robot interacted traversability analysis[END_REF], the authors propose an analysis of the traversability of the scanned area and an improved K-means algorithms to reduce the dimension of frontiers in the map and improve the efficiency of the results but, contrary to our case, the exploration is carried out by ground robots. A solution, also based on frontiers, to guide a team of UAVs with embedded vision in a three-dimensional environment is then presented in [START_REF] Mahdoui | Cooperative frontier-based exploration strategy for multi-robot system[END_REF]. The novelty of this work is to propose the exchange of only frontier points between robots instead of the entire grid map to reduce communication bandwidth.

In complex three-dimensional environments, constantly relying on frontiers has however some important disadvantages. The computation cost for efficient frontier allocation algorithms becomes very high and selecting the optimal observation point for a given frontier is often not trivial. To overcome these problems, the solution that we propose for this case has the intent to significantly limit the use of frontiers thanks to the combination with the decentralized local optimization algorithm.

A Common Problem Formulation for Coverage and Exploration

The goal of the multi-robot planning problem considered in this paper is to autonomously acquiring information of an unknown surface S in a 3D environment with a team of UAVs. This problem can be formally defined upon the tuple (N, S, P, C) where:

-N is the number of UAVs r ∈ {1, . . . , N }; -S is the unknown set of 3D points (voxels), denoted as s ∈ S, representing the surface of interest; -P is the set of all robots' positions up to the current time T , {P (r) t } r=1,...,N t=1,...,T , with Let us also define the portion of the surface S observed at time t by the robot r as O (r) t (S). A point s is considered as observed by a UAV if: 1) the point and the robot are connected by a line-of-sight and 2) they are at a distance smaller than the sensing range. Since the environment is initially unknown, this information can be uniquely retrieved on-line based on the acquired measurements.

P (r) t . = [x (r) 1,t , x (r) 2,t , x
Optimal static coverage and exploration are two very related problems that can be formulated by defining two slightly different optimization functions, J C and J E respectively, both in terms of this tuple (N, S, P, C). The goal of these optimization problems is then to find the UAVs' trajectories {P (r) t } r=1,...,N t=1,...,T that, by safely navigating through the three-dimensional environment, maximize the acquired information and lead the team to i) either find the static configuration that maximizes the surface visibility in the coverage case or ii) minimize the time in which the entire surface is observed for exploration.

Formally, the coverage problem can be defined by introducing the optimization function J (r) C as follows:

J (r) C (t) = s∈O (r) t 1 |µ t (s)| (1)
where µ t (s) is the subset of UAVs that can observe the point s at time t, i.e.

µ t (s) = {r ∈ {1, . . . , N } | s ∈ O (r) t (S)} . (2)
In other words, the terms |µ t (s)| has the role of penalizing redundant coverage by counting how many robots are observing the same element s at the same time, equally distributing the rewards among them and, as a result, forcing the team to spread out over the environment. Note that the objective function [START_REF] Amigoni | Multirobot exploration of communicationrestricted environments: A survey[END_REF] does not depend on the past states, but only on the robots configuration at a given time t.

In the exploration problem, being the environment completely unknown, the optimal solution which minimizes the mission time cannot be found offline. The exploration needs instead to be guided by suitable heuristics based on the information on the environment acquired on-line during the mission. To this end, a local optimization function similar to (1) can be defined for each robot r to maximize the portion of observed environment at each time step. Before to introduce this function, let us firstly define the surface E t (S) explored up to time t as a function of the observations O (r) k collected at every time k by each robot r:

E t (S) = N r=1 t k=1 O (r) k (S) (3)
and, for each voxel s, the set of robots e(s) that observed it for the first time:

e(s) = {r ∈ {1, . . . , N } | s ∈ O (r) t * s (S), s ∈ E t * s (S)\E t * s -1 (S)} , (4)
with t * s being the time of first observation of the voxel s. We can so define the objective function for the robot r as follows:

J (r) E (t) = s∈Et(S) δ r e(s) |e(s)| where δ r e(s) = 1 if r ∈ e(s) 0 if r / ∈ e(s) (5)
The difference with (1) is that now, by maximizing this function, the robots try at each iteration to maximize the quantity of newly observed points s (δ term), while minimizing the redundant observations (i.e. simultaneous exploration of the same areas by more than one robot). This is due to the term |e t (s)|, analogous to the term |µ t (s)| in [START_REF] Amigoni | Multirobot exploration of communicationrestricted environments: A survey[END_REF]. Note that the sum over all the robots of the functions J

E (t) = N r=1 s∈Et(S) δ r e(s) |e(s)| = |E t (S)| . (6)
Once these objective functions have been defined, our scope is to show that the optimization problems can be approached with the same kind of stochastic optimization algorithm. Additionally, we show that, being the optimization functions a combination of N single-robot functions, the solution can be completely decentralized and each robot can decide its next step only based on the information in its possession.

Local Stochastic Approximation-based Optimization

A first problem to solve is that the optimization functions J C and J E do not have a known analytical expression because of their dependency on the unknown environment to observe and, as a result, they cannot be directly optimized. However, each robot can obtain their values at a given time t, based on the data collected on-line and the information shared with the other teammates. By exploiting these data, these functions can be approximated and the result used to tackle the optimization problem [START_REF] Spall | Stochastic optimization[END_REF]. To this end, we consider a decentralized version of a recent stochastic optimization algorithm, the Cognitive-based Adaptive Optimization (CAO) algorithm. This algorithm was originally developed by Kosmatopoulos in [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF] and has been then proposed in various, mostly centralized, versions to tackle multi-robot planning problems, such as: multi-robot deployment [START_REF] Renzaglia | Multi-robot threedimensional coverage of unknown areas[END_REF], [START_REF] Scaramuzza | Vision-controlled micro flying robots: from system design to autonomous navigation and mapping in GPS-denied environments[END_REF], multi-robot localization and mapping [START_REF] Kapoutsis | Real-time adaptive multi-robot exploration with application to underwater map construction[END_REF] and target detection [START_REF] Tseng | Near-optimal probabilistic search via submodularity and sparse regression[END_REF] .

The algorithm is composed of two main steps: i) an approximation function is computed at each iteration to have a local estimation of the objective function J; ii) the obtained local approximation is maximized, respecting the problem constraints, to generate the next robot position. These two steps are detailed below.

Objective Function Approximation

At each time t, J is approximated according to:

Ĵ(r) t x (r) 1,t , x (r) 2,t , x (r) 3,t = ϑ T t φ x (r) 1,t , x (r) 2,t , x (r) 3,t (7)
where ϑ t is a vector of parameter to estimate and φ the nonlinear vector of regression functions. In this work, we fix these functions in φ to be polynomials of the state variables up to the third-grade, i.e. the vector components φ i are all the possible combinations given by:

x (r) 1 u x (r) 2 v x (r) 3 z s.t. u, v, z ∈ {0, 1, 2, 3}, u + v + z ≤ 3 . ( 8
)
The choice of the basis functions is not unique and other solutions can also be adopted. However, it is worth remarking that this approximation is strictly local and it has no ambition to reproduce the original function in the entire exploration space. Its role is exclusively to provide an estimation of J to allow each robot to select its next position, so it needs to be reliable only in the proximity of the current robot's position. The parameter estimation vector ϑ t , which has dimension L = 20, is calculated at each time t based on a limited set of past measurements according to:

ϑ t = argmin ϑ 1 2 t-1 = t J (r) -ϑ T φ P (r) 2 (9)
where t = max{0, t -L -T h } with T h being a user-defined non-negative integer and J (r) is the value of the objective function for the robot r at time based on the collected data. The problem in ( 9) is an ordinary -unweighted -least square problem linear in the parameters θ i . Its solution can be for instance obtained by considering the Moore-Penrose inverse of the matrix r) , where φ i are the components of the vector φ defined in [START_REF] Dornhege | A frontier-void-based approach for autonomous exploration in 3d[END_REF]. Note that using only a limited set of past values has the advantage to significantly reduce the computation time required to obtain the parameter vector ϑ t but the resulting function is a reliable approximation only locally and cannot be used for global optimization.

X i = φ i P (

New State Selection

Once the approximation Ĵt is built according to ( 7) and ( 9), the next robot position is obtained performing a local random search. A set of admissible new positions, i.e. which respect the problem constraints, is generated by randomly perturbing the current state and tested on Ĵt . More formally, the set of M possible states is constructed as follows:

P (r,m) t = P (r) t + α t ζ (r,m) t , ∀m ∈ {1, . . . , M } , (10) where ζ
(r,m) t is a zero-mean, unit random vector with dimension equal to the dimension of P (r) t and α t represents the exploration step. The standard condition that α t must satisfy to guarantee the convergence to a local optimum (see e.g. [START_REF] Bertsekas | Gradient convergence in gradient methods with errors[END_REF]) is:

lim t→∞ α t = 0, ∞ t=1 α t = ∞, ∞ t=1 α 2 t < ∞ . (11)
It is worth noting that this conditions must be satisfied only for the coverage problem, where the robots need to converge to a local optimum of the problem. On the contrary, in the exploration case, we do not have the necessity of a timedependent optimization step since there is no interest in moving less than the maximum feasible step and there is no convergence until the total environment is explored. The candidates that correspond to non-feasible positions, i.e. that violate the constraints C, are neglected and the new robot position is P (r)

t+1 = P (r,m * ) t
, where m * is:

m * = argmax m∈{1,...,M } Ĵ(r) t P (r,m) t . (12)
The random generation of the candidates is crucial for the efficiency of the algorithm because it ensures that the space is correctly sampled, without any persistent bias. In other words, it allows gathering the necessary information to reconstruct an unbiased approximation function and so guarantees that Ĵt is a reliable estimate of the unknown function J (see [START_REF] Kosmatopoulos | An adaptive optimization scheme with satisfactory transient performance[END_REF], [START_REF] Ioannou | Robust adaptive control[END_REF] for more details). We finally remark that in this work we focus only on high-level planning, without explicitly considering any dynamics constraints. However, more restricting assumptions, including the robot's dynamics, and so reducing the space of feasible positions would affect only the definition of C but not the algorithm.

Overcoming Local Optima

The optimization algorithm presented in the previous section can be employed by each robot to maximize the information gathered on-line on the environment and, depending on the choice of the optimization function, achieve the coverage or exploration task. However, in both cases, the proposed approach has the drawback to be strictly local and be subject to the limitations due to local optima. To overcome this problem, we propose two different strategies for the coverage and exploration tasks, allowing the team to respectively finding a better static configuration or having the guarantees of completing the exploration of the environment.

Static Coverage

In the coverage problem, the multi-robot team has to find the best static configuration which allows a maximum observation over the environment. If, as in our case, this search is exclusively conducted online, without the help of global information, the system can converge to unsatisfactory sub-optimal configurations. In this scenario, a critical choice that significantly influences the final result is the initial deployment from which the robots start their optimization.

In many practical situations, some prior information on the environment to cover, even if very limited and uncertain, can be exploited to compute initial suitable configurations.

In this section, we intend to study the effect of such an initialization phase and, in particular, we propose the use of the constrained Centroidal Voronoi Tessellation [START_REF] Okabe | Spatial tessellations[END_REF], computed over an approximating surface, to achieve this initial step. The optimization process is thus modified by the addition of a new layer to exploit the possible presence of partial prior information. The resulting method, exploiting the complementarity between stochastic and geometric optimization, significantly improves the performance of the previous approach by reducing the probability of a far-to-optimal final solution. Furthermore, the number of iterations needed by the online algorithm to converge, and so the time required to achieve the mission, is also reduced.

To obtain the initial deployment, we assume that few real or predicted points of the surface to cover are available. Based on this information, standard methods for surface reconstruction can be applied to build an approximating surface S1 . To obtain the surface S , we propose to use a piecewise linear interpolation of the available points {q * i } n i=1 , as shown in Fig. 2. The next step is to compute a partition on this approximating surface that allows us to have a balanced distribution of the UAVs over the environment and define their initial positions. To do so, we propose to exploit the centroidal Voronoi tessellation.

Voronoi tessellation is a fundamental concept in Locational Optimization theory [START_REF] Okabe | Spatial tessellations[END_REF]. Given a bounded set Ω ∈ R n and N points {z i } N i=1 , the Voronoi tessellation of Ω generated by

{z i } N i=1 is {V i } N i=1
, where:

V i = {u ∈ Ω : |u -z i | < |u -z j | ∀j = i}, i = 1, . . . , N . (13)
The tessellation is then called Centroidal Voronoi Tessellation (CVT) if z i = z * i , ∀i = 1, . . . , N , where z * i is the centroid of the i-th region, defined as:

z * i = Vi u du Vi du for i = 1, . . . , N , (14)
A commonly adopted solution to converge to a CVT, starting from an arbitrary set of generators, is the Lloyd's algorithm [START_REF] Lloyd | Least squares quantization in PCM[END_REF], [START_REF] Cortés | Coverage control for mobile sensing networks[END_REF]. This algorithm is very efficient for simply planar regions but not particularly suitable for complex scenarios since it requires the computation of the Voronoi partition corresponding to the set of points at every iteration. This step is often too difficult for arbitrary surfaces in 3D and probabilistic algorithms are more efficient. In this paper, we adopt a modified version of the MacQueen's algorithm, firstly proposed in [START_REF] Ju | Probabilistic methods for centroidal Voronoi tessellations and their parallel implementations[END_REF] and then extended to general surfaces in [START_REF] Du | Constrained centroidal Voronoi tessellations for surfaces[END_REF]. Each iteration of this algorithm is composed of the following steps: 1. Fix a positive integer q and the constants {β i , γ i } i=1,2 such that:

β 2 > 0, γ 2 > 0, β 1 + β 2 = γ 1 + γ 2 = 1 ; 2.
Choose a set of N random points {z N i=1 } ∈ S and set j i = 1 for i = 1, . . . , N ; 3. Select randomly q samples {y r } q r=1 ∈ S ; 4. For r = 1, . . . , q define z i * r as the closest {z} N i=1 to y r ; 5. For i = 1, . . . , N define the set W i as the set of all samplings y r closest to z i , compute the average y * i of the set W i and define:

ẑ * i ← (β 1 j i + γ 1 )z i + (β 2 j i + γ 2 )y * i j i + 1 , j i ← j i + 1 ; (15)
a projection operator is then used to impose the surface constraint, i.e. z * i = proj(ẑ * i ), and the set of {z * i } becomes the new set {z i } N i=1 . These steps are repeated until a convergence criterion is fulfilled. An example of a CVT of a non-planar surface is shown in Fig. 3. It is worth noticing that such a configuration does not provide any guarantee in terms of visibility of the environment and so it does not represent an optimal solution to the visual coverage problem. Its role is exclusively to generate a first deployment from which the optimization algorithm will start. An alternative approach could be to find the global optimum solution corresponding to the approximated surface. However, we do not consider this possibility for two reasons: firstly, this initialization step can be also performed online as part of a bigger mission, and the computational cost and time required to solve the global optimization problem would be excessive for the limited resources usually available in these scenarios; secondly, this solution would still be a sub-optimal solution of the original problem because computed over a coarse approximation of the real surface.

In Section 6, the performance of the proposed initialization process will be evaluated in simulation and compared with other possible solutions to show the effectiveness of this additional layer in achieving coverage tasks.

Exploration

In the exploration task, a solution only based on a local optimization approach cannot guarantee to lead the robots to explore the entire environment. In other words, as long as new regions are explored as the result of a robot trajectory, the algorithm can provide a direction to follow, but when a robot remains completely surrounded by previously explored regions, the optimization algorithm only generates random movements and needs external input to identify new unexplored areas. To overcome this problem, we propose to combine the stochastic optimization solution with a frontier-based approach. In practice, whenever the gain of information ceases to increase, the set of remaining frontiers in the environment is considered and the robot is driven towards the closest one. This additional step allows the robot to efficiently overcome the deadlocks and restart the local optimization from the new position. The resulting decentralized algorithm, called the Frontier-based Cognitive Adaptive Optimization (FCAO) algorithm, has the goal to overcome the limitations that these two approaches have separately while maintaining their strengths. In particular, alternating local explorations with frontier assignments leads to guarantee the completeness of the exploration as in any frontier-based method. The structure of the FCAO is reported in Algorithm 1, where the function CAO is called until the local optimization is no more able to improve the objective function J E . In this case, the function Closest Frontier selects a new goal among the list of available frontiers. The exploration is completed when this list is empty.

To build the list of frontiers, each robot makes use of its own global map that is updated by merging the collected data with local versions of the maps received from other robots (see [START_REF] Yamauchi | Frontier-based exploration using multiple robots[END_REF] for a more thorough description of this architecture). The fundamental advantage of this strategy is that, even though some information needs to be shared, the decision is completely decentralized. This makes the system robust with respect to single robots' failures because, if a robot stops transmitting information, the other teammates can continue the mission until the entire environment is explored. Additionally, being the system asynchronous, robots do not need to wait for other robots to take their decisions.

The limitation of this method for frontiers allocation is the possibility of redundancy in the exploration, i.e. the same frontier may be assigned to more than one robot. This drawback is also present in our approach but it is limited by the fact that frontiers are used only rarely and asynchronously by the robots.

Algorithm 1: The FCAO Algorithm

P t+1 ← Pt + αζ (m * ) t ← t + 1 until J E (t) = J E (t -1); function ClosestFrontier(Pt, E) F = CandidateFrontiers(E) return SelectClosestFrontier(F , Pt)
6 Performance Analysis for Coverage

The solution presented in Section 5.1 for the visual coverage has been tested and evaluated in simulation. Our primary objective is to show the importance of an initial positioning as a way to overcome the limits of online local optimization. To this end, after testing this approach with generic surfaces, we provide comparative results adopting a different or none initialization. Other advantages of using the Voronoi-based initialization, such as achieving faster convergence and increasing the robustness with respect to the environment's characteristics and team capabilities, are here also discussed. Afterwards, a study on the dependency on the quantity of initial information on the surface is provided. Finally, a more realistic simulation, in terms of both UAVs model and environment to cover, has been carried out in Gazebo to show the feasibility of this strategy in realistic scenarios.

Experimental Settings

In all the results here reported, the UAVs have sensing capabilities constrained by a limited sensing distance d max . The UAVs' trajectories are defined as sequences of way-points and we assume the presence of an internal controller that allows following them. The environments are assumed to be arbitrary and two distinct scenarios have been taken into account: i) smooth terrains generated as mixtures of Gaussians and ii) an outdoor environment composed of several buildings. These environments are assumed to be unknown except for a few points, which have been randomly selected on the surfaces. In both cases, the coverage objective function ( 1) has been computed by discretizing the environment and considering the number of points belonging to the surface visible from the UAVs. To facilitate the comparison between different scenarios, these objective functions have then been normalized with respect to the total surface to cover, and are so expressed as the percentage of environment covered by the UAVs. The parameters of the CAO algorithm have been set as proposed in [START_REF] Renzaglia | Multi-robot threedimensional coverage of unknown areas[END_REF].

Illustrative Example

Fig. 4 presents an illustrative example where the main steps of our approach are shown. In this scenario, 9 UAVs with a sensing range d max = 35m have the goal of covering an environment for which a limited set of few points is the only initial information available. The real surface is so approximated by a piecewise linear interpolation based on these points. The approximating surface is then used to compute the constrained CVT generated by the UAVs positions and the obtained configuration is the initial state from which the stochastic optimization algorithm starts. From the behavior of the coverage level during this last optimization step, it is possible to see that UAVs can finally monitor almost the whole environment. It is important to note that the initial positions provided by the CVT can be below the real surface due to surface approximation errors. When this happens, the UAVs will converge to the closest admissible location avoiding any collision with the ground.

Comparative Results

In order to show the advantages of a CVT-based initialization, we compare the results corresponding to the coverage obtained: i) without any initialization, i.e. with the UAVs starting in a corner of the environment, and ii) with randomly initialized sets of positions generated with a uniform distribution over the region to cover. As a first result, we analyze the coverage level achieved during the online optimization, after the three different initialization strategies, in an environment similar to the one showed in Fig. 4. From these results, reported in Fig. 5, we can see that the local optimization algorithm without any initialization is not able to find a good solution and the UAVs remain stuck in a bad local optimum. A better solution is then obtained with a random positioning of the resources. However, the CVT-based solution starts from an even higher coverage value than the final levels achieved by the two other solutions. In particular, in this case the initialization is very close to a local optimum and the online algorithm can only slightly improve the initial coverage.

To have more representative results, we carry out the same comparison but varying the number of robots and the characteristics of the environment to cover. For every fixed number of UAVs, we consider fifty different scenarios. The environments are generated as mixtures of seven Gaussians centered in random locations, with varying height and fixed variance. To generate the CVTs, a set of thirty randomly selected points is used to approximate the 1 Average and standard deviation of percentages of covered environment achieved over 50 scenarios, corresponding to terrains randomly generated in a 100x100m 2 square, and considering four different team sizes. The sensing range dmax is 25m. For the CVT and the random initialization, also the initial coverage is reported.

surface. We analyzed the performance corresponding to four different team sizes: 5, 8, 12 and 20 UAVs. All the simulations are bounded with 500 iterations, which ensure the convergence of the optimization algorithm. In Table 1, we reported for each scenario the average initial coverage (except for the non-initialized case) and the average maximum coverage achieved, all with respective standard deviations σ.

These values show that the optimization started from the Voronoi-based configuration always achieves the highest coverage level and with the smallest standard deviation. More specifically, the CVT initialization produces an improvement with respect to the random initialization comprised between 5% and 12%. Comparing then the coverage levels obtained after the initialization step, it is even more clear the advantage with respect to the random initialization. Note that these improvements also translate in a significant reduction of iterations, and so time, required to accomplish the task.

It is also interesting to analyze these performances in terms of team size. We can notice that the improvement obtained with the CVT initialization is more relevant in the intermediate cases, i.e. with 8 and 12 UAVs. This can be motivated by the fact that they represent more challenging scenarios where the optimal configuration is harder to find. In fact, with few UAVs, it is sufficient to achieve a good spread of the robots over the environment to reach configurations with non-intersecting fields of view. In such configurations, the team is close to its maximum sensing capability and so they correspond to close-to-optimal solutions. On the other side, in large teams, a high level of overlapping observations is inevitable and several distinct configurations can lead to the same total coverage of the environment (in the extreme case of a very large number of robots, even a random deployment may lead to full coverage). However, an initial dispersion of the UAVs can still be beneficial to avoid that some of them remain stuck in their initial location reducing the total sensing capability. This can be seen in the case of 20 UAVs, where the coverage without initialization is significantly lower than the other two.

Dependency on Prior Information

We now want to stress the importance of the CVT to uniformly distribute the robots over the environment and provide a good initialization even when the information on the surface to cover is extremely limited. For one scenario reported in Table 1 (12 robots, third line), we perform the same study over the fifty environments but now considering the approximating surface based on just the four points defining the boundaries of the region plus one point in the middle of the environment. The resulting approximating surface is in this case simply a pyramid. The CVT initialization for this extreme case provides an almost uniform distribution of the UAVs over the environment. The coverage level obtained after the initialization with this surface is 0.78 with a standard deviation σ = 0.03, while the final maximum coverage is 0.91 with a standard deviation of σ = 0.02. As expected, the first result shows that the initial value is lower than the one obtained with a more informative approximation, but it remains a better solution than a random initialization. Then, the online optimization is able to drive the UAVs to achieve a final coverage very close to the result of the more informative initialization. This highlights the important effect of the Voronoi tessellation to spread the team and generate a suitable initialization.

Urban Environment

To evaluate our approach in a more complex scenario, we present here the results obtained by carrying out simulations in Gazebo [START_REF]Gazebo robot simulator[END_REF] and realized using the ROS package tum simulator [START_REF]TUM ROS simulator[END_REF]. These simulations include more realistic models for both the environment and the system, showing the employability of our approach in realistic applications. In the proposed scenario, four simulated UAVs with a maximum visual range of 8m are employed to cover a 24mx24m area with several buildings. The steps illustrated in Fig. 6, from top left to bottom right, are: the initialization, computed on an extremely simplified reproduction of the environment; the final positions reached after the on-line optimization performed by the CAO algorithm; the covered environment represented as a point cloud; the coverage level achieved during the mission, here also compared with the result obtained from a different initialization (two UAVs starting close to each other on one side of the environment and the other two similarly on the opposite side). This comparison clearly shows that the initialized case begins from a higher value and converge significantly faster to the final positions, not far from the initialization.

Performance Analysis for Exploration

Considering now exploration tasks, the proposed FCAO algorithm has been tested in simulations and compared with possible alternative solutions. As for the coverage results presented in the previous section, the surfaces to explore are arbitrary and, in a first study, created as random mixtures of Gaussians with different heights and variances. Fig. 7 presents an illustrative example with four UAVs and their respective trajectories. In all the results reported in this section, the robots' initial locations can slightly vary but they are always in a corner of the environment, close one each other.

To evaluate the performance of the proposed FCAO algorithm, it has been compared with two alternative approaches:

1. Closest frontier: a decentralized frontier-based algorithm, where each robot selects its closest frontier point independently from the others. This algorithm corresponds to the standard solution presented in [START_REF] Yamauchi | Frontier-based exploration using multiple robots[END_REF]. 2. Greedy: a centralized solution that 1) divides frontiers in clusters and considers their centroids as candidate frontiers to improve the distribution among robots and reduce redundancy in the exploration (as done for instance in [START_REF] Tang | An autonomous exploration algorithm using environment-robot interacted traversability analysis[END_REF]), and 2) assigns clusters' centroids sequentially to each robot to avoid multiple assignments of the same frontiers to different robots. Note that this greedy assignment algorithm is a variation of [START_REF] Burgard | Coordinated multi-robot exploration[END_REF] considering constant frontiers utilities2 .

Fig. 8 presents a first result with eight UAVs starting from the same set of initial locations. As a first remark, it is possible to see that our solution significantly outperforms the closest-frontier algorithm. This clear effect can be explained with the high level of redundancy in the exploration only considering the closest frontiers selection. On the contrary, the local optimization efficiently spreads out the robots, avoiding important overlap in the observations. Then, it is interesting to notice that, even though the decision is decentralized, our approach obtains a performance very close to the centralized greedy solution.

Considering the same scenario, we then analyzed the computational time required to find the next robot position at each iteration. One of the strong advantages of our method is its low computational cost and it is clearly shown in Fig. 8(Right). This factor is indeed not only significantly below the time needed by the greedy algorithm but it is even faster than the simple closest frontier selection. The main computational difference here is the necessity for frontier-based algorithms to compute each time the distances between each robot and all the frontier points, an operation that becomes very expensive in three dimensions. This effect is visible even in our analysis where we do not compute the exact geodesic distance between two points, but instead we consider only the straight line connecting the two points and adjusting the altitude to pass above the surface, including a safety distance from the terrain. We can also see that the run time of the local optimizer is constant (except for the very few steps where frontiers are exploited) while the others are varying during the time because of the dependence on the number of frontiers in the map at each iteration. Moreover, being decentralized, its cost is also independent of the number of robots, contrary to the centralized greedy algorithm. The only Fig. 8 The FCAO approach is compared with two alternatives: a decentralized closest frontier algorithm and a centralized greedy selection of frontiers. The results correspond to a scenario with 8 robots in an environment similar to the one shown in Fig. 7. Left: Portion of explored environment as a function of time. Right: Time required to identify the next way point at each iteration. relevant cost of the local optimization comes from solving the problem in [START_REF] Du | Constrained centroidal Voronoi tessellations for surfaces[END_REF], which depends only on the dimension of the vector φ and the number of past considered measurements, which are both constant during the exploration.

Finally, we compared these three possible solutions in 100 randomly generated environments (eight peaks with the same variance and centered in randomly selected, uniformly distributed locations). It is worth noticing that, even though these environments are relatively open, the sensing constraints compared with the peaks elevation force the robots to fly among them and not just over them, making the complexity of the environment and the task higher. Furthermore, as for the coverage case, our study is not only limited to very complex environments but also to more generic open outdoor terrains, where the redundancy problem is one of the most relevant to reduce the performance of computationally light solutions. The number of iterations needed to complete the task and the averages for each algorithm are reported in Fig. 9. Similarly to the previous results, the FCAO algorithm performs very closely to the centralized greedy algorithm, and they both significantly outperform the closest-frontier solution. We can also notice that the latter, besides showing a worse average behavior, also presents a greater variance in the results, while the other two are more robust with respect to changes in the environment.

Urban Environment

Similarly to the results of the previous section, we finally tested our approach in more complex and realistic simulations executed using ROS and Gazebo simulator. A rich urban environment including buildings, complex structures, and vehicles has been here considered. The UAV's dynamics have been here defined following the model presented in [START_REF] Martin | The true role of accelerometer feedback in quadrotor control[END_REF]. A controller, based on feedback linearization (for more details see e.g. [START_REF] Jaulin | Mobile robotics[END_REF]), has been implemented allowing the translation of the new positions generated by the optimizer in rotational velocities for the rotors. The components have been developed to be as physical representative as possible, especially regarding the aerodynamics, which has been modeled after the specifics of the Intel Aero quadrotor3 . This allows us to reproduce realistic trajectories for the drones and to show that the way-points are generated in a way that leads to feasible motion.

The UAVs obtain at each position a point cloud representation of the environment included in a down-looking, 20m-radius half-sphere. The map of the observed environment and, if necessary, a list of frontier points are generated based on this point cloud thanks to the ROS package Octomap server [START_REF] Hornung | Octomap: An efficient probabilistic 3d mapping framework based on octrees[END_REF]. The result of the exploration conducted by three UAVs is presented in Fig. 10, where the initial environment, the final 3D map, and the UAVs trajectories are shown 4 . It is possible to see that the solution leads to a complete 3D reconstruction, i.e. without leaving any area unobserved. Furthermore, we can notice that the UAVs, after taking-off from close positions, are driven by the algorithm toward different regions to minimize the information overlap, and so resulting in a reduced exploration time. For a more quantitative analysis, in Fig. 11 we report the number of octree nodes explored over time. The result obtained in Fig. 10 is here compared to the same exploration conducted by a team of five UAVs. From this plot, we can see how the UAVs are able, in both cases, to continuously improve their knowledge of the environment, avoiding getting stuck in already explored regions. As expected, the team composed of five UAVs obtains a better result, requiring about 40% less time to complete the exploration.

Final Discussion

In this paper, we considered the cooperative static coverage and the autonomous exploration problems for three-dimensional environments with a team of aerial robots. The main subject of our study was the possibility to express both problems with a unique framework, showing their similarities with a structure based on the observations gathered online and trying to overcome the usual gap present in their solutions. Formally, these problems maintain a substantial difference in their nature. The solution of the coverage problem is indeed represented by the set of static positions that maximizes the total observation, while in the exploration the entire robots' trajectories are the objective of the problem. However, in our work, we showed that a common structure of single-robot objective functions can be defined and that it allows exploiting for both problems a measurement-based stochastic optimization approach to locally maximize the online acquired information. Furthermore, we showed that this approach, contrary to most of the existing solutions, leads to decentralized and light solutions, particularly suitable for computationally limited platforms.

Additional specific strategies, reflecting the still existing differences of these problems, have been then proposed in order to overcome the well-known limitations of local optimality and obtain more efficient solutions. In particular, for the coverage problem, we presented an algorithm, based on centroidal Voronoi tessellations, to generate suitable initial configurations from which the online optimization leads to significantly better final results in fewer iterations. For the exploration task, we instead proposed the FCAO algorithm which combines the local online optimization with a sparse exploitation of frontier points to obtain a light and efficient approach. In fact, this solution maintains all the advantages of the local stochastic optimization method, avoids the necessity of computing the full distance matrix between robots and frontiers at each iteration, and still guarantees the final total exploration of the environment.

The proposed methodology has been tested in simulation, considering different classes of environments and comparing the obtained results with possible alternative strategies. These results showed the efficiency of our solution for both coverage and exploration. In particular, the FCAO algorithm for cooperative exploration, despite being completely decentralized and less computationally requiring, provided results comparable with a more complex centralized approach and clearly outperformed a standard decentralized alternative.

It is worth mentioning that another significant advantage of a unified formulation for these two tasks is that it eases the process of finding new advanced solutions for problems that lie within the proposed framework, i.e. new combinations between exploration and coverage problems. Additionally, the use of a common framework, based on similar information and on the same optimization approach, could also significantly ease the execution of the two solutions in scenarios where switching between two tasks proves necessary. While the present work focused on independently solving the two problems, we plan to investigate these possible directions as part of future studies.

We also believe that both problems present further interesting research lines to improve our specific solutions. For the coverage problem, we intend to pursue the study on more informative coverage criteria, including also the quality of acquired visual data in the objective function. This will affect the generation of the initial partition since new metrics could be considered instead of standard distances for the Voronoi tessellation. For the exploration problem, we consider that relevant questions remain open on frontier-based approaches for 3D environments, such as finding fast and efficient strategies to define the best viewpoints once a frontier is identified. studying the effects of uncertainties generated by localization errors and the robustness of the proposed strategies will also be an important point to analyze.

  ] being the position of the UAV r at time step t; -C(P (r)t )0 is the set of constraints that each UAV position P (r) t must satisfy, e.g. robot-to-obstacle and robot-to-robot collision avoidance, minimum and maximum height of flight.

E

  (t) is simply the total explored environment up to time t, i.e.,

Fig. 2

 2 Fig. 2 Surface approximation: the real surface (red points) is approximated by a piecewise linear approximation based on 50 randomly selected points (black points).

Fig. 3

 3 Fig. 3 Constrained CVT with six generator points.

Fig. 4

 4 Fig. 4 Coverage obtained with 9 robots having a sensing range dmax = 35m. From top left to bottom right: 1. The real surface (in red) is approximated by a piecewise linear interpolation of few points (black dots). 2. CVT on the approximating surface. 3. Initial and final positions are in black and green respectively. The surface in red is covered by at least one UAV. 4. The portion of covered surface as a function of CAO iterations.

Fig. 5

 5 Fig. 5 Covered environment during the online optimization with a team of 8 UAVs. In red the UAVs start close one each other without initialization step, in blue their positions are initialized randomly with a uniform distribution and in green the CVT-based initialization is used as the first step of the coverage approach. The maximum coverage levels obtained by the three strategies are: 0.55, 0.75 and 0.84 respectively.

Fig. 6

 6 Fig. 6 Four UAVs cover a simulated environment with several buildings. From top left to bottom right: 1. The initial positions are computed exploiting an approximation of the environment. 2. These positions are optimized based on collected data to obtain the final configuration. 3. The map of the environment obtained by the team. 4. The coverage levels obtained after the CVT initialization and a non-initialized configuration.

Fig. 7

 7 Fig. 7 Four UAVs exploring an unknown terrain.

Fig. 9

 9 Fig. 9 Number of steps necessary with the three algorithms to complete the exploration for each of the 100 different randomly generated environments. Dashed lined represent the means for each approach.

Fig. 10

 10 Fig. 10 Environment to explore (top) and final 3D map (bottom) obtained with three UAVs guided by the FCAO algorithm and their respective trajectories.

Fig. 11

 11 Fig. 11 Number of octree nodes explored over time for the scenario presented in Fig. 10. The result with 3 UAVs is in blue, with 5 in orange.

Univ. Grenoble Alpes, Inria, 38000, Grenoble, France •

INSA Lyon, CITI Lab, Inria, Chroma, Lyon, France (e-mail: alessandro.renzaglia@inria.fr, jilles.dibangoye@inria.fr, vincent.le-doze@inria.fr, olivier.simonin@inria.fr)

As we show in Section 6, even without an informative prior it is still possible to generate an initialization able to improve the final result.

The utilities update step as in[START_REF] Burgard | Coordinated multi-robot exploration[END_REF] is already computationally expensive in two dimensions and for

3D environments would make the computational burden of this algorithm too high to be considered in this paper, where we focus on light solutions.

Our simulator is available online at https://gitlab.inria.fr/chroma/ drone-simulator.

A video of this simulation is also available at https://team.inria.fr/chroma/files/ 2020/05/FCAO_2020.mp4