Louis Massucci

Fabien Lauer

Marion Gilson

Structural Risk Minimization for Switched System Identification

This paper deals with the identification of hybrid dynamical systems that switch arbitrarily between modes. In particular, we focus on the critical issue of estimating the number of modes. A novel method inspired by model selection techniques in statistical learning is proposed. Specifically, the method implements the structural risk minimization principle, which relies on the minimization of an upper bound on the expected prediction error of the model. This so-called generalization error bound is first derived for static switched systems using Rademacher complexities. Then, it is extended to handle non independent observations from a single trajectory of a dynamical system. Finally, it is further tailored to the needs of model selection via a uniformization step. An illustrative example of the behavior of the method and its ability to recover the true number of modes is presented.

I. INTRODUCTION

Hybrid systems are dynamical systems that include both interacting continuous and discrete dynamical behaviors. This results in systems that switch, according to the value of the discrete variables, between different operating modes with continuous dynamics. This paper focuses on switched linear systems, for which the continuous dynamics are linear and the switching mechanism is arbitrary, i.e., the discrete variable arbitrarily triggering the switch from a mode to an other is an unobserved external input. The identification of such systems, i.e., their estimation from input-output data, is a complex problem whenever the active mode associated to every data point is unknown [START_REF] Lauer | On the complexity of switching linear regression[END_REF]. A number of efficient methods have been devised over the last 15 years, but a critical issue remains: the estimation of the number of modes. Indeed, as reviewed in [START_REF] Lauer | Hybrid system identification: Theory and algorithms for learning switching models[END_REF], switched system identification methods can be classified into two groups: those that work with a fixed number of modes (such as [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF], [START_REF] Lauer | A continuous optimization framework for hybrid system identification[END_REF], [START_REF] Lauer | Estimating the probability of success of a simple algorithm for switched linear regression[END_REF]), and those that estimate the minimal number of modes satisfying a predefined threshold on the error (such as [START_REF] Bemporad | A bounded-error approach to piecewise affine system identification[END_REF], [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF], [START_REF] Ozay | A sparsification approach to set membership identification of switched affine systems[END_REF] or [START_REF] Ohlsson | Identification of switched linear regression models using sum-of-norms regularization[END_REF] which uses a penalized minimization form of this principle rather than a strict constraint). Therefore, in all cases, the methods require the value of a hyperparameter that directly implies the number of modes for the model.

In this paper, we propose to tackle the model selection problem of estimating the number of modes with a structural risk minimization (SRM) approach inspired from statistical learning [START_REF] Vapnik | Statistical Learning Theory[END_REF]. This approach relies on the minimization of an upper bound on the expected error of the model, called the risk or generalization error. Much of learning theory is LM is with the Université de Lorraine, CNRS, LORIA, CRAN, F-54000 Nancy, France louis.massucci@univ-lorraine.fr FL is with the Université de Lorraine, CNRS, LORIA, F-54000 Nancy, France fabien.lauer@loria.fr MG is with the Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France marion.gilson@univ-lorraine.fr devoted to the derivation of such generalization error bounds, however most often with an independence assumption on the data, which is not compatible with the system identification context where data points come from a single (or a few) system trajectory. Alternatively, most works related to error bounds for dependent data are based on measures of dependence called mixing coefficients [START_REF] Bradley | Basic properties of strong mixing conditions. a survey and some open questions[END_REF] and the independent block sequence construction due to [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF]. These ideas were exploited to produce non-asymptotic error bounds for system identification in [START_REF] Weyer | Finite sample properties of system identification of arx models under mixing conditions[END_REF], [START_REF] Vidyasagar | A learning theory approach to system identification and stochastic adaptive control[END_REF]. In line with more recent work in learning theory, [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF] derived error bounds for dependent data based on Rademacher complexities.

In this paper, we combine the general framework of [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF] with recent results on Rademacher complexities of static switched models [START_REF] Lauer | Error bounds for piecewise smooth and switching regression[END_REF] to derive the first generalization error bound for switched system identification. Then, we derive an SRM approach to estimate the number of modes on the basis of this bound. This requires to produce another uniform error bound suitable for the practical setting of model selection for most switched system identification algorithms. To the best of our knowledge, this is the first time that such statistical learning techniques are applied in this switched system identification context.

Paper organization: Section II introduces the hybrid system framework and formally exposes the identification problem. Section III gives the necessary background on learning theory and derives the error bound for switched systems. Then, Section IV details the proposed model selection approach, which is illustrated on a numerical example in Section V. Section VI concludes the paper and discusses open issues.

Notation: For any positive integer n, [n] denotes the set of integers from 1 to n. An upright bold letter with a subscript, e.g., t n , denotes a sequence (t i ) 1≤i≤n of length n, which should not be confused with the vector t i of index i. Capital letters are used for random variables. Given two sets X and Y, Y X denotes the set of functions from X to Y.

II. GENERAL FRAMEWORK

In this paper, we focus on discrete-time Single Input Single Output (SISO) Autoregressive with external input (ARX) hybrid systems of the form

y i = x T i θ qi + e i (1)
with y i ∈ R the output, x i ∈ R d the regression vector, q i ∈ {1, . . . , C} the discrete state or mode, C the number of modes, θ j with j = 1, . . . , C the parameter vector of the jth mode and e i ∈ R a noise term. The regressor x i ∈ R d , d = n a + n b , with the model orders n a and n b , is given by:

x i = [-y i-1 , . . . , -y i-na , u i-1 , . . . , u i-n b ] T , (2)
where the u i-k 's denote the delayed inputs. The goal of switched system identification is to find a model f = {f j } C j=1 made of C component submodels f j that estimate the continuous behaviors of the system in the different modes. The problem can be set as follows.

Problem 1. Given a data set D = {(x i , y i )} n i=1 and a set of possible submodels F, estimate the number of submodels C, the submodels

f = {f j } C j=1 , f ∈ F C , with F C = F × F × • • • × F and the switching sequence q = (q i ) 1≤i≤n ∈ [C] n .
In this paper, we concentrate on the model selection subproblem of determining the number C of submodels. More precisely, we consider the general scheme in Algorithm 1, in which we assume that we have access to a generic switched system identification algorithm working with a fixed C to estimate f ∈ F C and q ∈ [C] n . In this scheme, the generic algorithm is applied for all number of modes C within a predefined range. Then, the "best" model is selected on the basis of a criterion J(C). 1 Here, we develop a model selection method based on a criterion J(C) inspired by the SRM principle in statistical learning. This criterion, derived below, will basically take the form of a generalization error bound.

