Martin Rosalie
email: martin.rosalie@uni.lu

Grégoire Danoy
email: gregoire.danoy@uni.lu

Serge Chaumette
email: serge.chaumette@labri.fr

Pascal Bouvry
email: pascal.bouvry@uni.lu

Chaos-enhanced mobility models for multilevel swarms of UAVs

Keywords: ACO, Swarm of UAVs, Mobility model, Chaotic dynamics, Global optimization, First return map PACS: 05.45.-a, 02.90.+p

published or not. The documents may come

Introduction

The development and usage of Unmanned Aerial Vehicles (UAVs) has quickly increased in the last decades, mainly for military purposes. Nowadays, this type of technology is also used in non-military contexts, for instance for environment protection, by search and rescue teams, by fire fighters and police officers, or for environmental scientific studies. In order to increase their potential, swarms of UAVs are now envisioned. In this context the payload and sensors are shared between UAVs to address the limitations encountered when using a single UAV, such as autonomy, operation range and resilience. Although the technology for operating a single UAV is now mature, additional efforts are still necessary for efficiently taking advantage of UAVs swarms.

In this paper we address the problem of area coverage with a swarm of UAVs collecting data via sensors. This problem is addressed in the framework of the ASIMUT project funded by the European Defence Agency (EDA). The military aspect implies that an additional constraint has to be considered: the unpredictability of the trajectories eliminating sweeping or patrol algorithms. Although this paper is dedicated to UAVs, the problem remains similar for other unmanned systems like ground, surface or underwater vehicles. To solve this problem two main techniques have been proposed: online and offline path planning. Offline path planning consists in precomputing the flight plan of the UAVs. The main asset of this approach is that the trajectories of the UAVs are easily monitored from the Ground Control Station (GCS). However, this technique is not adaptive to any change of configuration during the flight: the scheduled path can be irrelevant by the time the UAVs execute it. On the other hand, online methods compute the trajectories of the UAVs at runtime. The advantages and drawbacks of online methods are the opposite of those of offline approaches: they are flexible and resilient but it becomes impossible to predict the trajectories of the UAVs. As a consequence, in this paper we will propose a mobility model that combines the assets of both online and offline methods.

The latter method relies on the Ant Colony Optimization method (ACO) introduced by Dorigo [START_REF] Dorigo | Optimization, learning and natural algorithms[END_REF], and more precisely on the work of Kuiper & Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] who adapted the ACO algorithm to the coverage problem for UAVs. Kuiper's mobility model uses repulsive pheromones to guide the UAVs over the area they have to cover. The UAVs share a map of virtual pheromones that indicates recently visited areas when high pheromone concentrations are present. The UAVs then have a higher probability to move to the least recently visited areas.

Our contribution proposed in this paper consists in an ACO-based mobility model for UAV swarms that uses a chaotic dynamical system. It addresses part of the questions raised by a survey on algorithm dynamics and complexity [START_REF] Zelinka | A survey on evolutionary algorithms dynamics and its complexity -mutual relations, past, present and future[END_REF]. Chaotic dynamics are the solution to a deterministic system with the following properties: the solution is bounded, globally time invariant and sensitive to initial conditions, and consequently, unpredictable on a long-term. Using such solutions in a UAV mobility model would thus permit to obtain deterministic but unpredictable trajectories. Our objective is here to study in detail the impact of different chaotic systems and parameters on the area coverage by the swarm. Recently, Zang et al. [START_REF] Zang | Applications of chaotic dynamics in robotics[END_REF] provide a detailed review of various applications of chaotic dynamics for mobile robots: several discrete and continuous dynamical systems are used to generate chaotic dynamics. We choose to use the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] as a basic system to explore a combination of chaotic behaviours from Ordinary Differential Equations (ODE) with Ant Colony Algorithm. This system is a reference in the literature, as one of the first ODE systems with a simple chaotic mechanism. The Ma system [START_REF] Cai | A new finance chaotic attractor[END_REF] has also been considered to illustrate the transition from the random part of an ACO algorithm to a chaotic one.

This paper is an extension of some of our previous work [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. More precisely, we additionally provide a complete description of the chaotic dynamics used, as well as a detailed analysis of the methodology applied to build a mobility model from first return maps. We extended our study with two additional chaotic mobility models, detailing the properties of the periodic orbits and their impact on the efficiency of the coverage. We also analyze in detail the influence of the periodic orbits on the exploring patterns they can generate. Finally the performance of the new mobility models is also studied via the metrics and compared to the models from literature.

This article is organized as follows. We first introduce the context and the problem definition. Then we describe the related works regarding the two main topics upon which we elaborate our contribution: Ant Colony Optimization and chaotic dynamics. In the third section, we present our chaotic mobility models based on differential equations systems. The next section contains a description of the integration of our best chaotic mobility model into an Ant Colony Optimization algorithm. In the two last sections, we describe our experimentations including the metrics and the statistics we produced. We finally give a conclusion and describe our future work.

Context and problem definition

In this section we first present the context of the work, i.e. the ASIMUT (Aid to SItuation Management based on MUltimodal, MUltiUAVs, MUltilevel acquisition Techniques) project. We then describe the tackled problem of area coverage with a swarm of UAVs.

ASIMUT project

The purpose of the ASIMUT project, supported by the European Defence Agency (EDA), is to improve the situation awareness of an operator through area coverage and detection of threats based on multi-sensor and multi-source data fusion [START_REF] Bouvry | ASIMUT project: Aid to situation management based on multimodal, multiuavs, multilevel acquisition techniques[END_REF] 1 . Information is delivered by heterogeneous swarms of autonomous UAVs flying at different altitudes. One of the objectives of ASIMUT addressed in this article is the efficient surveillance of an area by means of a swarm of UAVs (these UAVs collect data with their embedded sensors). It should be noted that the operation takes place in a military context where the unpredictability of the trajectories is mandatory to prevent interception. The Fig. 1 details all the components of the system including the heterogeneous swarms (High Level Coordination Swarm and Low Level Coordination Swarm) and the mission management entities. A real deployment of this model is planned in our roadmap.

Problem definition

This work focuses on the mobility management of a swarm of autonomous UAVs to maximize the coverage of a squared geographical area. In addition to the unpredictability constraint induced by the military context, the trajectories of the UAVs still need to be monitored from the Ground Control Station (GCS) located on the middle of one the area edges. It is indeed mandatory for supervisors and users of this type of system to know and anticipate the positions of their UAVs.

To summarize, the objective is to "maximize the coverage while ensuring unpredictable trajectories". Additionally, we intend to provide an adaptive method resilient failures or losses of UAVs. As a consequence, the problem we address here is at the edge of the path planning for UAVs and the autonomous distributed coverage. This is not a path planning problem because of the resilience constraint. We are thus clearly explore an original problem.

Thus, two characteristics have to be considered:

Coverage consideration As the main purpose of the swarm is to cover a given area, the UAVs have to synchronize their exploration by preventing the other UAVs to explore already visited areas. However, we do not consider that revisiting an already visited area is forbidden but it is not profitable to satisfy the coverage objective.

Unpredictability consideration Because of the non predictability constraint, articles giving the optimized solution of a coverage problem [START_REF] Galceran | A survey on coverage path planning for robotics[END_REF] have not been considered. Indeed, in such approaches the UAVs coverage pattern is too explicit to be used in a military context. We propose to formulate the problem as follows. 10 UAVs evolves on a square area with positions given by a couple of real number (x, y). The surface is discretized in a 100 × 100 grid. The UAVs can move according to three directions: ahead (A), 45 • on the left (L), or 45 • on the right (R). During the resolution, when a UAV reaches a cell, the latter is considered as visited and this, during the whole simulation. The objective function is to minimize the number of steps required to visit the whole area:

min t max traj t (x) f (traj t (x)) ( 1
)
where t is the number of steps, x = [(x 1 , y 1 ), . . . , (x 10 , y 10 )] the position of the 10 UAVs, traj t (x) the trajectories of the UAVs from 0 to t and f the percentage of area covered depending on the trajectories of the UAVs using the 100 × 100 cells.

Related work

In this paper we address the following problem: an area has to be visited regularly by UAVs in order to collect information. One of the constraints is that an observer should not be able to anticipate the reconnaissance pattern of the swarm. A solution proposed by Kuiper and Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] was to introduce a random process to prevent the UAVs trajectories from being predicted. In the same article the authors also proposed ACO algorithm and compared it to the random approach. These two approaches are presented in detail hereinafter. Similarly to the random walk, the chaotic walk is introduced by Iba & Shimonishi [START_REF] Iba | The origin of diversity: Thinking with chaotic walk[END_REF]. It only uses the logistic map to determine the next angular direction of a UAV. A map is an iterative application that neither diverge, nor converges to a point. This logistic map (x n+1 = αx n (1x n )) converges to a dense set of points, between 0 and 1. This map generates a chaotic dynamic when α = 4 and it produces x n ∈ [0 : 1] as an output. We already tested this chaotic mobility model for coverage purpose and the results were not satisfactory [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF].

Even though we do not consider the path planning approach, some research work on this topic is worth presenting because the authors proposed original techniques to include chaotic dynamics. We remind the reader that the drawback of such an approach is that requires to have some initial knowledge about the area to cover. As an example for UAVs, Curiac and Volosencu [START_REF] Curiac | Path planning algorithm based on arnold cat map for surveillance[END_REF] provide the waypoints for the UAVs based on the Arnold Cat Map for such chaotic path planning. In order to generate waypoints for autonomous robots, Pimentel-Romero et al. [START_REF] Pimentel-Romero | Chaotic planning paths generators by using performance surfaces[END_REF] propose to use chaotic attractors as Random Number Generators (RNGs) using Chua, Lorenz and multi-scrolls attractors. Volos et al. [START_REF] Volos | A chaotic path planning generator for autonomous mobile[END_REF] provide a true random bit generator from a multi-scroll attractor to provide waypoints to the UAVs. In these articles, the coverage percentage is used to evaluate the performance of the presented methods while it is widely admitted that homogeneity and regularity for instance are key aspects. Finally, these methods provide waypoints as coordinate in the area to be covered contrarily to the Kuiper & Nadjm-Tehrani method where it is the direction (left, ahead or right) which is given at each step for the UAVs (see details hereinafter).

The rest of this section is organized as follows; we first provide details about Kuiper & Nadjm-Tehrani mobility models. Then chaotic dynamics are introduced to bridge the gap between the theoretical aspects and the practical use we introduced in this paper.

Random-based mobility model

Kuiper and Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] introduce a random mobility model to obtain an unpredictable mobility model for their UAVs. This model gives a direction to the UAVs: ahead (A), 45 • on the left (L), or 45 • on the right (R). In this random approach, the next direction of the UAV depends on its previous direction and is chosen using the probability rules given in Tab. 1. Each UAV is initially considered to fly "straight ahead". This mobility model result in a Markov chain where Tab. 1 is the transition matrix. The Alg. 1 details the mobility model. In the original paper [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF], the authors use this mobility model as a basic mobility model that is combined with an ant colony algorithm to increase the coverage efficiency. As our final purpose is to increase the performance of an ant colony algorithm we would like to first find a better mobility model than this random one and then we plan to combine it with an ant colony algorithm as well.

Algorithm 1 Random mobility model 1: procedure Random current state ← "ahead" 3: loop:

4: random ← value ∈ [0; 1] 5:
current state ← action(random, current state) # see Tab. move according to the current state

Pheromone-based mobility model

This method is a distributed pheromone repel mobility model introduced by Kuiper & Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF]. While UAVs fly, they also deposit virtual pheromones that indicate recently visited geographical areas and that evaporate over time. 100 units of pheromone are deposed on visited area (the maximum amount is 100 units per area) and it decreases by 1 unit per simulation step. If there are pheromones in its neighbourhood, a UAV chooses its next direction using the probability indicated in Table 2. This probability depends on the total amount of pheromones sensed and on the amount of pheromone sensed around the UAV. lef t is the amount of pheromone sensed on the left of the UAV, ahead is the amount of pheromone sensed in front of the UAV and right is the amount of pheromone sensed on the right of the UAV; finally total = lef t + ahead + right. If there is no pheromone to guide the UAV, it will use the random action based on its previous action (see Tab. 1). For clarity, we name this mobility model ACO UAV and Alg. 2 details it.

Table 2: Pheromone action table for ACO UAV mobility model. lef t is the amount of pheromone sensed on the left of the UAV, ahead is the amount of pheromone sensed in front of the UAV and right is the amount of pheromone sensed on the right of the UAV; total = lef t + ahead + right.

Probability of action Left

Ahead Right

p L = total-lef t 2×total p A = total-ahead 2×total p R = total-right 2×total Algorithm 2 ACO UAV mobility model 1: procedure ACO UAV 2:
current state ← "ahead" 3: loop:

4: random ← value ∈ [0; 1] 5:
if no pheromone sensed in the neighbourhood then move according to the current state

Chaotic attractor and topological analysis

In the following section, we survey the state-of-the-art concerning the use of chaotic dynamics for optimization problems. We then present the tools used to study and describe chaotic attractors solution to differential equations system. In most cases, the authors use the logistic map x n+1 = αx n (1x n ) to do it. We assume that the reason for that is the simplicity of this equation combined with an easy-to-use output. For instance Iba & Shimonishi [START_REF] Iba | The origin of diversity: Thinking with chaotic walk[END_REF] use it to define a chaotic walk (referring to random walk ). The shape of the trajectories depends on the number of decimals to compute the next iteration. This underlines how chaotic processes are sensitive to the initial conditions as well as the implementation. Gong & Wang [START_REF] Gong | Chaos ant colony optimization and application[END_REF] also introduce a chaotic process in an ACO algorithm by testing various parameters for the logistic map to improve the algorithm's results. This map is also successfully used in an Ant Colony Algorithm in order to find best UAV path between threats [START_REF] Zhou | Chaotic differential evolution approach for 3d trajectory planning of unmanned aerial vehicle[END_REF][START_REF] Zhang | UAV path planning based on chaos ant colony algorithm[END_REF].

Although the logistic map is omnipresent, this chaotic map has some drawback. Arroyo et al. [START_REF] Arroyo | On the inadequacy of the logistic map for cryptographic applications[END_REF] detail some of the properties of the logistic map (including probability distributions) while the parameter α is varied (see also [START_REF] Boeing | Visual analysis of nonlinear dynamical systems: Chaos, fractals, self-similarity and the limits of prediction[END_REF] for visual representation of these chaotic dynamics). They underline that this equation exhibits some periodic behavior between the values of the parameters used in [START_REF] Gong | Chaos ant colony optimization and application[END_REF]. The same ranges of parameters are used by Li et al. [START_REF] Li | An improved chaotic ant colony algorithm[END_REF] to solve the Traveling Salesman Problem (TSP). These works highlight the fact that the logistic map provides chaotic dynamics but for a small range of parameter values. Thus the logistic map is a good start but does not permit to explore all the capabilities and richness of chaotic dynamical systems. Some recent work have investigated further the integration of chaotic systems in optimization algorithms. Ahmadi et al. use Ordinary Differential Equations (ODE) systems to combine chaos and the weed invasive algorithm [START_REF] Ahmadi | Chaotic invasive weed optimization algorithm with application to parameter estimation of chaotic systems[END_REF]. Gandomi and co-workers use various chaotic maps (tent map, Lozi map, etc.) to explore the performance of their firefly [START_REF] Gandomi | Firefly algorithm with chaos[END_REF], Particle Swarm Optimization (PSO) [START_REF] Gandomi | Chaos-enhanced accelerated particle swarm optimization[END_REF] and Gravitational Search algorithms (GSA) [START_REF] Mirjalili | Gravitational search algorithm with chaos[END_REF]. Finally, several maps including the Lozi, Burger and Sinai maps have been used to tune a PSO algorithm with six different chaotic systems used as Chaotic Pseudo Random Number Generators (CPRNGs) [START_REF] Pluhacek | Chaos particle swarm optimization with eensemble of chaotic systems[END_REF]. The drawback of the last methods with chaotic maps is that there is no advanced dynamical tool to compare these chaotic dynamics from maps as it is the case for chaotic attractors from ODE with the topological analysis. Here we thus focus on chaotic dynamics of chaotic attractors without the aim to build CPRNG.