III. ERROR BOUNDS FOR SWITCHED SYSTEM IDENTIFICATION

A short introduction to learning theory and regression error bounds is presented in Sect. III-A, before the exposition of error bounds dedicated to static switched systems in Sect. III-B. Then, we discuss the non-independent case in Sect. III-C and finally derive a bound applicable to switched dynamical systems in Sect. III-D.

A. Preliminaries

Let X be the input space and Y the output space with Y = [-M, M ] for some M > 0. A relationship between inputs and outputs is characterized by an unknown probability distribution of a random pair Z = (X, Y ) ∈ X × Y = Z, of probability density function p(x, y). Given a realization of a sample Z n = (Z i ) 1≤i≤n = ((X i , Y i )) 1≤i≤n of n independent copies of Z, the aim of regression is to learn the model f : X → Y that minimizes, over a certain model class F, the generalization error (or risk)

L(f ) = E X,Y (f, X, Y ), (3)
defined as the expected value

E X,Y (f, X, Y ) = X ×Y (f, x, y)p(x, y) dxdy (4)
of the loss function (f, X, Y ). This loss function measures the pointwise error in predicting f (x) instead of y and is typically of the form

(f, X, Y ) = |y -f (x)| p with p ∈ {1, 2}.
As the probability distribution of Z is unknown, we cannot compute the expected value and many methods minimize instead the empirical risk

L(f ) = 1 n n i=1 (f, X i , Y i ), (5)
for a realization of Z n = z n . So-called generalization error bounds are upper bounds on the risk (3) that typically involve the empirical risk [START_REF] Lauer | Estimating the probability of success of a simple algorithm for switched linear regression[END_REF] and hold in the non-asymptotic case with high probability and uniformly over the model class F, i.e., bounds of the form:

P ∀f ∈ F, L(f ) ≤ L(f ) + (n, F, δ) ≥ 1 -δ. (6)
A large part of learning theory is devoted to the characterization of the tightest confidence interval (n, F, δ) in ( 6), which typically depends on the capacity of the model class F. More precisely, the capacity of the loss class

L = { (z) = (f, x, y), f ∈ F} (7)
must be considered, and can be measured for instance with the Rademacher complexity, as initiated by [START_REF] Koltchinskii | Empirical margin distributions and bounding the generalization error of combined classifiers[END_REF], [START_REF] Bartlett | Rademacher and Gaussian complexities: Risk bounds and structural results[END_REF].

Definition 1 (Rademacher complexity). Given a sequence

Z n = (Z i ) 1≤i≤n of random variables Z i ∈ Z, the empirical Rademacher complexity of a class L of functions from Z to R is defined as RZn (L) = E σn sup ∈L 1 n n i=1 σ i (Z i ) Z n , (8)
where σ n = (σ i ) 1≤i≤n is a sequence of Rademacher variables, i.e., random variables uniformly distributed in

{-1, +1}. The Rademacher complexity of L is R n (L) = E Zn RZn (L). (9)
A general bound based on the Rademacher complexity is the following: Theorem 1 (Theorem 1 in [START_REF] Mohri | Foundations of Machine Learning[END_REF]). Let L be a class of functions from Z into [0, B] and Z n = (Z i ) 1≤i≤n be a sequence of independent copies of the random variable Z ∈ Z. Then, for any fixed δ ∈ (0, 1), with probability at least 1 -δ, uniformly over all ∈ L,

E Z (Z) ≤ 1 n n i=1 (Z i ) + 2 RZn (L) + 3B log 2 δ 2n . (10)
To derive risk bounds from Theorem 1, the loss function must be bounded. Given a bounded output Y ∈ [-M, M ], this can be obtained by clipping the model. Definition 2 (Clipping). For any M > 0 and t ∈ R, we define the clipped version t of t as

t =    -M, if t < -M t, if t ∈ [-M, M ] M, if t > M. (11)
The clipped version f of a function f : X → R is obtained by clipping its output: ∀x ∈ X , f (x) = f (x). And F denotes the clipped function class {f : f ∈ F}.

Indeed, since for any (x, y) ∈ X × Y, (f, x, y) ≥ (f , x, y), the risk (3) of the clipped model, L(f ), is always smaller than L(f ) and we can consider that the final result of the estimation is f instead of f and derive upper bounds on L(f ).

Example 1 (Linear regression). For linear regression in X = R d , we can consider the model class

F = {f : f (x) = w T x, w ∈ R d , w ≤ R w } (12)
and the loss class (induced by the clipped model)