Before going into details, we first describe the topological characterization method [START_REF] Gilmore | The topology of chaos: Alice in stretch and Squeezeland, 2nd Edition[END_REF]. From a chaotic attractor obtained by solving an ordinary differential equations system, several tools dedicated to nonlinear analysis are used: Poincaré section and first return map; these will be defined later. The output of the topological characterization method is a template describing the topological properties of the chaotic attractor. In one of the intermediate steps of this method, we obtain a first return map describing the chaotic dynamics with a step by step process. We will use this first return map to replace the random part of a mobility model. In the following section we also describe the properties of chaotic attractors in order to use them to build efficient mobility models.

Introduced in 2001 [START_REF] Ma | Study for the bifurcation topological structure and the global complicated character of a kind of nonlinear finance system (i)[END_REF][START_REF] Ma | Study for the bifurcation topological structure and the global complicated character of a kind of nonlinear finance system (ii)[END_REF], the Ma system is a three differential equations system which purpose is to model financial interactions: it details the links between the interest rate (X), the investment demand (Y ) and the price exponent (Z). A dynamical analysis of this system [START_REF] Cai | A new finance chaotic attractor[END_REF] permits to show that for a given set of parameters, there is a chaotic attractor as a solution of this system:

     Ẋ = -aX + Z + XY Ẏ = 1 -bY -X 2 Ż = -cZ -X . (2)
The structure of this dynamical system is similar to the Lorenz one [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF]. In order to highlight the symmetrical structure of the Ma system, we reorganized the variables X, Y and Z to have a structure close to the Lorenz system with respectively y, z and x:

     ẋ = -cx -y ẏ = x -ay + yz ż = 1 -bz -y 2 .
(3) As a consequence, the symmetry is defined as follows:

γ • f ( x) = f (γ • x) , (4)
where γ is a rotation symmetry

γ =   -1 0 0 0 -1 0 0 0 1   . (5)
For the parameter values a = 0.00001, b = 0.00001 and c = 1.063, the symmetry contained in this system leads to a particular structure where the solution is the chaotic attractor A (Fig. 3). The solution is obtained numerically using a Runge-Kutta 4th order method. To satisfy the convention orientation used to analyse chaotic attractors bounded by genus-1 torus [START_REF] Rosalie | Systematic template extraction from chaotic attractors: I. genus-one attractors with an inversion symmetry[END_REF][START_REF] Rosalie | Systematic template extraction from chaotic attractors: II. genus-one attractors with multiple unimodal folding mechanisms[END_REF], we performed a rotation to obtain a flow evolving clockwise. The attractor A is thus displayed accordingly in Fig. 3. We already performed a detailed analysis of a similar chaotic attractor for the same system [START_REF] Rosalie | Composition laws of branched manifold for description of chaotic attractors bounded by high genus torus[END_REF]. Compared to this previous study, there is a change in the parameter values where c = 1.03 instead of c = 1.063. The purpose of this previous study was to highlight the symmetrical structure of this attractor by performing its topological analysis. This method permits to detail the topological properties of an attractor using nonlinear tools and topological invariants (see [START_REF] Gilmore | The topology of chaos: Alice in stretch and Squeezeland, 2nd Edition[END_REF] for a description of this method).

We now present the tools used to performed this topological characterization. First of all, the Poincaré section is a tool used to obtain a discrete representation of a continuous flow. Such a section is defined as follows:

P A = {(y n , -z n ) ∈ R 2 | -x n = 0.7, -ẋn < 0} . (6)
This Poincaré section provides a sequence of values describing the solution values after a topological period (Fig. 3). From this sequence, it is possible to build a first return map with ρ n ∈ [0; 1] that is a normalized value of the solution in the Poincaré section (see [START_REF] Rosalie | Systematic template extraction from chaotic attractors: I. genus-one attractors with an inversion symmetry[END_REF][START_REF] Rosalie | Systematic template extraction from chaotic attractors: II. genus-one attractors with multiple unimodal folding mechanisms[END_REF] for conventions). For the attractor A, the associated return map contains three branches ordered as follows: one decreasing, one increasing and one decreasing (Fig. 4). The extra structures on the left side with decreasing and increasing branches are due to a projection limitation. These extra branches do not influence the dynamical analysis (see [START_REF] Rosalie | Composition laws of branched manifold for description of chaotic attractors bounded by high genus torus[END_REF] for details). This small modification of the parameter value c does not influence the main results obtained during the topological characterization [START_REF] Rosalie | Composition laws of branched manifold for description of chaotic attractors bounded by high genus torus[END_REF] because the dynamical properties are the same according to the first return maps (same periodic points representing orbits with the same topological properties). For the topological characterization method, each branch represents a set of points with the same dynamical behavior. Thus a label (the symbolic dynamic) is assigned to each point of the solution in the Poincaré section according to its branch in the first return map. Using the first return map, it is possible to extract some specific periodic points. These points represent periodic orbits in the flow of the attractor A and they are considered as the skeleton of the attractor A. Even if there is an infinite number of orbits, we present only the orbits with a period lower than four (Fig. 4) these orbits are more frequently visited. One of the properties of the flow of a chaotic attractor is that the trajectory in the phase space evolves by visiting alternatively all the orbits of the attractor. As these orbits are unstable, the flow follows them for a moment and moves to another one. This behavior explains why there are periodic points in the first return maps even if the solution is chaotic. We remind the reader that by definition chaotic can be use to define a deterministic dynamic with the following properties: aperiodic behavior, sensitivity to initial conditions and globally time invariant. These orbits with a low period are more often visited that orbits with an higher period. We obtain these periodic points by collecting their position in the first return maps when the solution passes through the Poincaré section; if it requires p crossovers through the Poincaré section to come to the same point, thus it is a period p orbit.

Here stops the analogy with the topological characterization method because in the next sections we will only use the transition between the branches without the topological implication induced to the flow of the attractor A. We will also use the periodic points to build specific mobility models knowing that these points are regularly visited with a given order. In the present paper, we will use the Poincaré section, the first return map and the periodic points to replace the random process with chaotic dynamics.

Developing mobility models from chaotic dynamics

In this section, we introduce a new methodology to build mobility models from chaotic dynamics. We not only present the models adapted for the coverage problem but we also detail the mechanisms to build them depending on the addressed mission. We will first introduce three new mobility models that use chaotic dynamics only to replace the random part of a ant colony algorithm (i.e. when there is no pheromone to guide the UAVs). We then propose the CACOC (Chaotic Ant Colony Optimization for Coverage) algorithm that combines the best previously introduced chaotic mobility models with an ant colony algorithm.

Mobility models based on the Ma system

We choose to use the Ma system because the classical partition of its return map permits to have similar transition rules as those of the Kuiper & Nadjm-Tehrani random mobility model. In the remainder of this article, both terms mobility model and algorithm are used interchangeably. To build our mobility model, we consider that the three branches of the Fig. 5 are equivalent to the three possible directions the UAVs can take (Tab. 1 p. 7). From this first return map, we propose to build a mobility model named MAMM (Ma Mobility Model).

First partition for the Ma system: MAMM

The partition of the first return map (Fig. 4) mainly corresponds to the classical partition made in the topological characterization method where there are as many parts as the number of branches:

• if ρ n < 0.15, there is a decreasing branch on the left associated with the symbol L;

• if 0.15 ≤ ρ n < 0.8, there is an increasing branch in the middle associated with the symbol A;

• if 0.8 ≤ ρ n , there is a decreasing branch to the right associated with the symbol R.

The pseudo-code that implements MAMM is detailed in Alg. 3. current state ← "ahead" 3: loop:

4:
ρ ← next value in the first return map (Fig. 5)

5:
if ρ < 0.15 then current state ← "left"

6:

else if ρ < 0.8 then current state ← "ahead"

7:

else current state ← "right"

8:

move according to the current state As the first return map is ρ n+1 = f (ρ n ), for each point of the map with the symbol L (i.e. with an abscissa lower than 0.15) the ordinate value is lower than 0.8 (Fig. 5). In that case, this means that ρ n+1 is lower than 0.8. Thus the two next symbols that can be possible after the symbol L are L and A. With the same reasoning, a symbol A can be followed by L, A and R and finally, the symbol R can be followed by A and R. Thus, the partition made for MAMM permits to reproduce the transitions allowed by the Random model (Tab. 1).

We previously mentioned that the periodic points are associated with periodic orbits. As we have now assigned a symbol to each part of the first return map we are able to translate these periodic orbits in terms of sequences of symbols. The periodic points of Fig. 5 with a period lower than four define:

• three orbits of period 1: L, A and R;

• two orbits of period 2: LA and RA;

• four orbits of period 3: LLA, LAA, RRA and RAA.

Here we remind that solving the Rössler system give the trajectory evolving on the attractor in the phase space. This trajectory collapses momentarily to an orbit before diverging to another, thus the symbols of the reached orbit will be repeated successively several times. For instance, when the orbit of period 1 with the symbol A is reached during the simulation, it will result in a series of A (AAAAA . . . ), that is to say, the UAV will fly following a straight line. This mechanism is the same for each orbit. In terms of behavior for our UAVs, these series of symbols are associated to mobility patterns that will occur during the simulation:

• straight line: A;

• sharp turns: L and R;

• turns: LA and RA;

• wide turns: LAA and RAA;

• loops: LAL and RAR.

According to the possible patterns that can be obtained using MAMM we suppose that it will not perform well in terms of coverage because of the numerous ways to make turns and loops (see section 5 for quantitative results). These turns might prevent a good exploration by letting the UAVs hover in the same area.

Second partition for the Ma system: MAMM2

To overcome the bad coverage issue of the MAMM model, we propose another mobility model based on the same first return map but using a different partition. MAMM2 (Alg. 4) obeys the following partition of the first return map (Fig. 6):

• if ρ n < 0.1, then the symbol is L; Compared with Fig. 5, the periodic point labelled R, A and R (Fig. 5) are respectively associated to the symbols R, A and A. This orbit is thus labeled RAA and its corresponding pattern for the UAV movement is a large turn.

• if 0.1 ≤ ρ n < 0.9, then the symbol is A;

• if 0.9 ≤ ρ n , then the symbol is R.

These values are chosen with respect to the positions of the periodic points.

Algorithm 4 MAMM2 mobility model

1: procedure MAMM2 2:
current state ← "ahead" 3: loop:

4:
ρ ← next value in the first return map (Fig. 6)

5: if ρ < 0.1 then current state ← "left" 6:
else if ρ < 0.9 then current state ← "ahead"

7:

else current state ← "right" move according to the current state We chose these separation values (0.1 and 0.9) because it includes the periodic points corresponding a periodic orbit with a period equal to 1 in the part where the symbol A, for ahead, is assigned. Indeed, this new partition permits to remove the sharp turns and the loops from the list of patterns that the MAMM2 mobility model can produce. These changes are summarized in the Tab. 3. They are definitively more suitable for the coverage purpose (see section 5 for quantitative results).

Table 3: Patterns induced by the partition of the first return maps between MAMM (Fig. 5) and MAMM2 (Fig. 6). The partition of a return map influences the topological entropy of the system (i.e. the number of distinct sequences obtained for a given partition with the same number of symbols in each sequence). The topological entropy of MAMM2 is lower than the entropy of MAMM. The entropy is maximum when the partition matches the branch orientation partition.

MAMM

Mobility model based on the Rössler system

The three previous models (Random, MAMM and MAMM2) are symmetric because of the properties assigned to the symbols L and R: the probability are the same for L and R moves. In order to improve the coverage and to increase the unpredictability of the trajectories of the UAVs, we introduce an asymmetric mobility model. As we have done with the Ma system, we consider the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] to introduce another mobility model, that we call CROMM (Chaotic Rössler Mobility Model).

The Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF] is a set of three ordinary differential equations:

     ẋ = -y -z ẏ = x + ay ż = b + z(x -c) , (7)
where x, y, and z are the variables and a, b and c are the parameters. The solution of this system is a chaotic attractor B for the parameter values a = 0.1775, b = 0.215 and c = 5.995 (Fig. 7). In a previous paper [START_REF] Rosalie | Templates and subtemplates of Rössler attractors from a bifurcation diagram[END_REF] we performed the topological analysis of this attractor. Here we provide a summary of this process to obtain the mobility model. First of all, we build the Poincaré section

P B = {(y n , z n )|x n = 0, ẋn > 0} . ( 8
)
so as to obtain a first return map with ρ n ∈ [0; 1] the normalized value of y n in the Poincaré section. This first return map (Fig. 8) is made of two branches, an increasing branch followed by a decreasing branch. As there is less branches than the number of symbols we want to assign, we decided to split the first return map into equal parts. Thus, we made the following partition:

• if ρ n < 1/3 then the symbol is R;

• if 1/3 ≤ ρ n < 2/3 then the symbol is L;

• if 2/3 ≤ ρ n then the symbol is A.

This partition has been selected after an empirical evaluation of these symbols' combination on this first return map with equal partition. The CROMM mobility model pseudo-code is given (Alg. 5).

Algorithm 5 CROMM mobility model

1: procedure CROMM 2:
current state ← "ahead" 3: loop:

4:
ρ ← next value in the first return map (Fig. 8)

5: if ρ < 1 3 then current state ← "right" 6:
else if ρ < 2 3 then current state ← "left"

7:

else current state ← "ahead"

8:

move according to the current state

0 1 3 2 3 1 0 1 3 2 3 1 R L A R L A ρ n+1 ρ n period 1 period 2 period 4 L A R A Figure 8
: First return map of the Rössler attractor (Figure 2). This map is partitioned in three parts indicating the UAV directions for the CROMM mobility model: L (left), A (ahead) and R (right). Orbits of period 1, 2 and 4 respectively illustrate the patterns A: straight line, AR: turn and RALA: serpentine trajectory transition taking place between the periodic points highlighted in dark grey.

The shape of the first return map presented in Fig. 8 indicates the possible direction changes of the UAVs. From a previous ahead direction (A), a transition to A, L and R is possible. From a previous L direction, transition to A is the only possible one. Finally, from a previous R direction, only transition to L or A is possible.

The first return map (Fig. 8) contains the periodic points of the periodic orbits of the attractor B. These periodic points are associated with the following labels and patterns:

• period 1: A for straight lines;

• period 2: RA for turn;

• period 4: LARA for serpentine trajectories.

Based on some experiments that we have conduced, we believe that this absence of symmetry combined with the serpentine trajectory and straight line behaviors are in favor of a better coverage method (see section 5 for quantitative results).

CACOC model: Chaotic Ant Colony Optimization for Coverage

To satisfy all our problem constraints (good coverage and unpredictable trajectories), we combine the pheromone method [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] with CROMM (the best chaotic mobility model among those that we have proposed as it will be shown in the following section: experimental results). In the ACO UAV mobility model, random processes are used for any iteration: either the UAV randomly chooses the direction because there is no pheromone, or it includes a probabilistic process to chose a direction accordingly to the pheromone sensed. We propose to replace the use of a random process by the first return map value (Fig. 4 where ρ n+1 = f (ρ n ) is displayed) as we have done for CROMM.