L = { ∈ [0, 4M 2 ] Z : (z) = (y -f (x)) 2 , f ∈
Using a contraction argument (see [START_REF] Mohri | Foundations of Machine Learning[END_REF]), the Rademacher complexity of L can be related to that of F, with RZn (L) ≤ 4M RXn (F), which in turn can be bounded using standard computations for Rademacher complexities [START_REF] Bartlett | Rademacher and Gaussian complexities: Risk bounds and structural results[END_REF] as

RXn (F) ≤ R w n i=1 X i 2 n . (13)
Thus, Theorem 1 leads to the following error bound for linear regression: for any δ ∈ (0, 1), with probability at least 1 -δ, for all f ∈ F as in [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF],

L(f ) ≤ L(f ) + 8M R w n i=1 X i 2 n + 12M 2 log 2 δ 2n . (14)

B. Switching Regression Bounds

When estimating arbitrarily switching systems (in the static context of switching regression or for hybrid dynamical system identification), the selection of the active submodel among the C submodels f j is embedded in the loss function:

(f , x, y) = min j∈[C] |y -f j (x)| p (15)
with p ∈ {1, 2}. By minimizing the error thus defined, we ask that at least one of the submodels f j ∈ F accurately approximates the output y.

Error bounds for switching regression in the static case were obtained in [START_REF] Lauer | Error bounds for piecewise smooth and switching regression[END_REF] through the decomposition of the Rademacher complexity of the loss class based on [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF] in terms of the one of F:

RZn (L) ≤ p(2M ) p-1 C RXn (F). ( 16
)
Example 2 (Switching linear regression). For switching linear regression with X ⊆ R d , F as in [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF] and the loss [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF] with p = 2, Theorem 1 and ( 16) combined with [START_REF] Weyer | Finite sample properties of system identification of arx models under mixing conditions[END_REF] guarantee that, for any δ ∈ (0, 1), with probability at least

1 -δ, for all f ∈ F C , L(f ) ≤ L(f ) + 8M CR w n i=1 X i 2 n + 12M 2 log 2 δ 2n .
Note that, due to the independence assumption in Theorem 1, this bound cannot be used directly for switched system identification, where x i depends on lagged outputs y i-k due to (2).

C. Bounds for Dependent Data

The following assumes that the sample Z n is taken from a stationary β-mixing process.

Definition 3 (Stationarity). A sequence of random variables

Z = {Z t } ∞
t=-∞ is said to be stationary if for any t and non-negative integers m and k, the random vectors (Z t , . . . , Z t+m ) and (Z t+k , . . . , Z t+k+m ) have the same distribution.

The time index t does not affect the distribution of a variable Z t in a stationary sequence. Definition 4 (β-mixing). Let Z = {Z t } ∞ t=-∞ be a stationary sequence of random variables. For any i, j ∈ Z ∪ {-∞, ∞}, let σ j i denote the σ-algebra generated by the random variables Z k , i ≤ k ≤ j. Then, for any positive integer k, the β-mixing coefficient of the stochastic process Z is defined as

β(k) = E B∈σ 0 -∞ sup A∈σ ∞ k |P[A|B] -P[A]| . ( 17
) If β(k) → 0 as k → ∞, then Z is said to be β-mixing.
For a sequence of independent variables, P[A|B] = P[A] in [START_REF] Koltchinskii | Empirical margin distributions and bounding the generalization error of combined classifiers[END_REF] and β(k) = 0 for all k ≥ 1. For β-mixing processes, the dependence between two events separated by k time steps weakens as a function of k. For more details on mixing processes, refer to [START_REF] Bradley | Basic properties of strong mixing conditions. a survey and some open questions[END_REF]. In this paper, we simply rely on Definition 4 to characterize the degree of dependence between data.

Rademacher complexity bounds for dependent data were obtained in [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF] by applying the independent block sequence construction of [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF]. In short, the idea is to see the whole data set as blocks of length a and by considering only a subset of µ blocks instead of the set of n data points (see Figure 1). By doing so, for β-mixing processes, concentration inequalities can be applied to a sequence of independent blocks distributed as those in the odd subset while controlling the error thus introduced. This leads to the following adaptation of Theorem 1.

Theorem 2 (Theorem 2 in [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF]). Let L be a class of functions from Z into [0, B] and Z n = (Z i ) 1≤i≤n be a sequence drawn from a stationary β-mixing distribution. For any µ, a > 0 with 2µa = n and δ > 4(µ -1)β(a), with probability at least 1 -δ, uniformly over all ∈ L,

E (Z 1 ) ≤ 1 n n i=1 (Z i ) + 2 RZµ (L) + 3B log 4 δ 2µ , (18)
where δ = δ -4(µ -1)β(a) and Z µ = (Z 2a(i-1)+1 ) 1≤i≤µ is a sample of length µ as in Fig. 1.

Compared with Theorem 1, Theorem 2 shows a confidence interval that decreases with the "effective number of data" µ = n/2a instead of the original number n. In addition, the confidence interval also slightly increases due to the use of δ instead of δ.

D. Bounds for Switched System Identification

We now have all the ingredients needed to derive a new generalization error bound for switched system identification. Theorem 3. Let F C be a vector valued function class with C components from F as in [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF] and the loss be as in [START_REF] Mohri | Rademacher complexity bounds for non-i.i.d. processes[END_REF]. Then, for any sample Z n = ((X i , Y i )) 1≤i≤n ∈ (R d × Y) n drawn from a stationary β-mixing distribution, and for any µ, a > 0 with 2µa = n and δ > 4(µ -1)β(a), with probability at least 1 -δ, the following inequality holds for all f ∈ F C :

L(f ) ≤ L(f ) + 2p(2M ) p-1 CR w µ i=1 X 2a(i-1)+1 2 µ + 3(2M ) p log 4 δ 2µ
,

where δ = δ -4(µ -1)β(a).

Proof: Apply Theorem 2 and bound the Rademacher complexity of the loss class based on ( 15) with ( 16) and [START_REF] Weyer | Finite sample properties of system identification of arx models under mixing conditions[END_REF].

IV. MODEL SELECTION

We now turn to the model selection issue of estimating the number of modes, for which we derive an SRM approach. The SRM principle consists in tuning the hyperparameters by minimizing a generalization error bound, which in our case is given by Theorem 3. However, before we can proceed with the practical implementation in Sect. IV-B, this bound will be tailored to the specific needs of model selection in Sect. IV-A.

A. Uniform Bounds

Since the general scheme in Algorithm 1 requires to compute the bound for all C ∈ [C] with the same data set, we need a bound that holds uniformly over all C. In addition, the bound in Theorem 3 holds with a predefined radius R w for the model class [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF], whereas most practical algorithms for switched system identification do not impose constraints on the parameter vectors w j . To fill this gap between theory and practice, we would need a bound in which the radius R w could be computed a posteriori as

R w = max j∈[C] w j (19)
from the estimated w j 's. This can usually be obtained via a uniformization over a well chosen discretization, as, e.g., in Appendix F of [START_REF] Lauer | Error bounds for piecewise smooth and switching regression[END_REF]. However, here, additional terms in Theorem 3 due to the dependence of the data breaks the possibility to apply the same technique and we have to rely on a predefined grid of K discretized values,

G = {R 1 , . . . , R K }, with R K = R max .
Then, we can derive a bound in which the radius R w is replaced by its discretized value

Rw = min k∈[K] R k , s.t R k ≥ R w , with R w as in (19).
(20)

Theorem 4 (Uniform Bound). Given a constant R max , let F be as in [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF] with R max instead of R w . Then for a sample Z n of size n drawn from a stationary β-mixing distribution, and for any µ, a > 0 with 2µa = n, and δ > 4CK(µ -1)β(a), with probability at least 1 -δ:

∀C ∈ [C], ∀f ∈ F C , L(f ) ≤ L(f ) + 2p(2M ) p-1 C Rw µ i=1 X 2a(i-1)+1 2 µ + 3(2M ) p log(CK) + log 4 δ 2µ , (21)
with δ = δ -4CK(µ -1)β(a) and Rw as in (20).

Note that the cost of the uniformization is merely a log(CK) term within the second square root.

Proof: Let F k be defined as in [START_REF] Yu | Rates of convergence for empirical processes of stationary mixing sequences[END_REF] with R k instead of R w . For any fixed C, R k and δ 0 > 0, let δ 0 = δ 0 + 4(µ -1)β(a) and

(C, R k , δ 0 ) = 2p(2M ) p-1 CR k µ i=1 X 2a(i-1)+1 2 µ + 3(2M ) p log 4 δ 0 2µ .
Then, Theorem 3 gives

P ∃f ∈ F C k , L(f ) > L(f ) + (C, R k , δ 0 ) ≤ δ 0 .
Thus, by the union bound and CK applications of Theorem 3 with confidence δ 0 , we have

P ∃C ∈ [C], k ∈ [K], f ∈ F C k , L(f ) ≥ L(f ) + (C, R k , δ 0 ) ≤ C C=1 K k=1 P ∃f ∈ F C k , L(f ) ≥ L(f ) + (C, R k , δ 0 ) ≤ CKδ 0 .
Then, for any δ = δ + 4CK(µ -1)β(a) with δ > 0, set δ 0 = δ /CK. This leads to δ 0 = δ/CK and thus

P ∀C ∈ [C],k ∈ [K], f ∈ F C k , L(f ) ≤ L(f ) + (C, R k , δ 0 ) = 1 -P ∃C ∈ [C], k ∈ [K], f ∈ F C k , L(f ) > L(f ) + (C, R k , δ 0 ) ≥ 1 -CKδ 0 = 1 -δ (22) with (C, R k , δ 0 ) = 2p(2M ) p-1 CR k µ i=1 X 2a(i-1)+1 2 µ + 3(2M ) p log(CK) + log 4 δ 2µ
.

B. Proposed Method

The proposed model selection method for switched system identification is now detailed. In Sect. IV-A, we derived a uniform generalization error bound that holds uniformly over any number of modes C ≤ C. Thus, a criterion J(C) can be devised on the basis of this bound to work in the general model selection framework of Algorithm 1. More precisely, the structural risk minimization principle amounts in this case to selecting the "best" number of modes Ĉ as the one that minimizes the bound (21) on the generalization error. By leaving aside constant terms that do not depend on C nor on the estimated model f in the bound, the model selection procedure can be written as Ĉ = argmin

C∈[C] J(C), (23)
where

J(C) = L(f ) + (f , C), (24)
(f , C) = 2p(2M ) p-1 C Rw µ i=1
x 2a(i-1)+1 2 µ with a and µ such that 2µa = n. Note in particular that the precise value of β(a) in the bound of Theorem 3 does not influence the model selection procedure and the values of J(C). In this section, we illustrate the method proposed in Sect. IV on an example taken from [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF]: the identification of a switched system composed of C = 3 linear subsystems of orders n a = n b = 2 with parameter vectors

V. NUMERICAL EXAMPLE

Ln (f )

θ 1 = [-0.4 0.25 -0.15 0.08] T , (25)
θ 2 = [1.55 -0.58 -2.1 0.96] T , θ 3 = [1 -0.24 -0.65 0.30] T .
This system is used to generate a data set of n = 400 000 points with (1)-( 2) under the following conditions. The excitation input u i is a zero-mean Gaussian signal of unit variance. The noise e i is a white Gaussian noise whose magnitude is such that the Signal to Noise Ratio (SNR) is equal to 10 dB with respect to the output signal. The active mode q i is uniformly distributed in {1, 2, 3}. Algorithm 1 is applied with J(C) computed as in (24) for p = 1, a = 2 and C = 10. The K-LinReg algorithm [START_REF] Lauer | Estimating the probability of success of a simple algorithm for switched linear regression[END_REF], which remains particularly efficient on large data sets with a satisfactory accuracy, is used to estimate a model f with C modes at each iteration. Figure 2 illustrates the behavior of Algorithm 1 and J(C) together with its two constituents: the empirical error Ln (f ) and the term (f , C) inherited from the bound (21). On the one hand, as submodels are added to f and the number C of modes is increased, the error naturally decreases (top plot). On the other hand, the model complexity increases with C, which implies an increase of (f , C) (middle plot). Finally, the minimum of J(C) at Ĉ = 3 (bottom plot) offers the optimal trade-off between the two terms and leads to the correct estimation of the number of modes.

To evaluate the robustness of the method over several data sets, another set of experiments is conducted. One hundred new data sets are generated under the same conditions as before, but with colored noise ẽi such that ẽi -0.35ẽ i-1 = e i + 0.5e i-1 and a signal-to-noise ratio of 5dB. Figure 3 shows the histogram of the estimated number of modes Ĉ. As can be seen from these results, the method correctly estimates the number of modes (with Ĉ = 3) in 65% of the trials and remains close to the true value otherwise, even in this situation with a large and colored noise.

VI. CONCLUSIONS

This paper investigated a model selection method inspired by the structural risk minimization principle from statistical learning to estimate the number of modes in switched system identification. For this method, new generalization error bounds for switched dynamical systems were derived.

Preliminary results showed that the proposed approach can recover the true number of modes in the considered experimental setting. However, the error bounds on which our model selection method is based inherently depend on the number of data. Currently, the method is effective only when this number is sufficiently large in order to make the confidence interval of the bound of the same order of magnitude as the empirical risk L(f ) and correctly balance the two terms. This is in part due to the dependence in the data, which degrades the confidence interval through the use of the number of blocks µ = n/2a instead of the number of data n. Tightening the bound thus appears as an important issue to make the method applicable to a broader range of problems, especially for the case p = 2.

Future work will also concentrate on the characterization of the β-mixing coefficient for switched systems, for instance by extending results available for linear systems [START_REF] Weyer | Finite sample properties of system identification of arx models under mixing conditions[END_REF], [START_REF] Vidyasagar | A learning theory approach to system identification and stochastic adaptive control[END_REF]. Though the value of this coefficient does not directly influence the proposed model selection strategy, its computation would be needed in practice to apply our generalization bounds as guarantees on the prediction error.

Algorithm 1

 1 General model selection scheme Require: The data set D = {(x i , y i )} n i=1 and a maximum number of modes C for C = 1 to C do Run the generic algorithm to estimate a model f with C modes Compute a criterion J(C) end for Select the "best" number of modes Ĉ = argmin C∈[C] J(C) return the selected model with Ĉ modes

Fig. 1 .

 1 Fig. 1. Illustration of the block sequence with n = 20, a = 2 data points per block and µ = 5 odd blocks (white) and µ even blocks (grey). The sample Zµ = (Z 1 , Z 2a+1 , Z 4a+1 , . . . , Z 8a+1 ) contains the first point of each odd block. The dependence between the odd blocks (and the points of Zµ) decreases with the length a of the even blocks separating them.

For

  any C ∈ [C], we can rewrite F C with radius R max as F C = k∈[K] F C k ; and any f ∈ F C also belongs to F C k with, by (20), k such that R k = Rw . Thus, (22) is equivalent to (21) and the statement is proved.

Fig. 2 .

 2 Fig. 2. Top: the empirical error Ln(f ) decreases as a function of the number C of modes. Middle: The second term in J(C) (24), i.e., (f , C), increases with C. Bottom: the sum of the two terms, i.e., J(C), has a minimum at Ĉ = 3 which coincides with the true number of modes.

Fig. 3 .

 3 Fig. 3. Histogram of the estimated number of modes Ĉ over 100 trials with true C = 3.

J(C) also depends on the estimated model f ∈ F C , but we keep the notation short to put the emphasis on the number of modes C.


Dorin Rautu
email: dorin.rautu@enseeiht.fr

Riadh Dhaou
email: dhaou@enseeiht.fr

Emmanuel Chaput

Initial Placement Optimization for Multi-channel UAV Networks

Keywords: Unmanned Aerial Vehicle (UAV), Drone, Deployment, Channel switching, Ad-hoc network

Unmanned Aerial Vehicles (UAV) can be used to deploy communication networks by acting as access points for ground users. Taking advantage of the lightness and the high maneuverability of drones, such a network can be implemented quickly and inexpensively in situations where network infrastructures are damaged or overloaded (emergency situations), or nonexistent (wild life observation). To mitigate these issues, an off-loading network based on UAVs carrying radio access points was proposed in our previous work. The goal is to temporarily provide multiple services, voice, video, data, etc., over a specific zone.

The design of the aerial network was formulated as a self-deployment method built on a Coulomb's law analogy where users and UAVs act as electrical charges. In this paper, we go beyond the proposed scheme by considering a multi-channel model taking into account the interference. We set up association and channel switching schemes that boost the overall performance of the network.

Introduction

The use of flying aircraft such as unmanned aerial vehicles (UAVs), also known under the name of drones, is fast growing in a wide range of networking applications. In particular, with their ingrained attributes such as mobility, adaptive altitude and flexibility, UAVs concede several key applications in wireless systems.

Benefiting of these characteristics, future generation networks can integrate UAVs to improve the quality of service provided as well as implementing novel functionalities. Drones can be used to extend and support terrestrial networks in information dissemination, [START_REF] Zeng | Wireless communications with unmanned aerial vehicles: opportunities and challenges[END_REF]. D2D (Device-to-Device) networks offer an efficient solution to alleviate terrestrial networks by offloading some traffic, but their benefits are limited as the communications are short-ranged. Drones can play a major role, as they can offer a rapid dissemination platform. As suggested in [START_REF] Orsino | Effects of heterogeneous mobility on D2D-and drone-assisted mission-critical MTC in 5G[END_REF], UAVs can take a major part in vehicular networks as they can facilitate the information spreading by reducing the number of links needed at ground level.

With the rapid expansion of the IoT (Internet of Things) market, the network operators have to rethink how conventional networks operate, to incorporate the massive number of IoT devices (smart-city sensors, health-care wearables, smartphones, vehicles). In poorly covered areas, IoT devices may have a hard time sending messages, as the device energy constraints limit the transmission power. In [START_REF] Pang | Efficient data collection for wireless rechargeable sensor clusters in harsh terrains using UAVs[END_REF], UAVs could be deployed to act as base stations and to provide an energy efficient link for these kind of communications. Thanks to the line-of-sight communications and variable altitude, the signal attenuation can be reduced and the coverage area can be increased.

Natural catastrophes, apart of the devastating material destruction, can bring forth massive communication disruptions as the terrestrial communication networks can be damaged or destroyed. In such events, reliable public safety communications are needed to facilitate first responders deployment, victim search and rescue operations. The use of drone based aerial networks can be a promising solution as for the fast deployment, high coverage and flexibility [START_REF] Merwaday | UAV assisted heterogeneous networks for public safety communications[END_REF].

All of these applications, heterogeneous as may seem, have an important common issue : the drone placement. In order to work efficiently, the UAVs have to be at the right spot at the right time.

We believe that a versatile solution needs to be dynamic (as the users could move) distributed (in order to be resilient to any loss of device or link capability) and independent of users position (as we believe that the UAV fleet cannot know the actual users position). The aim of this paper is then to propose an efficient solution. Our technique, introduced in [START_REF] Rautu | Crowd-based positioning of UAVs as access points[END_REF], is based on a Coulomb's law analogy. Of course, in dense areas, a single UAV cannot provide network access to a large number of users. The IEEE 802.11 standard, used in our network, offers several orthogonal channels to mitigate interference. So in this paper, we improve the coverage ratio with the use of multiple channels. For this purpose, we have to tackle two new issues: how can a drone select the channel to use? and how can a user select the drone to use as access point?

The remainder of this paper is organized as follows: Sect. 2 introduces a brief state of the art for channel allocation and optimization problems, Sect. 3 presents the system description, Sect. 4 introduces the interference management mechanism and its associated algorithms, Sect. 5 introduce the performance evaluation and Sect. 6 concludes the paper and proposes some potential extensions.

Related Works

Several works have already been done to take up similar issues.

In [START_REF] Han | Optimization of manet connectivity via smart deployment/movement of unmanned air vehicles[END_REF][START_REF] Senturk | Relay placement for restoring connectivity in partitioned wireless sensor networks under limited information[END_REF] the goal is to restore the network connectivity. Drones are used to fill the gaps serving as bridges between the disrupted infrastructures. UAVs are sent over the affected areas to interconnect terrestrial networks by relaying the messages between them. In [START_REF] Han | Optimization of manet connectivity via smart deployment/movement of unmanned air vehicles[END_REF], the authors use Delaunay triangulation to improve connectivity and in [START_REF] Senturk | Relay placement for restoring connectivity in partitioned wireless sensor networks under limited information[END_REF], they use a game theory approach to interconnect partitions of a network by using drones to drop relays in pre-computed spots. None of them can be used in dynamically changing networks due to the complexity needed to determine the positions of the temporary access points.

Several propositions exist in the literature as for how to optimally assign channels to reduce or eliminate the interference. Authors in [START_REF] Jahanshahi | On channel assignment and multicast routing in multi-channel multi-radio wireless mesh networks[END_REF] have studied the assignment of channels in the Multiple Radios, Multiple Channels Wireless Mesh Networks (MRMC WMN), and proposed a cross-layer mathematical formulation of joint channel assignment and multicast tree construction designed to minimize the total number of links by forming a multicast tree and thus minimize the total interference.

All these propositions cannot be used in our network as the approaches are centralized. A drone has to make decisions based on local information. In [START_REF] Nargesi | Efficient multicast and channel assignment in multi-channel wireless mesh networks[END_REF], the authors have proposed a distributed algorithm called Efficient Wireless Multicast (EWM) which builds a tree in which channels are assigned to the transmission for next hop in function of the used channels.

In [START_REF] Pang | Efficient data collection for wireless rechargeable sensor clusters in harsh terrains using UAVs[END_REF], the authors efficiently collect data and recharge sensors by the aim of drones. The network is separated in multiple clusters. Unfortunately, this publication offers an overview over a static sensor network and does not look into the optimal deployment of the UAVs.

The optimal placement of UAVs in order to cover targets on-the-ground is already researched form different perspectives.

In [START_REF] Wang | Wsn04-1: deployment for information oriented sensing coverage in wireless sensor networks[END_REF][START_REF] Martínez | Optimal sensor placement and motion coordination for target tracking[END_REF], the authors assume that the devices evolve in a 2D space. Therefore their problem is simplified because the coverage radius is fixed for each mobile devices. In the former article, the authors consider a mathematical model to maximize the amount of information collected based on a greedy approach and in the latter, the authors present a decentralized model for optimal positioning of sensors in order to track a target.

In [START_REF] Zorbas | Optimal drone placement and cost-efficient target coverage[END_REF], the authors aim to find UAV positions in order to minimize the number of drones used ensuring the surveillance of all the targets, by defining a mixed integer non-linear optimization models.

Nonetheless, the difference between the works presented above and ours is related to the constraints. Indeed, we consider that not all users have to be covered by a drone and introduce it as a constraint in our problem.

System Description

Scenario

The purpose of our work is to use a fleet of UAVs to set up a backbone network which provides communication means for a particular event such as a public gathering or a disaster situation when the traditional infrastructure is overloaded or wiped out. This network should be quick to deploy and inexpensive to implement, so that implicated parties (organizers, firefighters, public services) could use it rapidly. We will assume that the UAVs can communicate with each other, eg with the help of directional antennas. The primary objective of this work is to implement the best coverage with the minimal number of drones and with the most extended battery life possible for the UAVs.

Users position cannot be known with precision, so we will assume that we only know some specific points close to which users are more likely to be found (eg checkpoints).

Our main objective is to define the initial position of each UAV in order to maximize the number of associated users. A user can be associated to a drone if the reception power is strong enough and if the UAV can provide the required service.

Modeling Users Position

The scenario we focus on is depicted in Fig. 1. Let D be the set of available drones and U the set of users to cover. Each user u ∈ U has a fixed position at (x u , y u , h u ), where x u and y u represent the positions in 2D plane and h u the altitude of the user, fixed at 0. Each UAV, d ∈ D is free to move into a 3D plane, having the coordinates (x d , y d , h d ). The altitude of the drone is fixed at h d = 100 m. We assume that each drone has a fixed sensing range s. A user u found in this range is then detected and can be associated to a drone. The radius of the sensing range is computed based on the UAVs altitude and its directional antenna half beamwidth, in our model, s = 50 m.

We consider that a drone d ∈ D covers a user u ∈ U if the distance (1) between the drone and the user is

L u d ≤ s. L u d = (x d -x u ) 2 + (y d -y u ) 2 + (h d -h u ) 2 (1)
As already stated, we have already proposed an initial UAV placement strategy based on a statistical knowledge of users positions [START_REF] Rautu | Crowd-based positioning of UAVs as access points[END_REF]. For this purpose, we have introduced a Point of Interest or POI as a point close to which users are more likely to be found. We have also introduced p, the probability that a user is at most at a distance d of one of the N p POIs, named P 1 , . . . , P N p .

In Fig. 1, N p = 3 and p = 0.75. Distance dist has been arbitrarily chosen to match the sensing range of a drone.

Modeling UAV Behavior

UAVs are supposed to get as close as possible to a maximum number of users while preserving a minimal distance to each other (in order to mitigate interference). We have thus chosen to represent their interactions with the help of a model inspired by Coulomb's law [START_REF] Rautu | Crowd-based positioning of UAVs as access points[END_REF]. In this model, a user is described as a positive electric charge and a UAV as a negative electric charge. Drones are then attracted by users within their sensing range.

On the other hand, UAVs using the same channel have to repel one another to avoid interference. As can be seen in Fig. 2, if γ is the ratio between a user charge and a UAV charge, a low value of γ will induce a high ratio of "interfered" users (users within the range of several UAVs, that could then suffer from hidden terminal situations). Of course a high value for γ is more difficult to implement and could lead to a lower association ratio. Users are not affected by this force, being able to move freely. In our simulations, UAVs are activated above POIs and aim at reaching an initial position with an optimal coverage. Such a position must be reached as soon as possible to improve battery life.

Traffic Model

In this paper, we focus on an initial UAV positioning, seeking an optimal coverage. We will not study traffic scheduling, so we will not implement a sophisticated data traffic. To bring our model closer to reality, we will however assume that each user has the same communication needs, with limited resources available on each drone. Each UAV has a fixed capacity equal to κ = 50 users. In these conditions, more UAVs are needed to cover users concentrated in the same region. A typical example, when p values are high, users will gather tightly in a small area around the POIs.

Observed Metrics

Our objective is to provide the best network access to a large number of users. We will thus measure the following parameters. The coverage ratio is the ratio of users with a high enough reception power from at least one UAV. The interfered ratio is the ratio of users covered by more than one UAV using the same channel. The associated ratio is the rate of non-interfered users, associated with a UAV. The reception power is evaluated for each user with the help of a propagation model. The battery capacity will be described with a simple model.

Path Loss Model

UAVs are situated above the users at a constant altitude h = 100 m, and we assume that they are in line-of-sight so, we used the Friis transmission formula [START_REF] Friis | A note on a simple transmission formula[END_REF] to calculate the received power P R :

P R = P T G T G R c 2 (4πRf ) 2 (2)
Equation 2 allows acquiring a magnitude of radio power sensed by a receiver located at a certain distance of a transmitter, in free space. P T represents the transmission power, G T represents the transmitting antenna gain, and G R represents the receiving antenna gain. R and f symbolize the distance between transmitters and the used frequency respectively. The simulation parameters used are displayed in Table 1.

Energy Management

A simple energy model was also integrated into the model. The purpose is to estimate the mean lifetime of a drone participating in the deployment and operation of the network. To compute the energy consumption, we took into consideration only the main equipment embarked on an UAV such as the propeller motors, the CPU (central processing unit) and the Wi-Fi antennas.

A precise model of motor consumption depends on weather conditions (temperature, atmospheric pressure, winds), propeller diameter and pitch, internal resistance and motor efficiency. We choose to use the measurements presented in [START_REF] Dietrich | An empirical study on generic multicopter energy consumption profiles[END_REF]. The authors did an extensive study on the energy consumption in various discrete movement states of a UAV. The Wi-Fi consumption is also based on other experiments done in [START_REF] Zhao | Energy-efficient data gathering in high-voltage transmission line monitoring system[END_REF]. In this paper, the authors measured, in detail, the energy consumption for wireless nodes. The last piece of equipment took into consideration in our energy model is the CPU. We imagined that the UAVs would embark an on-board computer, like Raspberry Pi, that will run our model. In [START_REF] Kaup | Powerpi: measuring and modeling the power consumption of the raspberry pi[END_REF] the authors measured the power consumption of a Raspberry Pi based on different CPU utilization. A summary of different power consumption utilized in our simulation is given in Table 2. To determine the battery lifetime we used the following equation: T r = C b C e ǫ, where T r is the remaining battery lifetime in hours, C b the battery capacity in mAh, C e the current in the drone's load. ǫ takes into account the external factors that can affect the autonomy of the battery, being equal to 0.7 in our simulations.

A UAV replacement mechanism was already studied in [START_REF] Rautu | Maintaining a permanent connectivity between nodes of an air-to-ground communication network[END_REF]. When a drone's battery is about to be depleted, it will send a replacement request over the control plane. A replacement UAV will come alongside the depleted drone permitting the routing algorithm to adapt the routes. When traffic is rerouted, the end-life drone will depart to the control center to recharge its batteries.

Implementation Concerns

As mentioned, we will use IEEE 802.11 as our communication protocol between drones and between UAVs and users. To adapt our method to real life, Wi-Fi beacons and probe requests come in handy. Two different discovery approaches are proposed by the standard, passive or active scanning. In passive scanning, a station will scan all possible channels, one by one and listen to beacons. A beacon frame is sent by every access point to announce its presence.

For active scanning, a station still goes through each channel in turn, but instead of passively listening to the signals on that channel, It will send to the broadcast address a Probe Request management frame asking what network is available on that channel.

We use these two frames to implement the interactions in our network. As the access points mounted on UAVs send beacons to the broadcast address, other drones will receive them. By receiving the beacon on one of the directive antennas the UAV can approximately determine the direction of the neighboring drone. Moreover, as each UAV is equipped with GPS (Global Positioning System), we can use the beacon's timestamp to roughly determine the relative distance between drones.

Regarding the user detection, we use the probe requests that user equipment send to detect the presence in the drone's sensing range. After that, by using the information on signal and noise level, the UAV will try to determine the relative position of the user.

Interference Management

As presented in [START_REF] Rautu | Crowd-based positioning of UAVs as access points[END_REF], our network performs very well without having the exact users positions. But when taking into account the interference that can be produced by neighboring drones, we notice a slight decrease in association ratio. If two or more drones are needed to cover the same area, user communications are disturbed by the nearby drone, if they use the same radio channel. As already stated, it is possible to reduce interference by increasing γ (minimal distance between drones), but in the same time, the association ratio decreases as some users will not be covered anymore.

Introducing Multi-channel

Figure 3 depicts the percentage of associated, non-interfered users for several numbers of channels. We notice that our assumptions where founded. Of course, using the same frequency on each drone, creates massive interference with more than 50% of users interfered when p is high, as UAVs are trying to cover as many users as possible.

Fig. 3. User association ratio based on the number of channels

The overlapping can be managed by altering the ratio between repulsive and attraction forces. In fact, this allows us to state a trade-off between covering users and reducing interference.

Interference Management and Channel Switching

One of the improvements is related to the usage of multiple radio channels. As explained before, the system becomes interfered when multiple UAVs try to cover the same area. The main purpose is to offer an interference-free environment reducing the risk of frame collisions and reduced throughput. For this, we added a multi-channel reuse scheme to reduce interference.

Luckily, physical specifications for IEEE 802.11 standard allow for simultaneous operation of several orthogonal (non-overlapping) channels. As an example, in the 2.4 GHz Wi-Fi band, three channels can be used concurrently. As to IEEE 802.11a standard, introducing the 5 GHz band, a total of twelve orthogonal channels are set forth.

By deploying antennas which allow the usage of multiple channels and affecting different non-overlapping channels when they are located in proximity, the UAVs can provide services to user on-the-ground simultaneously with minimum interference. Therefore, the association capacity of our proposed network can be increased.

Drones are collaborating, when in the immediate vicinity, to change the used channel. Each one exchanges, with its neighbors, information on the used channels. Accordingly, they adapt their own to eliminate local interference. In the beginning, every drone starts with the same Wi-Fi channel. At each iteration all UAVs will search for potential neighbors. Once found, it will choose a channel not used by another UAV in the vicinity.

As simple as it may look, it is very efficient, as it provides complete decentralized channel management. The flip-flopping between channels is mitigated as only the newest interfering drone will change the channel, bringing the neighborhood to a stable state. Besides that, as only one drone can send messages over the same Wi-Fi channel at one time, there cannot be two or more UAVs changing the channel at the same time as they have to notify the neighbors and the associated clients.

Association Strategy

As a user may be within the range of several UAVs (each one using a different channel), the question of the channel (and thus the drone) to use arises.

When a new user is covered by several UAVs (still with available capacity), then it can be associated to any of them. We have studied three different UAV association strategies:

• The simplest strategy is to use any of them, for example the first one from which a beacon is received. Let us call this the random strategy. • The first real strategy is to implement a very simple load balancing algorithm.

Users can be transferred between UAVs to reach a fair share.

• The second one is to transfer a user from a UAV to another one only if this improves the received power for this user.

If the drone capacity can be clearly stated and depends only on the number of associated users, then the random strategy is relevant. However, if the back-haul link provided by the UAV is a bottleneck, then the simple load balancer seems more appropriate. Finally, if the bottleneck is the Wi-Fi capacity, then the last one can be thought as suitable as it aims at improving transmission conditions.

Through an extensive set of simulations, we have noticed that these metrics are not really sensitive to the choice of the algorithm. As depicted in Fig. 4, the ratio of covered users is merely the same. The reception power for the users, described in Fig. 5 is not really affected neither, mainly because of the UAV altitude (in our model, the reception power depends on the distance from the UAV, and thus on its altitude). Finally, we can see in Fig. 6 that for high values of p, the random solution is more energy efficient. The reason is that it avoids some UAV movements that would follow a client transfer. In this section, we evaluate the performance of our proposed schemes with the help of a simulator that we developed. The simulation parameters are presented in Table 3. As depicted in Fig. 7, using only one channel produces massive interference when the users are gathered around the POIs. In this cases, frame collisions and a reduced throughput can be expected. By taking into account these limitations we refined our simulator to include these new features. The interference that a user can sustain when multiple drones, operating on the same channel, are in its vicinity is also analyzed. To consider a user as interfered, the power received from the interfering AP (Access Point) must be higher than the receiving threshold, in our case T R = -65 dBm. All the following results are taking into account the interference. Only non-interfered users are counted as associated.

In Fig. 7 we compared the association rates for 500 users when using one, two and three channels. For low values of p, users are disseminated and the association rate is about 20%. For high values of p, because of the high user density, one channel is not enough. We see that the use of multiple channels dramatically enhances the performance for high values of p. We can also notice that for p close to 1, the number of channels needed is roughly equal to (3) where N u is the total number of users, N p is the total number of POIs in the network and κ is the maximum number of users that a drone can cover.

N C = N u N p * κ (3)
For example, using the values from Table 3 we obtain 2 channels, N C = 2. Even with this number of channels, we cannot reach a 100% user association rate for the following reasons:

• The users are not evenly distributed among POIs so that if a point of interest has 110 users, 10 of them will not be covered, with only 2 drones. • Two POIs can be close enough so that 2 channels are not sufficient, some users being interfered by neighbor UAVs.

Taking our model to a higher scale, we analyzed the performance when a larger number of users is present. We looked into the association rates by varying the number of channels available when n = 2000 and have to be covered with only 40 UAVs. In Fig. 3 we depicted the association ratio when p = 0, 25, p = 0, 5, p = 0, 75 and p = 1. When disposing of only a small number of channels it is better when the users are more dispersed than gathered around a POI. As users are concentrated in a small area, p is high, the number of drones that can cover them without creating interference will be limited by the number of channels. When the number of channels is sufficient to cover everyone, as users are more dispersed, the number of drones needed to provide a 100% association rate grows to be inefficient to deploy that many resources to increase the association rate.

Table 4 shows the energy consumption of a drone in our network. Based on the model described before, we compared the mean battery lifetime and the mean energy consumption. We evaluate the impact of the p value and the number of channels on the energy consumption. The number of channels does not impact the battery lifetime. However, as expected, the users position has a tremendous effect on battery depletion. For p = 1, i.e., users are all gathered around POIs, UAVs do not have to move a lot, which saves the battery lifetime, offering approximately 16 minutes of flight time. As users are more dispersed in the given area, the drones have to travel longer distances to find them, using more power. When p = 1.0, 30% more energy is consumed by the drone decreasing the mean lifetime, as we expected. We believe that it worth adding the antennas to the UAV, battery lifetime not being drastically impacted. Our proposed network performs remarkably well, the interference management and association schemes improve drastically the ground association rates offering better quality, not interfered transmissions and better resource allocation between users.

Conclusion and Future Works

We are currently working on optimizing drone placement taking into account our constraints. This optimization will allow us to observe the efficiency and performance of our self-deploying distributed network. User mobility is another critical point in our study. Being distributed and easy to implement, our solution should be efficient in a mobile scenario.

This paper proposed several novel enhancements to our network of UAVs used as access points. The network deployment is based on Coulomb's law, with users being represented as positive charges and drones as negative charges, being attracted by users. After detailing this model, we introduced two mechanisms that boost the user association ratio of the network, an interference management, and several association mechanisms. Allowing each drone taking part in the network to choose in a distributed manner the best channel to use by taking into account local interference, we point out that the non-interfered association rate of users increases remarkably. Furthermore, an analysis of energy consumption is done. Adding extra weight to drones (antennas, additional on-board computer) does not impact the energy consumption, overall battery consumption remaining reasonably acceptable, contributing to the feasibility of our proposed model in real life.

In this paper, users are evenly distributed among the POIs, they all have the same communication needs and all the UAVs have the same capacity κ. We plan to run more simulations to study the impact of these parameters.

We have shown that p is an important parameter, but its value is probably unknown in actual scenarios. We can imagine the use the UAVs discovering users position to dynamically compute an estimator of p that could be helpful to improve their behavior. Improving the performance for low values of p is also a challenge.

Finally, this work was dedicated to the initial placement of UAVs. The next step is to study how our system behaves with moving users.

Fig. 1 .

 1 Fig. 1. Our scenario's model

Fig. 2 .

 2 Fig. 2. Association ratio for high and low values of γ

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig. 4. Association ratio for different association strategy

Fig. 7 .

 7 Fig. 7. User association ratio comparison when adding multiple channels

Table 1 .

 1 Scenario parameters

	Characteristics	Value
	UAV transmission power P T d 27 dBm
	User transmission power P T u 27 dBm
	Transmission gain G T	5.2 dBi
	Reception gain G R	2 dBi
	Transmission frequency f	5150 Mhz
	Received power threshold T R -65 dBm

Table 2 .

 2 Energy consumption parameters

	Characteristics	Value
	One motor -hoovering	3 Amps
	One motor -in move	5 Amps
	One Wi-Fi antenna -idle/RX 300 mAmps
	One Wi-Fi antenna -TX	400 mAmps
	CPU -50% utilization	600 mAmps
	Battery capacity C b	6200 mAh

Table 3 .

 3 Simulation parameters

	Characteristics Value
	card(u)	500
	card(d)	10
	κ	50
	dist	50 m
	N p	5
	(L × W )	1 km × 0.8 km
	γ	2
	h	100 m

Table 4 .

 4 Mean energy consumption and mean battery lifetime

	p	Battery lifetime Energy consumption
	p = 1.0 16 mn	16 amps
	p = 1.0 12 mn	24 amps