When there is no pheromone to guide the UAV, the ACO UAV mobility model uses a random process defined in Random model (see Tab. 1 p.7). We replace this random process by CROMM to explore a non visited area by our CACOC model. If there are some pheromones the ACO UAV mobility model, the UAV chooses its next direction with probabilities. For CACOC, we propose to replace the random choice with the next value in the return map using the probabilities of Tab. 2 (p L , p A and p R ) as partition for the first return map.

• if ρ n < p R then direction is right;

• if p R < ρ n < p R + p L then direction is left;
• else the direction is ahead.

As detailed in Alg. 6 as well as in Fig. 9, the successive return map values (ρ) are used to choose the next direction for each step in the CACOC model. The use of pheromones permits to have a system that is resilient to failures or losses of UAVs: if a UAV is out of the system, it will no more deposit pheromones.

Algorithm 6 CACOC mobility model

1: procedure CACOC 2:

current state ← "ahead" 3: loop:

4:
ρ ← next value in the first return map (Fig. 8)

5:
if no pheromone sensed in the neighbourhood then 6: current state ← CROM M (ρ) # see Alg.

if ρ < p R then current state ← "right" # see Tab. 2 for p R 9: else if ρ < p R + p L then current state ← "left" # see Tab.

Experimental results

The keystone of our work is the replacement of the random part of an algorithm by chaotic dynamics so as to increase its performance. First, we have introduced three new mobility models based on chaotic dynamics and associated tools coming from nonlinear analysis. Second, we have picked the best one to obtain an ant colony algorithm in which we have replaced the random parts: in direction selection process when there are pheromones to guide the UAV or not. In this section, we will first present the tools used to evaluate several aspects of the coverage performance of the algorithms. The setup used for the experiments, including the different metrics, is presented in the next subsection, and the experimental results are presented and analyzed thereafter. The second part of this section is dedicated to the analysis of the mobility models without ACO. From these experiments, we select the best mobility model and include it in a mobility models (CACOC). Finally, the last part details the analysis of the mobility models with ACO including random process or chaotic dynamic.

Experimental setup

The simulation area is a 100 m × 100 m square, divided in square cells of 1 m × 1 m. The UAVs have a constant speed of 1m/s and they all depart from the base station located in the middle of the bottom edge of the area, i.e. position (50,0). At each simulation step, each UAV can do one of the following three actions: (1) go ahead: the UAV keeps the same direction; (2) go left: the UAV turns left with a 45 • angle; (3) go right: the UAV turns right with -45 • angle. To prevent collisions between the UAVs they all have non equal flight altitudes [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF].

Each of the six mobility models has been evaluated based on a swarm of 10 autonomous UAVs equipped with wireless communication capabilities. We ran 30 independent simulations of 7000 steps each to obtain significant, results and we applied statistical tests in order to compare the four models. The setup information are summarized in Tab. 4. The mobility models using chaotic dynamics are deterministic. As a consequence, if we run them 30 times to compare the metrics, the results will be the same. We then need to introduce variability to perform a statistical analysis of our models. We thus change the initial conditions of the systems. We add a random value ε ∈ [0 : 0.1[ to obtain the following initial conditions: Ma system: (x 0 , y 0 , z 0 ) = (-1 + ε + 0.1 × Id, 1, 1) Rössler system: (x 0 , y 0 , z 0 ) = (-0.4 + ε + 0.1 × Id, 0, 0)

where Id = 1, 2, . . . , 10 for the ten UAVs. This mechanism prevents having the same initial condition for all UAVs. This change is enough to obtain non identical maps for each UAV due to the sensitivity of chaotic dynamical systems to the initial conditions. Even though, the map shape will be the same but the sequence of visited points will be different. Consequently, the series of actions of the UAVs will be different while following the same chaotic dynamical properties. Finally, we implement our models using Python programming language and we provide some examples of patterns obtained from these mobility models (Fig. 10). For MAMM, there are a lot of loops. The partition modification in MAMM2 permits to avoid loops but there are a lot of straight lines. For CROMM, the purple trajectory is a straight line followed by a serpentine pattern while the green line ended by a serpentine pattern followed by a large turn.

Metrics used for the evaluation of the mobility models

In order to evaluate the performance of the different models in terms of area coverage we have used the following metrics from [START_REF] Danoy | Connectivity stability in autonomous multi-level UAV swarms for wide area monitoring[END_REF]: coverage, recent coverage ratio, fairness. These are detailed hereinafter and Fig. 11 shows the linear regression we introduce to compare these metrics. This figure also underlines the importance of the linear regression used by also showing the evaluation of the metric for one single run of each mobility model we implemented.

Coverage. The coverage is the portion of the total area visited during the whole simulation. The coverage value change during the whole simulation. To have a representative value of the coverage, we compare the coverage value after 7000 steps for each model. This indicates how efficient the models are to visit the total area. We also want to evaluate the first steps of each model in order to compare their initial behaviour while there is no guideline for UAVs. This will be done by extracting the slope of a linear regression (a × x) considering the 500 first steps.

Recent coverage ratio. This metric, introduced in [START_REF] Danoy | Connectivity stability in autonomous multi-level UAV swarms for wide area monitoring[END_REF], represents the percentage of coverage during the last 100 iterations. These 100 steps correspond to the pheromones evaporation time used in ACO algorithms. We exclude the results of the 100 first iterations of the simulation to compute the mean value of the recent coverage. For 10 UAVs over a 100 × 100 grid with pheromone duration equal to 100, the optimal value is 10% because (10 × 100)/(100 × 100) = 0.1.

Fairness. The fairness evaluates if all cells are regularly and equally visited. This is measured by the standard deviation of their respective number of visits [START_REF] Schleich | UAV fleet area coverage with network connectivity constraint[END_REF]. To evaluate the fairness during the whole simulation, we perform a linear regression (a × x + b) using the last 6500 steps. This balances the coverage initial slope that only evaluates the initial UAV trajectories. formance over the 4000 first steps (we remind that we used 7000 steps for the experimentation). Linear regressions are used to evaluate the algorithms: the 500 first steps for the slope of coverage and the rest of the simulation (6500 steps) for the slope of fairness. The average value of the recent coverage is used for the whole simulation except the 100 first steps.

Comparison of mobility models without ACO

The experimentations are performed for the four models using the above four metrics. We summarized the collected data in Fig. 12 where a box plot and a smooth density evaluation underline the statistical properties of the model over the metrics. For the metrics: "Coverage after 7000 steps", "Slope of coverage" and "Recent coverage", the best value is the highest. Conversely, the lower the "Slope of the fairness" the better the result.

For the four metrics studied in this work, the 30 results do not follow normal distributions. As a consequence we performed a Kruskal-Wallis rank sum test [START_REF] Kruskal | Use of ranks in one-criterion variance analysis[END_REF] to make an unpaired multiple comparison of our model. All the results for the four metrics have statistically significant differences with 95% confidence between any two algorithms (i.e., with a p-value < 0.05).

The first result that can be noted is that the MAMM is worse than MAMM2 except for the "Recent coverage" metric. The average value of the latter metric is 8% of the total area (we remind that the optimal value is 10%) while it is 7.8% for the random mobility model. Using the chaotic dynamics (partition of first return map) seems to be efficient to increase the revisit time even if the model is not better regarding the other metrics. On top of that, the standard deviation of the three other metrics is twice as much as it is for the random model. Thus, MAMM fails to fulfill our main objective that was at least to reproduce the performance of the random model by using chaotic dynamics. However, if we compare the Random model with MAMM2, we realize that, in average, the latter outperforms the random one for the four metrics (MAMM2 also outperforms MAMM for all metrics). However, the standard deviation of the metrics of Random and MAMM2 mobility models are roughly the same. MAMM2 permits to obtain an average value of coverage higher than 95% for most of the simulations while it is below this 95% threshold for the random model. The recent coverage is higher with an average value around 8.5%. Our hypothesis concerning the impact of the coverage the reduction of small turns and the increase of straight lines is thus validated. Finally, CROMM outperforms MAMM2 on all the metrics except for the "Slope of the fairness of the coverage". However, for the latter, the standard deviation is low compared to the others. This standard deviation is also lower for the other metrics especially for the coverage, where the average value is about 97.5%, and the recent coverage with an average value equal to 8.8%.

As a consequence, the best chaotic mobility model candidate to replace the random mobility model is CROMM (Chaotic Rössler Mobility Model). We also obtain some results for the coverage problem from the chaotic dynamics point of view. First, from the same first return maps, the reduction of the number of sequences associated with turns increases the coverage performance. Second, our hypothesis concerning the asymmetry is confirmed by the excellent performance of CROMM.

Comparison of mobility models with ACO

Using ACO permits to increase the coverage performances because it takes advantage of the pheromones. With the same setup, we realized experimentations using CACOC and ACO UAV mobility models to evaluate and compare their performances. The settings are the same as reported in Tab. 4. As the CACOC mobility model is deterministic, we use the Rössler initial conditions of Eq. ( 9). All the results for the four metrics have statistically significant differences with 95% confidence between any two algorithms (i.e. with a pvalue < 0.05). The results are summarized Fig. 13. Even if this figure does not contain the results of Fig. 12, it can be noted that the metrics evaluation of ACO UAV and CACOC are better than the CROMM ones.

First of all, regarding the coverage value at the end of the simulation, the average value of ACO UAV is 97% and CACOC is 1% better with an average value of 98%. We also observe that the standard deviation is lower for CACOC. Regarding the slope of coverage, the result are better for our CACOC mobility model (the average value is 15% better). Indeed, the way we used chaotic dynamics to obtain the next direction according to the quantity of pheromones (Alg. 6) is very efficient. We remind the reader that it takes into account the pheromones sensed around the UAV to make the partition of the first return map. This combination permits to spread the UAVs at the beginning of the simulation.

Concerning the slope of fairness obtained by linear regression, the average values are equal. However, the standard deviation of the CACOC mobility model is lower by 45%. This permits to say that from a global point of view, the fairness of our coverage has less variability. Finally, considering the recent coverage, the results show that with an average value higher than 9.2% of the total surface, CACOC is better than ACO UAV (only 9.05% of total surface). As the optimal value that is 10% of the total surface, CACOC represents 92.2% of the optimal value that is better than the 90.5% of ACO UAV.

From these metrics, we can conclude that the results obtained for our CACOC mobility model are better than those of the ACO UAV algorithm provided by Kuiper & Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF]. It increases the performance while it reduces the variability of the algorithm. By maintaining the unpredictability of the trajectories using chaotic dynamics we also suppress a drawback of the random processes: the variability of the results.

Conclusion

In this article we have tackled the problem of area coverage by a swarm of UAVs in a military context. We have exploited the deterministic properties of chaotic dynamics to produce mobility models that permit to monitor the swarm behavior from a Ground Control Station while generating UAVs movements that are unpredictable for any external observer. We have first empirically demonstrated that our Chaotic Rössler Mobility Model (CROMM) improves the system behaviour compared to the Kuiper & Nadjm-Tehrani model [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] when there is no pheromone to guide the UAVs.

We also studied the usage of another chaotic dynamical system (the Ma system) to highlight the capabilities of partitioning a first return map used to define the movement patterns of the UAVs using its periodic points. We have empirically shown that reducing the number of orbits associated to turns and increasing the number of orbits associated to straight lines permits to enhance the coverage performance of the algorithm. We proposed ways to tune the partitioning of the first return map of a dynamical system resulting from Ordinary Differential Equation (ODE) systems: by tuning the proportion of the partition in accordance to the periodic points it contains and their relative organisation. As the periodic points are the skeleton of the dynamical properties of the system, their positions in the map have to be identified to build efficient exploratory movements from a first return map.

Then, we have proposed the Chaotic Ant Colony Optimization to Coverage (CACOC) algorithm that combines an Ant Colony Optimization approach (ACO) with a chaotic dynamical system. CACOC has provided better results than the state-of-the-art approaches on the coverage problem while satisfying the following constraints: the mobility is decided online in a deterministic and reproducible way but the trajectories of the UAVs remain unpredictable. CACOC is an original approach that combines the assets of offline and online methods. To summarize, we have shown in this paper that the usage of an appropriate and well-parameterized chaotic first return map to define the possible UAV direction changes permits to improve the global swarm coverage performance. In addition, when combined with a nature-inspired techniques, i.e. ant colony optimization, it permits to improve the state-of-the-art.

As future work, in the frame of the ASIMUT project, we plan to use both chaotic mobility models for multilevel swarms of UAVs: CROMM for the HLCS and CACOC for the LLS. In terms of algorithm perspectives, we propose two maps and three partitions for the coverage problem. The same methodology can be applied to replace the random part of other optimization algorithms such as Firefly Algorithm or Particle Swarm Optimization algorithms and to compare these partitions using topological analysis tools. Finally, we plan to investigate if it is possible to find another chaotic mechanism (e.g. return maps shape) that will make our system even more efficient. As a consequence, we are also intend to use even more complicated dynamics from well-known ODE systems such as Lorenz, Chen or Chua systems. We also plan to compare the performance of our algorithm, varying a parameter of the Rössler system, with the corresponding partition of the return maps [START_REF] Rosalie | Templates and subtemplates of Rössler attractors from a bifurcation diagram[END_REF].

Figure 1 :

 1 Figure 1: ASIMUT components where swarms of UAVs (High Level Coordination Swarm & Low Level Swarms) collect data and send them to the GCS (Ground Control Station) to be handled by the Intelligence Center.

Figure 2 :

 2 Figure 2: Ten UAVs using pheromones to cover an area of 100 × 100 cells. The pheromones disappear over time and are represented using a green scale for each cell.

1

 1

6 :

 6 current state ← action(random, current state) # see Tab.

1

 1

8 :

 8 current state ← action(random, pheromone sensed) # see Tab.

2

 2

Figure 3 :

 3 Figure 3: A is the attractor solution of the Ma system (3) for the parameter values a = 0.00001, b = 0.00001 and c = 1.063. The solution is presented in the phase space (-x, y, -z).

Figure 4 :

 4 Figure4: First return map (ρn+1 = g(ρn)) to the Poincaré section[START_REF] Cai | A new finance chaotic attractor[END_REF] for the attractor A. The periodic points indicate how the periodic orbits structure the flow of the attractor A. In dark grey, we highlight the period 3 orbit composed by three periodic points with the following structure x = (ρn, ρn+1): 1 = (0.945, 0.637), 2 = (0.637, 0.845) and 3 = (0.845, 0.946). The transition using the bisector lines detail the links between the periodic points.

Figure 5 :

 5 Figure5: First return map for MAMM with the partition associated to each branch: L for the first decreasing branch, A for the increasing branch and R for the second decreasing branch. Compared with Fig.4, the periodic points labelled 1,2 and 3 are respectively associated to the symbols R, A and R. Thus this orbit is labeled RAR and its corresponding pattern for the UAV movement is a loop.

Figure 6 :

 6 Figure6: First return map for MAMM2. The partition of the map does not correspond to the branches. Compared with Fig.5, the periodic point labelled R, A and R (Fig.5) are respectively associated to the symbols R, A and A. This orbit is thus labeled RAA and its corresponding pattern for the UAV movement is a large turn.

Figure 7 :

 7 Figure 7: Chaotic attractor B solution of the Rössler system for the parameter values a = 0.1775, b = 0.215 and c = 5.995. The Poincaré section (8) is indicated with the orientation to obtain ρn from the inside to the outside to fulfil the orientation convention.

5

 5

2 for p R and p L 10 :Figure 9 :

 109 Figure 9: Flowchart of CACOC algorithm (Alg. 6) where ρn refers to the value obtained from the Rössler system; it symbolises the chaotic dynamic.

Figure 10 :

 10 Figure10: Patterns obtained from the different mobility models with the same scale ratio. For MAMM, there are a lot of loops. The partition modification in MAMM2 permits to avoid loops but there are a lot of straight lines. For CROMM, the purple trajectory is a straight line followed by a serpentine pattern while the green line ended by a serpentine pattern followed by a large turn.

Figure 11 :

 11 Figure11: Example of metrics for each model for one single simulation computed to evaluate the performance over the 4000 first steps (we remind that we used 7000 steps for the experimentation). Linear regressions are used to evaluate the algorithms: the 500 first steps for the slope of coverage and the rest of the simulation (6500 steps) for the slope of fairness. The average value of the recent coverage is used for the whole simulation except the 100 first steps.

Figure 12 :

 12 Figure 12: Metrics comparison for the mobility models without ACO: Random, MAMM, MAMM2, CROMM. The smooth density estimation is given in addition to the box plot.

Figure 13 :

 13 Figure 13: Metrics comparison for the chaotic mobility models using ACO: ACO UAV and CACOC. The smooth density estimation is given in addition to the box plot.

Table 1 :

 1 Action table for the random mobility model.

			Probability of action
	Last action	Left	Ahead	Right
	Left Ahead Right	0.70 0.10 0	0.30 0.80 0.70	0 0.10 0.30

Table 4 :

 4 Main experimental parameters. Possible UAV actions ahead, 45 • left, 45 • right Initial UAVs position middle of the bottom of the map

	Parameter Name	Parameter Value
	Simulation area
	Geographical Area Number of cells UAV Autopilot	100m × 100m 100 × 100
	UAVs speed	1 m/s
	Experiments
	Algorithm Number of UAVs Simulation steps Independent runs	[Random, MAMM, MAMM2, CROMM, ACO UAV, CACOC] 10 7000 30

More information on ASIMUT is available on the project's website: https://asimut.gforge.uni.lu

Acknowledgments

The authors acknowledge the support of the ASIMUT project A-1341-RT-GP, which is coordinated by the European Defence Agency (EDA) and partially funded by 8 contributing Members (Austria, France, Germany, Italy, Luxembourg, The Netherlands, Poland and Sweden) in the framework of the Joint Investment Programme on Innovative Concepts and Emerging Technologies 2. The ASIMUT project consortium is composed of Thales, Fraunhofer IOSB, Fly-n-Sense, University of Bordeaux and University of Luxembourg (see https://asimut.gforge.uni.lu). The experiments presented in this paper were carried out using the HPC facilities of the University of Luxembourg [36] (see http://hpc.uni.lu).

$ The authors acknowledge the support of the ASIMUT project A-1341-RT-GP, which is coordinated by the European Defence Agency (EDA) and partially funded by 8 contributing Members (Austria, France, Germany, Italy, Luxembourg, The Netherlands, Poland and Sweden) in the framework of the Joint Investment Programme on Innovative Concepts and Emerging Technologies 2. The ASIMUT project consortium is composed of Thales, Fraunhofer IOSB, Fly-n-Sense, University of Bordeaux and University of Luxembourg.


F Fulchini

W Nan

M Ghadiri
email: m.ghadiri@leeds.ac.uk

M Yazdan Panah

S Bertholin

B Amblard

A Cloupet

T Gauthier

CFD-DEM Analysis of Particle Attrition in a Jet in a Fluidised Bed

In fluidised bed processes, the solids are in vigorous motion and thus inevitably subjected to mechanical stresses due to inter-particle and particle-wall impacts. These stresses lead to a gradual degradation of the particles by surface wear, abrasion and body fragmentation commonly termed attrition. One significant contribution of attrition comes from the air jets of the fluidised bed distributor. Particles are entrained into the air jet, where they get accelerated and impacted onto the fluidised bed particles. The jet induced attrition only affects the part of the bed which is limited by the jet length, where the mode of attrition is largely collisional. The overall jet attrition rate is therefore the result of the combination of the single particle damage and the flux of particles entering into that region. The attrition behaviour of particles in the jet region is analysed by evaluating their propensity of breakage experimentally and by simulating an air-jet in a bed of particles by CFD-DEM. The frequency of collisions and impact velocities are estimated from which the attrition due to a single air-jet is predicted.

Introduction

Circulating fluidised beds are widely used for carrying out heterogeneous reactions in which particulate solids are needed to be regenerated in order to be reprocessed. The circulation of the solids provides plenty of chances for attrition, especially where the solids meet the most stressful conditions, i.e. the gas distributors jets, bubbling bed and cyclones [1]. The jetting region is one of the most significant contributors to particle breakage. According to Zhang et al. [START_REF] Zhang | Proc.World Congr. Part. Technol[END_REF] and Boerefijn et al. [START_REF] Boerefijn | [END_REF], particles can undergo different types of stresses ranging from rapid shear to impact causing surface abrasion, chipping and fragmentation depending on the orifice to particle size ratio. The orifice velocity is the influential factor for accelerating the particles and projecting them upwards, where they collide with the bed particles leading to their attrition [4]. For orifice to particle size ratios greater than about 7.5 the main mechanism of jet attrition is collisional [START_REF] Zhang | Proc.World Congr. Part. Technol[END_REF][START_REF] Boerefijn | [END_REF]. Massimilla [5] developed a hydrojet dynamic model of the flow pattern in a fluidised bed jet from which particle velocity and flux can be calculated. Coupling of such information with a single particle breakage model can give an estimation of the overall jet induced attrition as proposed by Ghadiri et al. [6] using FCC catalyst. They first evaluated the attrition propensity by single particle impact testing and then coupled it with the hydrodynamic model of the jet of Donsi et al. [7] in order to analyse the dependence of the attrition rate on the orifice velocity and size. Two fundamental aspects are recognised to contribute to the jet attrition: (i) the jet hydrodynamics which dictate the particle velocity in the jet and the rate of entrainment of the particles into it, and (ii) the particle physical properties which influence the dependence of attrition on the operating conditions. This work aims to elucidate this phenomenon by evaluating the single particle attrition propensity of Chemical Looping Combustion (CLC) oxygen carrier particles coupled with an analysis of a single submerged jet by CFD-DEM simulation. CLC is a process developed on the circulating fluidised bed design. The particulate solids are looped between two fluidised beds, referred to as risers, in which the solids alternate chemical oxidation to chemical reduction. CLC attrition propensity is evaluated by the single particle impact test [8] and by Scirocco disperser impact test [9] for a wide range of particle impact velocities. The outcomes of the two tests are joined and a single particle attrition model is presented based on the chipping model of Ghadiri & Zhang [10]. The CFD-DEM simulation is used to calculate process and their entrainment rate in a jet. Their coupling enables the overall attrition rate to be calculated.

Materials and Method

The material, provided by IFPEN, Solaize, France consisted of a batch of fresh natural ores manganese based CLC oxygen carriers. The particle size distribution (PSD) is very wide and the particle shape very irregular as shown in Figure 1. Particle density is 3300 kg/m 3 .

Figure 1. CLC particle analysed by Scanning Electron Microscopy (SEM)

was evaluated gravimetrically by sieving using BS sieves (as given in Table 1), so that the most representative sieve cuts of the distribution could be chosen and used to carry on single particle impact test and Scirocco disperser impact test. Therefore, the effect of particle size on breakage could be taken into account in the breakage model. The two impact tests are able to correlate the impact velocity with the extent of breakage. The extent of breakage is evaluated by gravimetric analysis: the material is collected post impact and sieved using a sieve which is two sizes below the lower limit of the feed particle cut size in order to separate the mother particles from the fines, as described below. The extent of breakage is calculated as the ratio of the mass of fines to mass of feed particles. Ghadiri and Zhang model for chipping for semi-brittle materials [10] is applied to the extent of breakage in order to determine the lumped parameter representing material properties. In parallel, the CFD-DEM simulation of a single jet by EDEM TM and ANSYS-Fluent is carried out to evaluate the fluid-particle interactions in the jetting region. The flow of entrained particles into the jet, the number of collisions of each particle and corresponding velocity are calculated and used to estimate the overall extent of breakage.

Particle Size Distribution

The particle size distribution was evaluated by mechanical sieving, and reported in Figure 2. The sample was first split to get a small representative sample, suitable for sieving, and avoiding errors due to segregation of the initial large batch [START_REF] Germaine | Grain Size Analysis[END_REF].

Single Particle Impact test

The single particle impact test is according to the procedure described in [8]. Particle impact velocity is monitored and the extent of breakage is evaluated as:

* D D m m R = 100% m + m u ( 1
)
where mD is the mass of fines and mm is the mass of mother particles. The velocity range achievable for this test, for CLC, goes from free fall, about 1.8 m/s, up to 26 m/s depending on particle size.

Scirocco Disperser Impact Test

Scirocco is the dry disperser unit of Malvern Mastersizer 2000. Particles are fed and dispersed into an educator by a high pressure nozzle operating in the range from 0.5 to 4 barg. Given the L shaped bend of the educator, particles impact and experience attrition. Following passage through the disperser the particles are presented to the laser diffraction for the size measurement. They are recovered downstream the laser diffraction by a cyclone and analysed following the same procedure of the single particle impact test for consistency. A more detailed description is given by Ali et al. [9]. Varying the nozzle pressure from 0.5 to 4 barg, it is possible to achieve a range of impact velocity from 18 to 62 m/s depending on the particle size.

CFD-DEM Simulation of a Single Jet

The single jet simulation was carried out on a very simple geometry such as a cuboid of 100 mm×10 mm×600 mm mm size. The orifice is a 5 mm×10 mm rectangular slot and placed at centre bottom of the bed. Air, at normal condition, is issued from it at 40 m/s. The particles are considered spherical, having the diameter of 925 μm and the density of 3300 kg/m 3 . The other DEM parameters such as coefficient of restitution, sliding friction and rolling friction, shear modulus and Poisson's ratio were specified as, 0.5, 0.5, 0.01, 0.1 GPa and 0.25, respectively, as given by default by the softwere. The total number of particles was set as 20000 for a total bed mass of 31.97 g.

The system was simulated for 3.4 s.

Results and Discussion

Single Particle Impact Test and Scirocco Disperser Impact Test

Results from the impact tests, for each size cuts, are reported below in Figure 3. The extent of breakage is plotted as a function of particle impact velocity. The extent of breakage is corrected to take account of the loss of debris and mother particles. A very small amount of Powders & Grains 2017 debris is lost in the Scirocco because the collecting efficiency of the cyclone is not 100%, so the extent of breakage is referred to as R + . Whilst in case of the single particle impact test, the operation of manual sieving causes some losses of mother particles. In this latter case the extent of breakage is R -. Given that the same gravimetric analysis performed for the two tests, the results are combined for the full range of particle velocities, as shown in Figure 3. The trend is parabolic, indicating that breakage is dependent on the incident impact energy. Besides, larger particles break more than the smaller ones which is consistent with the model of Ghadiri and Zhang [10], addressed below. Applying this model to the experimental data points a unification is achieved.

2 2 p p c H R = αη = α ρd V K (2)
where ρ is the particle density, dp is the particle size and Vp the particle velocity. The logarithm of the extent of breakage, log(R), is plotted as a function of log(ρdpVp 2 ) in Figure 4. The slope of the straight line is a lumped parameter, representing the mechanical properties of the particle, αH/Kc 2 , where α is a proportionality constant, H is the hardness and Kc is the fracture toughness. The chipping model describes very well the pattern of impact breakage. Interestingly, considering the intercept of the fitted line with abscissa, there is a minimum impact velocity for a given particle size below which there is little/no breakage. The intercept can be used to estimate the transition velocity from no breakage to chipping for any particle size. In our case, considering particles of 925 μm, the transient velocity is 0.4 m/s.

CFD-DEM Simulation of a Single Jet in a Fluidised Bed

The dynamics of particle motion in the jet region was analysed in terms of particle velocity, number of collisions per particle and flux of particles into the acceleration zone. The breakage model developed in the experimental part is applied in order to estimate the overall attrition in the jet region. Particles enter the channel, created by the air jet, from the sides of the bed, Figure 5.

According to the review of Boerefijn et al. [START_REF] Boerefijn | [END_REF], in the case of large orifice to particle size ratios (5.4 in our case), particles entering the high velocity channel accelerate till they collide with the particles above the jet, experiencing attrition. The acceleration zone is limited to the channel height (about z=20 mm) as it is clear from Figure 6. Therefore, the whole analysis of attrition is limited to the jet core, where particles experience attrition by collision after they are accelerated. The breakage model of Ghadiri and Zhang was applied taking account of the flux of

where R, fc and Q are calculated every 0.01 s within the jet core, using the average impact velocity (Vp), the number of particles (Np), the number of collisions (Nc), the mass of single particle (mp) and Δt equal to 0.01 s. In this condition, the breakage rate, i.e. the mass of debris formed by the jet per unit time, is calculated to be 0.025 g/s, which results in 0.084 g of mass of attrition product for the complete simulation time. At the end of the simulation, the bed inventory has therefore lost 0.26% of its mass due to the attrition of a single jet.

p p N m Q t ' (4)
c c p N f N (4)
From the mass of fines produced, the debris size distribution of fines is calculated as follows, assuming that each impact generates only one fine particle.

jet D B d = 3 p 1/3 c d Q f (5)
The calculated particle size distribution after attrition is thus shown in Figure 7.

Conclusions

The attrition propensity of CLC particles is evaluated experimentally using two impact tests: the single particle impact test and Scirocco disperser impact test. Using the model of chipping for semi brittle material of Ghadiri and Zhang, a relationship between the extent of breakage and impact velocity is established, taking into account particle size and mechanical properties. The extent of breakage follows the square of the impact velocity and varies linearly with particle size. Particle mechanical properties such as hardness and fracture toughness were determined as a lumped parameter by plotting the extent of breakage as a function of particle density, size and velocity squared. A CFD-DEM simulation of a single jet in a fluidised bed is then performed and used to evaluate the overall attrition using the single particle breakage model and coupling it with the number of collisions per particle within the jet region and the flow of particles entering the jet core. The analysis reveals that, under these conditions, 0.025 g/s of fines is produced by attrition and that the average debris size is about 60 μm. This analysis is a starting point for a more detailed and realistic approach for predicting the attrition rate, since the single particle attrition propensity can be established for any system at any conditions and coupled with CFD/DEM to estimate the attrition rate. Validation of the simulation is planned next in order to finalise the model for a single jet. Powders & Grains 2017

Figure 2 .

 2 Figure 2. PSD of CLC oxygen carriers

Figure 4 .Figure 5 .Figure 3 .DOIFigure 6 .

 4536 Figure 4. Extent of breakage: log(R + ) for Scirocco disperser impact test and log(R -) for Single Particle Impact test, as a function of log(dpVp ) showing data unification

Figure 7 .

 7 Figure 7. PSD before and after attrition

Table 1 .

 1 Representative size cuts used and sieve size for separating mother particles/fines for the breakage analysis Representative

size cuts (μm)

		180-	212-	250-	300-	355-
		212	250	280	355	400
	Cut off sieve size (μm)	125	150	180	212	250

Table 2 .

 2 Time averaged values for Q, fc and Vp

	Q [g/s]	fc [-]	Vp [m/s]
	35.39	5.30	0.83

Acknowledgements

The financial support from Total, France, is gratefully acknowledged.