José M Mendes

Eric Filho

David Filliat Lucet

Filliat

Experimental Validation of a Multirobot Distributed Receding Horizon Motion Planning Approach

come

Experimental Validation of a Multirobot Distributed Receding Horizon Motion Planning Approach

José M. Mendes Filho a, b, * , Eric Lucet a and David Filliat b

Abstract-This paper addresses the problem of motion planning for a multirobot system in a partially known environment where conditions such as uncertainty about robots' positions and communication delays are real. In particular, we detail the use of a Distributed Receding Horizon Approach that guarantees collision avoidance with static obstacles and between robots communicating with each other. Underlying optimization problems are solved by using a Sequential Least Squares Programming algorithm. Experiments with real nonholonomic mobile platforms are performed. The proposed framework is compared with the Dynamic Window approach to motion planning in a single robot setup. A second experiment shows results for a multirobot case using two robots where collision is avoided even in presence of significant localization uncertainties.

I. INTRODUCTION

A Multirobot System (MRS) is a system composed by multiple robots operating in the same environment. These systems have received a great deal of interest for their potential to accomplish complicated tasks in different application domains such as exploration [START_REF] Burgard | Coordinated multi-robot exploration[END_REF], [START_REF] Nestmeyer | Decentralized simultaneous multi-target exploration using a connected network of multiple robots[END_REF], transportation [START_REF] Maza | Multi-uav cooperation and control for load transportation and deployment[END_REF], [START_REF] Amato | Planning for decentralized control of multiple robots under uncertainty[END_REF], construction [START_REF] Augugliaro | The flight assembled architecture installation: Cooperative construction with flying machines[END_REF] and more.

An indispensable capability for MRS is Motion Planning (MP). MP can be broadly defined as the problem of generating a collision-free trajectory (sequence of poses and velocities) from an initial to a final pose. Although simple in its definition, this problem is PSPACE-hard [6]. For the past decades, many approximate approaches to solve this problem have been proposed [START_REF] Choset | Principles of Robot Motion: Theory, Algorithms, and Implementations[END_REF], [START_REF] Lavalle | Planning Algorithms[END_REF]. It remains, though, a major focus of robotics, especially when robots' workspaces are dynamic partially known environments -in which case the problem is shown to be intractable (NP-hard) [START_REF] Canny | New lower bound techniques for robot motion planning problems[END_REF]. Furthermore, far less work has been carried out for addressing the MP problem for MRS than for single robot systems, as pointed out in [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF].

An increasingly common approach for addressing this problem is the use of mathematical programming. It offers flexibility to explicitly accommodate multiple systems requirements simultaneously. In most cases, these requirements are a subset of the following list: kinematics, dynamics, collision avoidance and connectivity maintenance requirements [START_REF] Abichandani | Mathematical programming for multi-vehicle motion planning problems[END_REF].

Our previous work presented in [START_REF] Mendes Filho | Real-Time Distributed Receding Horizon Motion Planning and Control for Mobile Multi-Robot Dynamic Systems[END_REF] fits this group of mathematical programming methods by proposing a formulation of the MRS MP problem as distributed nonlinear a CEA, LIST, Interactive Robotics Laboratory, Gif-sur-Yvette, F-91191, France b U2IS, ENSTA ParisTech, Inria FLOWERS team, Université Paris-Saclay, 828 bd des Maréchaux, 91762 Palaiseau cedex France * corresponding author: jose.mendesfilho@cea.fr programming problems (NLPs) solved by each robot in the system. It then interleaves planning and execution in the partially known environment by using a Receding Horizon approach. It aims to generate dynamically feasible and collision-free trajectories for a MRS composed of nonholonomic ground vehicles modeled as unicycle-type robots. NLPs are numerically solved by using a Sequential Least Squares Programming algorithm (SLSQP) [START_REF] Kraft | A software package for sequential quadratic programming[END_REF], [START_REF] Johnson | The NLopt nonlinear-optimization package[END_REF]. The complete approach is referred to as DRHMP (Distributed Receding Horizon Motion Planning).

This paper builds on that work with three main objectives: i) to discuss strategies and implementation techniques in order to use the DRHMP approach in real mobile robots; ii) to compare the DRHMP to another MP approach for a single robot system; iii) to demonstrate the robustness of the DRHMP in presence of localization and tracking errors and real communication conditions between two robots.

In order to accomplish those objectives, the DRHMP approach was implemented on two Turtlebots 2 platforms [START_REF]TurtleBot2 -Open-source robot development kit for apps on wheels[END_REF] equipped with rangefinder sensors (structuredlight 3D scanners). Its performance was compared with that of a Dynamic Window approach (DWA) used by default in the ROS navigation stack.

The paper is thus organized as follows: the next section gives a brief summary of related work on MRS MP, Section III provides a presentation of the DRHMP approach. In Section IV the experimental setup for verifying that approach is detailed along with the most significant contributions of this work. The last two sections present results and conclusions.

II. RELATED WORK

Approaches for solving the MRS MP problem can be classified in many different ways. Refer to [START_REF] Hoy | Algorithms for collisionfree navigation of mobile robots in complex cluttered environments: A survey[END_REF], [START_REF] Mohanan | A survey of robotic motion planning in dynamic environments[END_REF] for comprehensive surveys on the subject. Here we aim only to review a few recent works that closely relate to ours.

In [START_REF] Zhou | A Real-Time and Fully Distributed Approach to Motion Planning for Multirobot Systems[END_REF] a mathematical programming and distributed receding horizon approach is used for MRS MP. Robots are required to communicate and exchange their current and most recent states. That information is used to predict robots trajectories assuming uniform linear motion. Then, those predictions are used to form collision avoidance constraints in the optimization problem. An interesting aspect of this approach is that it splits the planning horizon into two parts: during the first part collision avoidance and smoothness of trajectories are dealt with; in the second part only global target convergence is a concern. An incremental sequential convex programming (iSCP) algorithm for solving the optimization problems is used. Only results in simulation are shown.

Work presented in [START_REF] Filotheou | Robust decentralized navigation of multi-agent systems with collision avoidance and connectivity maintenance using model predictive controllers[END_REF] proposes a Decentralized Nonlinear Model Predictive Control (DNMPC) for addressing MRS MP where careful convergence and feasibility analysis are provided. Formation maintenance, avoidance of static obstacles and inter-robot collision avoidance are verified in simulation for unicycle robots but the approach could be generalized for other types of systems. The method requires though that the robots communicate sequentially and it is not clear how the underlying finite horizon optimization control problems for each robot are solved. This approach too remains to be tested in a real experiment.

Another recent relevant work on Distributed MPC for MRS is presented in [START_REF] Mehrez | Occupancy grid based distributed mpc for mobile robots[END_REF]. Instead of relying on complete predictions of other robots trajectories it uses occupancy grid data aiming for a reduction in the required communication means. The approach was tested on nonholonomic mobile platforms using an external motion capture system for localization of the mobile platforms.

Our approach closely relates to the one presented in [START_REF] Zhou | A Real-Time and Fully Distributed Approach to Motion Planning for Multirobot Systems[END_REF]. The main differences consist in how collision constraints are handled and how localization and tracking errors are modeled. In [START_REF] Zhou | A Real-Time and Fully Distributed Approach to Motion Planning for Multirobot Systems[END_REF] the problem of having non-differentiable constraints for obstacle avoidance is addressed by transforming them into smooth nonlinear constraints. Conversely, our work derives differentiable smooth constraints from sampled data (inflated occupancy grids given by the perception module) by doing local interpolations around sampled points in robots' planned trajectories. As for localization and tracking errors, [START_REF] Zhou | A Real-Time and Fully Distributed Approach to Motion Planning for Multirobot Systems[END_REF] assumes they are always inferior to a small constant while our work uses a probabilistic model and confidence regions to produce robust collision-free trajectories.

III. DISTRIBUTED RECEDING HORIZON MOTION PLANING

In the DRHMP approach, each robot in the MRS computes its own local trajectory. Analogous to a MPC, a prediction time-horizon T p and a computation timeslot T c are defined 1 .

T p is the time-horizon for which a local solution to the motion problem will be computed and T c is the timeslot during which a portion of that solution is implemented while the next plan -created for the next time-horizon T p -is being computed.

For each receding horizon planning problem, the following two steps are performed:

Step 1: Each robot computes its own intended solution trajectory by solving a constrained nonlinear optimization problem (NLP). In that problem all needed constraints are 1 The term receding horizon is favored over the more used term MPC [START_REF] Grne | Nonlinear model predictive control: theory and algorithms[END_REF] because the DRHMP uses more than just the first element of the optimal solution for computing the system's input. In other words, since T c is usually greater than the lower level controller time period T t , the first T c /T t elements of the solution are implemented.

included except coupling constraints, meaning constraints that involve solving a conflict between multiple robots such as collision or loss of communication. The solution found at this step is referred to as step1-generated trajectory.

Step 2: Robots involved in a potential conflict (risk of collision, or communication loss) update their trajectories computed at Step 1 by solving a second NLP that additionally takes coupling constraints into account. For evaluating those additional coupling constraints the robots asynchronously exchange their step1-generated trajectories. If a robot is not involved in any conflict, its step2-generated trajectory becomes equal to the one found at Step 1 avoiding solving the second NLP.

When a robot arrives closer to its goal, a termination procedure for reaching the goal is used. It takes the desired final configuration as a hard constraint in the optimization problem and uses the time for reaching the goal as the cost function to be minimized. Refer to works [START_REF] Mendes Filho | Real-Time Distributed Receding Horizon Motion Planning and Control for Mobile Multi-Robot Dynamic Systems[END_REF], [START_REF] Mendes Filho | Multi-Robot Motion Planning: a Modified Receding Horizon Approach for Reaching Goal States[END_REF] for more details.

An example of the resulting receding horizon optimization problem solved at Step 2 for a given robot R can be written as follows: min

q(t) q(τ k + T p ) -q goal , ∀t ∈ [τ k , τ k + T p ] (1)
subject to:

                               q(t) = f (q(t), u(t)), ∀t ∈ [τ k , τ k + T p ] q(τ k ) = q(τ k-1 + T c ) q(τ k ) = q(τ k-1 + T c ) v 2 (t) ≤ v 2 max , ∀t ∈ (τ k , τ k + T p ] ω 2 (t) ≤ ω 2 max , ∀t ∈ (τ k , τ k + T p ] a 2 (t) ≤ a 2 max , ∀t ∈ [τ k , τ k + T p ] α 2 (t) ≤ α 2 max , ∀t ∈ [τ k , τ k + T p ] d(O,t) ≥ ε o , ∀t ∈ (τ k , τ k + T p ], ∀O ∈ O d(R c ,t) ≥ ε r (R c ), ∀t ∈ (τ k , τ k + T p ], ∀R c ∈ C d(R d ,t) ≤ min(d com , d d,com ), ∀t ∈ (τ k , τ k + T p ], ∀R d ∈ D
with τ k the starting discrete time for the kth receding horizon optimization problem, f the kinematic model of the robot, q = (x, y, ψ) the vehicle's pose, In our previous work, when applying this approach, several assumptions about communication, perception and magnitude of different errors were made that do not hold in a real application.

u the system's input, v = [ ẋ ẏ] , ω = ψ, a = [ ẍ ÿ] , α = ψ,
First, when testing in simulation, computations of each robot were made by a different thread under a unique process. This allowed the use of shared memory for communicating among virtual robots which could be considered instantaneous within that context.

Another unrealistic assumption was about perception. Once an obstacle entered the perception range of a vehicle, a perfect knowledge of that obstacle's shape (polygon or circle) and pose was considered. It was possible then to solve simple algebraic expressions for computing the distance from a given point to the obstacle.

Finally, localization of simulated robots was precise and accurate and tracking errors were reasonably small. That allowed for the use of constant, small tolerance values for constraints concerning collision avoidance (i.e. constant, small ε o and ε r ).

IV. EXPERIMENT SETUP

This section details how each of the unrealistic assumption presented above -concerning communication, perception and errors -were compensated for in order to use the DRHMP in a real experiment. Two Turtlebots 2 [START_REF]TurtleBot2 -Open-source robot development kit for apps on wheels[END_REF] mobile robots were used for such purpose.

A. MRS architecture

The DRHMP approach was implemented as a plugin (named codrha local planner) adhering to the BaseLocalPlanner interface provided by the nav core ROS package. This allows for using other features already present in the ROS navigation stack [START_REF]ROS 2D navigation stack[END_REF] such as occupancy grids, global planners, etc.

Each robot launched its own ROS master node and all other nodes needed for navigation, including an instance of codrha local planner.

Communication between robots was conducted outside ROS systems for messaging (typically Topics or Services). Specific socket programming was implemented so exchange of information among robots could be as minimal as possible.

A mobile ad-hoc network (MANET) over Wi-fi (IEEE 802.11) was used. Ad-hoc networks represent a more robust architecture compared to infrastructure mode networks (which use central Access Point) and more flexible, scalable and cheap than architectures with fixed topology in general [START_REF] Giordano | Mobile ad hoc networks[END_REF].

The embedded computers on each robot had a dual core Intel Celeron N3060 @ 1.60GHz CPU. For better performance, priority and affinity configurations of threads had to be set properly. Communication-related threads were set with affinity to one of the cores (e.g. core 0) with higher priority than any other processes in the system (e.g. priority 10). In turn, all ROS-related threads were launched with affinity to the remaining core (e.g. core 1) with priority only lower than communication-related threads (e.g. priority 9).

Time synchronization among robots is a critical aspect of any MRS. In the DRHMP the need of a common global time is fundamental to the evaluation of conflict between vehicles. In a wireless network, the problem becomes even more challenging due to the possibility of collision of the synchronization packets on the wireless medium and the higher drift rate of clocks on the low-cost wireless devices.

In this experiment, time synchronization was achieved by using an alternative NTP client and server called chrony [START_REF]Chrony -an alternative NTP client and server[END_REF], designed for systems that are not permanently online. It is supposed to provide an accuracy typically in the tens of microseconds.

An architecture representing the most relevant components of a robot in the MRS can be seen in Fig. 1.

B. Localization and tracking error compensation

Localization for vehicles used in the experiments was done by using a particle filter [START_REF] Thrun | Probabilistic robotics[END_REF] merging data from gyroscopes, encoders, rangefinder sensors and static map. As well as for any localization system, there is inherent error associated with the robot's estimated configuration. That uncertainty can be characterized by the covariance matrix of the robot's configuration which is known for each robot at each instant τ k . Likewise, each robot can estimate its tracking error based on the planned reference trajectory and its configuration estimate.

Both covariation values and tracking error were sent as part of the information exchanged between robots along with their estimated configuration and planned trajectory. This enables a given robot R to compute a conservative "safety distance" from robot R c that is time dependent and can replace the constant ε r (R c ) in the NLP constraints.

The new value is computed according to the equation:

ε r (R c ,t) = ε tr (t) + ε loc (t) + ε tr, R c (t) + ε loc, R c (t) (2)
where ε tr is simply the euclidean norm of lateral and longitudinal tracking errors and ε loc is a confidence value related to position estimates. ε loc depends on the covariance matrix given by the robots' localization modules, and its expression is given by Eq. 3:

ε loc (t) = χ 2 v (2) λ 0 , λ 1 ∞ (3)
with λ 0 and λ 1 being the eigenvalues of the covariance matrix associated with the bivariate normal distribution of possible locations in the XY plane of a robot at an instant τ k . χ2 v (2) is the 2 degrees of freedom chi-square value for a (1v) confidence region (CR). 2 Similarly, ε o in the NLP constraints can be replaced by ε o (t), where:

ε o (t) = ε tr (t) + ε loc (t) (4)
Fig. 2 shows a graphical representation of errors ε tr (t) and ε loc (t) at instant τ k .

C. Use of costmap for obstacle avoidance

Instead of a complete geometric description of obstacles as used in previous works [START_REF] Mendes Filho | Real-Time Distributed Receding Horizon Motion Planning and Control for Mobile Multi-Robot Dynamic Systems[END_REF], [START_REF] Mendes Filho | Multi-Robot Motion Planning: a Modified Receding Horizon Approach for Reaching Goal States[END_REF], only 2D costmaps are available for the purpose of obstacle avoidance when using ROS navigation stack. These costmaps are based on occupancy grids generated from 3D sensor data and user defined inflation radius that is meant to take the robot's footprint into account.

Directly using those costmaps for computing obstacle avoidance constraints in the DRHMP is not straightforward. The solution that worked best can be described in two steps.

First, the occupancy grid is inflated according to a linear function that goes from the highest possible cost value at an occupied position in the grid, to zero at a position located at a distance equals to the radius of the robot.

Secondly, bicubic interpolations of the costmap around sample points taken along a trajectory candidate are performed. Those interpolations provide a locally defined, continuous and differentiable distance function for each sample during optimization. A representation of this approach for one sample is displayed in the Fig. 3.

The importance of having a continuous differentiable distance functions comes from the way the SLSQP solver searches for a solution: it uses finite differences to estimate first and second order derivatives of constraints. Without interpolations, those finite differences would usually evaluate to zero. That is because the differentiation step size h employed by SLSQP is much smaller than the occupancy grid resolution.

V. RESULTS

Experiments were carried out in order to investigate two aspects: how the DRHMP compares to another local MP in a "single robot avoiding an obstacle" situation and how well collision avoidance between two robots running the DRHMP is performed in face of real communication, perception and trajectory tracking issues. Across all experiments we used v max = 0.2 m/s, a max = 0.5 m/s 2 , T p = 3.8 s, T c = 0.3 s. A simple controller designed for tracking a reference admissible trajectory based on the kinematic model of the system was used [START_REF] Siciliano | Springer handbook of robotics[END_REF].

A. Experiment 1 -Single robot obstacle avoidance

A testbed as shown in Fig. 4 was used for comparing collision avoidance with a static obstacle using the DRHMP approach and the well known Dynamic Window approach (DWA) [START_REF] Fox | The dynamic window approach to collision avoidance[END_REF]. Although admittedly each planners' performance can be highly impacted by the configuration of its parameters, an effort was made to set those values so equivalent behaviors could be obtained. Velocity and acceleration limitations were set to same values for instance. As indicated in Fig. 4 the "robot" located near the XY origin had to reach the "target" point located about 2.3 m of its initial position. An obstacle of about 20 cm in diameter not known by the robot in advance was placed in the way.

Table I summarizes the performance of each planner regarding five criteria. Mean and standard deviation (Std) data were based on 10 trials for each algorithm. In none of the tests neither MP approaches failed to avoid obstacles. Compared to DWA, DRHMP can be seen as a less conservative approach. It keeps less clearance from obstacles and can produce a higher average linear velocity in order to minimize travel time. That behavior derives from the type of objective function used in the NLPs. Similarly, DWA's behavior derives from its scoring algorithm. It worth noticing that DRHMP presents an inferior standard deviation compared to DWA which suggests it is a more stable approach. Typical paths adopted by both algorithms can be seen in Fig. 4a.

B. Experiment 2 -Multirobot motion planning

Experiment 2 consisted of having two robots going successively from one target location in their shared workspace to another. Those targets positions were such that the robots would execute two different triangle-shaped loops that share a common side. Along this shared side the two robots (if the timing was right) would have to cross each other to reach their next target. The robots' trajectories in Fig. 5 illustrate this setup.

For better evaluating the DRHMP performance in that scenario two different sub-experiments were set. In sub-experiment 2a a simplification about tracking and localization errors was made. They were considered equal to zero by the DRHMP algorithm running in both robots (ε r = 0). It implied though that the physical robots had to be kept still during planning to prevent them from colliding. As we will see, the ε r = 0 assumption is far from realistic, provided the platforms we worked with. Sub-experiment 2b, on the contrary, makes no such simplification and ε r is based on the real information about the physical robots executing the planned trajectories. This second case shows how even with errors of about 50 cm the DRHMP can safely find collision-free trajectories for both physical robots.

Everything else, specially communication is done in the same way across both experiments.

1) Sub-experiment 2a:

The trajectories produced and the inter-robot distance along this entire sub-experiment can be seen in Figs. 5a to 5e. Planned trajectories would allow both robots to avoid each other with almost no clearance at four different moments as indicated by the inter-robot distance curve passing near zero in Fig. 5a.

Figs. 5b to 5e represent four snapshots of both robots planning processes around the first moment of collision avoidance. The green continuous lines in front of the robots represent the XY points of their planned trajectory for the horizon T p (i.e. step2-generated trajectories). The dashed line in front of a given robot R a represents what the other robot R b thinks R a will do (i.e. R a 's step1-generated trajectory as known by R b as a result of their communication).

When collision is not an issue, dashed lines may superpose almost exactly the green continuous ones 3 . In contrast, as collision becomes an issue, one can observe a greater 3 communication delays may still prevent them from being identical Following the time sequence and observing Fig. 5d, one can see that the actually followed paths after solving the conflict reflect a smaller deviation from the otherwise straight-line trajectory when compared to those planned green lines in Fig. 5c. This is due to the Receding Horizon nature of the approach which interleaves planning and execution and therefore reconsiders the avoidance trajectory adapting to the changes made by the other robot. Its implication, albeit conditioned to T c values, is that collision avoidance ends up being achieved with near the minimum clearance possible (zero in this case). Furthermore, a natural compromise between robots is achieved: they both deform their initial straight trajectories of comparable amounts.

That is precisely the behavior expected for the DRHMP in a multirobot scenario and observed in simulation in previous works [START_REF] Mendes Filho | Real-Time Distributed Receding Horizon Motion Planning and Control for Mobile Multi-Robot Dynamic Systems[END_REF].

2) Sub-experiment 2b:

In this sub-experiment the physical robots actually did navigate their workspace and the real observed information about their localization and trajectory tracking was used in the DRHMP.

Due to imprecision of sensors, actuators and low level controllers, ε r becomes meaningful in this second case as shown by the error bands in Fig. 5f. That has the effect of reducing the space of acceptable solutions in which the DRHMP searches for optimal trajectories. If, additionally, the workspace is very cluttered there may be no acceptable solutions left yielding optimization errors at the SLSQP algorithm level. Nevertheless, for the studied setup, acceptable solutions were always found and the inter-robot collision of the physical robots were prevented at all times throughout the experiment (as observed in the experiment video in [START_REF] Mendes Filho | Distributed Receding Horizon Motion Planning for Multirobot Systems[END_REF]).

From Fig. 5g to 5j the difference between black and colored circles representing the robots reflects the tracking error. Black circles use the mean information of robots' localization systems while colored ones use the planned poses. As before, dashed and green continuous lines in front of a given robot represent respectively the step1-generated trajectory known by the other robot and its own step2generated trajectory.

The utility of a nonconstant ε r taking robots localization uncertainty into account can be appreciated in Fig. 5f. It can be observed that when ε r takes smaller values (such as near time 45s compared to time 20s) the inter-robot distance can be reduced as well.

Table II summarize statistics of experiment 2 regarding communication and inter-robot distances. Comparing communication-related values between robots R a and R b and then between sub-experiments 2a and 2b shows indeed that communication conditions where very symmetric. The column titled "Obs" shows the number of observations and it is roughly equal to the duration of the experiment divided by T c (the period the robots exchange step1-generate trajectories). Furthermore, due to the asynchronism between robots' planning processes and to communication delays, it is common that part of the trajectory information received by a robot concerns a time interval of no interest to that robot. In other words, at Step 2 of the DRHMP the information about another robot's trajectory may be partially too old or planned for too far into the future. The percentage of the information that can actually be used is referred to as URI.

Overall, despite considerable communication delays (tens of ms), low quality localization (main responsible for high ε r values) and use of a simple kinematic controller that does not take sliding and actuator response times into account, the DRHMP manages to produce satisfying results with respect to multirobot collision avoidance.

VI. CONCLUSIONS

In this paper, a Distributed Receding Horizon Motion Planning approach for Multirobot Systems is presented and tested in real experiments. We detail its implementation on two mobile platforms as a mean of discussing solutions for dealing with real application limitations concerning communication and perception. Its performance is compared with the Dynamic Window approach in a single robot setup. Tests for multirobot system are also carried out showing the robustness of the DRHMP in presence of localization uncertainties and real communication delays.

Among the future research directions that could be explored based on this work, we particularly think that a generalization of the approach for different systems (not only unicycle-like vehicles) for supporting non-homogeneous multi-robot systems is of great interest and relatively straightforward.

Moreover, in this work we assume that robots can obtain information about other robots by directly communicating with them (as in a fully connected topology). This is an unrealistic assumption for robots in a mobile ad-hoc network where topology is dynamic and self-organizing. The use of an efficient and secure distributed ledger [START_REF] Ferrer | The blockchain: a new framework for robotic swarm systems[END_REF] for sharing information among all robots in the network would be an interesting improvement and enable tests with a greater number of robots without prohibitively increasing communication requirements.

  O the set of obstacles known by robot R at τ k-1 , C the set of collision candidates vehicles for robot R at τ k-1 , D the set of communication loss candidates for robot R at τ k-1 , d com the maximum communication range of robot R, d d,com the maximum communication range of robot R d , ε o the minimum distance allowed between robot R and any obstacle, and ε r (R c ) the minimum distance allowed between robots R and R c .

Fig. 1 :

 1 Fig. 1: Overview of a robot composing the MRS

Fig. 2 :

 2 Fig. 2: Localization and tracking errors representation. The 95% Confidence Region ellipse encompasses the central 95% of the probability mass of possible locations for a robot at instant τ k . The ellipse is dilated by the robot's radius.

Fig. 3 :

 3 Fig. 3: 3D representation of the costmap and local interpolation used by the codrha local planner. Cost values of 254 correspond to the detected surface of an obstacle

Fig. 4 :

 4 Fig. 4: Experiment 1. (a) shows a representation of the testbed using the same map information as the localization system of the real robots. (b) is an actual photo of the testbed

Fig. 5 :

 5 Fig.5: Experiment 2 results. (a) to (e) refer to sub-experiment 2a while (f) to (j) to sub-experiment 2b. "Planned" or "P" data shows inter-root distance computed based on the planned trajectories generated by the DRHMP approach. "Observed" or "O" data represents the same distance but based on observed/estimated actual position of the robots. In (f) the blue and purple error bands centered around the "Planned" line represent respectively the safety distance of Eq. 2 (ε r ) and its localization component (ε loc,R a + ε loc,R b ).

TABLE I :

 I Experiment 1 summary

		DRHMP	DWA
		Mean	Std	Mean	Std
	travel time (s)	12.637 0.084 18.790 0.385
	average linear speed (m/s)	0.193 0.000	0.136 0.002
	final position error (m)	0.014 0.004	0.092 0.013
	final yaw error (rad)	0.011 0.004	0.355 0.053
	min clearance from obstacles (m)	0.069 0.014	0.228 0.023

TABLE II :

 II Experiment 2 summary 's step1-generated trajectory as known by R b is shown as coming straight at R b . R a 's step2-generated trajectory, in green, is quite different from that and clearly avoids a future collision. That is because R a has already taken into account R b 's step1-generated trajectory into its Step 2.

			Min	Max	Mean	Std	Obs
	Planned IRD	0.014	-	-	-	-
		COD	0.000	0.403	0.016 0.044	370
	R a	USD	0.282	1.187	0.472 0.158	370
	2a	URI	64.506 92.586 86.920 5.249	370
		COD	0.000	0.265	0.017 0.038	323
	R b	USD	0.280	0.895	0.512 0.134	323
		URI	64.204 92.623 85.962 4.714	323
	Planned IRD	0.376	-	-	-	-
	Observed IRD	0.253	-	-	-	-
	ε r (m)	0.292	0.933	0.435 0.106	351
		COD	0.004	0.232	0.014 0.026	355
	2b R a	USD	0.073	0.973	0.408 0.103	355
		URI	61.112 98.089 88.865 3.608	355
		COD	0.000	0.186	0.018 0.036	351
	R b	USD	0.201	0.821	0.526 0.064	351
		URI	68.930 94.701 85.660 3.223	351
	RI = received information
	URI = % of the RI that is actually used by the DRHMP at Step 2
	COD = communication delay measured by the receiving end in seconds
	USD = delay between receiving and start using RI in DRHMP in seconds
	IRD = inter-robot distance in meters
	difference between those two types of trajectories. At Fig. 5c,
	R a

A 95% CR (i.e. χ 2 0.05 (2) = 5.991) was used in the experiments.


Gabriel Vasile
email: gabriel.vasile@gipsa-lab.grenoble-inp.fr

Pascal Bellemain

Eric Guy D'urso
email: guy.durso@edf.fr

De Oliveira
email: eric.de-oliveira@edf.fr

Real-Time Wideband Acoustic Signal Processing for Fish Counting

This paper proposes an IoT active acoustic system for fish monitoring with direct application in biodiversity preservation near hydropower plants and dams. The system employs a cost effective single board computer as central processing unit. The acoustic signal processing is performed at two levels: the first at the board level (embedded), to decide upon transmission of data to a server, and the server level to show the results of detection and the number of fish passages (acoustic signature analysis). The performance analysis of the proposed fish monitoring system is carried out during a four days experimental campaign in Chatou, France.

I. INTRODUCTION

Dam and hydropower plants are currently known to obstruct migration of fish species along the watercourses during certain phases of their lifetime. However, the 2030 Climate and Energy Framework, agreed in 2014, imposed to reduce carbon footprint by 40% [START_REF] Parlament | Proposal for a DIRECTIVE OF THE EUROPEAN PARLIAMENT AND OF THE COUNCIL on the promotion of the use of energy from renewable sources (recast)[END_REF], which makes this renewal energy more and more used. We are in the situation of trading the preservation of fish biodiversity against the production of clean energy and, hence, there is a stringent search for alternative measures to mitigate the impact on biodiversity. A popular approach is to construct channels where the fish by-pass the targeted hydraulic installations.

In this work, we present the signal processing algorithms embedded in our prototype IoT (Internet of Things) acoustic system for fish activity monitoring, with the objective to count the fish getting through any by-pass channel. In a separate paper proposed to the same conference, we have already presented the hardware configuration of the system. Illustrated in the Fig. 1, it contains a FPGA based single board computer for both signal acquisition and arbitrary function generation (RedPitaya STEMLab 125-14, RP), 2 transmitting ultrasound transducers (Ultran Group WS-100-1), 2 ultrasound transducers in reception (KS-100-1), a power amplifier (APEX Microtechnology PA107), a low noise amplifier(TI VCA 5807) and a server for real-time signal processing (Intel NUC). The sensors are calibrated by sphere-calibration method as described in [START_REF] Demer | Calibration of acoustic instruments[END_REF].

The paper is structured as follows. Section II illustrates in several steps the proposed software architecture of underwater fish monitoring system, while Section III presents both qual- itative and quantitative performance assessment. Section IV concludes the paper.

II. IOT SOFTWARE ARCHITECTURE

The proposed signal processing algorithm operates in realtime and it is composed of two stages:

• the IoT client -an embedded system capable of managing signal acquisition, real time processing (intelligent triggering option) and sending the recorded data to the server ; • the IoT server -a processing unit for data management, storage and final signal processing. The former is considered in taking decisions to transmit the acoustic signature (or not) to the server (energy saving) for further analysis, whereas the latter is used for detecting, recognising and counting the fish passages.

A. IoT client

The client software has been developed in C using the POSIX threads / sockets and the RedPitaya libraries on the embedded Ubuntu operating system. The software is continuously configuring and controlling emission, reception and data transmission in real time using parallel threading on the embedded coprocessor.

The RP board's arbitrary function generator is programmed to transmit Hamming weighted linear frequency modulation signals (chirps) [START_REF] Vasile | Transmission waveform analysis for active ultrasound fish monitoring systems[END_REF] in the band comprised between 350 kHz and 650 kHz. Each reception has 2 × 16384 samples stored on a single "signed short" vector. In order to achieve better transmission results, the last N = 100 recorded signals are locally stored in a circular buffer. Once this buffer is filled with data, it is send to the IoT server via a TCP socket. The flowchart of the embedded signal processing algorithm is presented in Fig. 2. At each data transmission iteration, we proposed an algorithm to select eventual passages within the acoustic channel: the RMS value is computed for each recorded signal and the 50% and 5% quantiles are computed. The data is transmitted to the IoT server only if their ratio is inferior to 3dB. Hence, the algorithm decides if the signal is to be sent for further signal processing to the NUC server. The entire acquisition process is real-time thanks to the possibility to program the RP board using POSIX thread and sockets libraries.

B. IoT server

Regarding the processing stage on the NUC server, the received chirps from the RP board are matched filtered (MF). The reference signals used in the matched filtering are computed as in [START_REF] Vasile | Online software nonlinearity correction for wideband active ultrasound monitoring systems[END_REF] over the most recent 200 received signals, making the fish detection application adaptive to the current underwater environment. After MF, the results are expressed in terms of RMS amplitudes. Due to the fish passage between two sensors (Tx/Rx), the RMS amplitude drops in value because of propagation channel attenuation. This way, it allows fish counting continuously. Finally, Fig. 3 illustrates the number of files (which corresponds to each detection) stored by the IoT server during a 4-days measurement campaign.

III. RESULTS AND DISCUSSION

This paper investigates the first results obtained using the proposed real-time wideband acoustic fish monitoring system experimented in February 2018 at the EDF R&D headquarters in Chatou, France. Several types of targets that were inspected (wood branches, rocks, live bait fish and salmon) are illustrated in Fig. 4. One of the main objectives of this campaign was the study of trout acoustic signatures, also. In order to investigate weather the spectral diversity can be able to distinguish different types of targets passing through an acoustic barrier, Fig. 5 shows the spectrum of each of the four target signatures. One can notice a clear decrease in energy around 450 kHz for the fish passages, only. This indicates the capabilities to distinguish different types of targets using wideband signal processing. Regarding the trout acoustic signatures, statistical analysis of spectra reveals the decrease in energy, which is linked to each fish passage, is consistent over several representative independent measurements. Fig. 6 illustrates the percentile analysis (median -50%, maximum -95 % and minimum -5 %). In case of a clear acoustic channel, each spectrum superposes almost perfectly, indicating thus a performant calibration algorithm was applied.

In our experimental setup, the validation of the fish acoustic monitoring measurements has been carried out using an underwater video camera which operates independently and continuously throughout the experiment. Fig. 7-(a) shows the expert visual validation procedure: the RMS signal is analysed in conjunction with the recorded video stream.

In order to be able to perform a more detailed analysis of the recorded acoustic signatures, sub-band decomposition has been applied.

The received signal is band pass filtered such as: For each band, the RMS value has been computed for each received ultrasound burst. Such RMS curves computed by subbands of 50 kHz, in the nominal bandwidth of the system, are illustrated in the Fig. 7-(b).

Several parameters are stored and analysed for recognition purposes (chirp rate and nonlinearity coefficients). The monitoring of the channel (2 × 1.5 × 10 m), with several trouts within, shows a clear pattern of the movement of fish during the 12 hours of monitoring (at night): detections in the first hours of deployment and less activity toward morning because the trouts accommodated to their new environment. These observations are consistent throughout all sub-bands, hence the reliability of the RMS based detection scheme has been demonstrated.

IV. CONCLUSION

In conclusion, this paper proposed a novel IoT active acoustic system for fish monitoring with direct application in biodiversity preservation near hydropower plants and dams.

The system uses a cost effective single board computer for data acquisition and it is processing acoustic signals at two levels: the first at the board level (embedded), to decide upon transmission of data to a server, and the server level to show the results of detection and the number of fish passages (acoustic signature analysis).

One future direction for improving the detection results is to employ a spatially adaptive clustering [START_REF] Vasile | Intensitydriven-adaptive-neighborhood technique for POLSAR parameters estimation[END_REF], [START_REF] Vasile | Imagerie Radar à Synthèse d'Ouverture interférométrique et polarimétrique[END_REF] or likelihood ratio testing [START_REF] Besic | Stochastic approach in wet snow detection using multitemporal SARdata[END_REF] of the acoustic signatures based on the total energy.

We also conclude that this solution is suited to large scale deployment given its relatively low cost and energy effectiveness. This is achieved by taking advantage of the last IoT technological capabilities both in terms of hardware and software.

ACKNOWLEDGMENT

This work has been developed in the framework of the ACOUEAU project, supported by the European Union Regional Development Fund, the SCANUS project, supported by the SATT LINKSIUM Grenoble Alpes, and the Scientific Equipment 2017 AO project, supported by the LabEx OSUG@2020. The authors would like to thank to Gelu Ionescu, Steeve Zozor, Amaury Nègre and Teodor Petrut, form the GIPSA-lab, and to Géry Hachet and Alexandre Girard, from the EDF R&D, for their precious support in hardware and software development.

Fig. 1 .

 1 Fig. 1. Setup for fish counting experiments.

Fig. 2 .

 2 Fig. 2. Flowchart of the embedded signal processing algorithm.

Fig. 3 .

 3 Fig. 3. Number of files stored on the IoT server during the Chatou measurement campaign.

Fig. 4 .

 4 Fig. 4. Different targets selected for the measurement campaign: (a) Branche -small wood branch, (b) Cailloux -small rock, (c) Gardon -live fishing bait and (d) Saumon -salmon.

Fig. 5 .

 5 Fig. 5. Spectrum of the acoustic signatures for each target.

Fig. 6 .

 6 Fig. 6. Statistical analysis of the acoustic signatures (spectrum): (left) for the clear acoustic channel and (right) for the trout target.

Fig. 7 .

 7 Fig. 7. RMS curves showing passages of trouts: (a) as validated by the underwater video camera and (b) one night results and the corresponding sub-band decomposition.


Nicolas Le Moine

A spectral approach to the Shallow Water Equations, and its implication for At-Many-stations Hydraulic Geometry

Keywords: Shallow Water Equations, Spectral modeling, River morphology, In-channel landform sequence, Backwater curve, At-Many-stations Hydraulic Geometry

A spectral approach to the Shallow Water Equations, and its implication for At-Many-stations Hydraulic Geometry

Introduction

Solving the Shallow Water Equations (SWE) in 1-D or 2-D geometry requires inputs such as along-stream channel and floodplain bathymetry, friction parameters, and reach boundary conditions (typically, upstream inflow hydrograph and downstream rating curve). Bathymetry and friction data are not always easily available at a sufficient spatial resolution, e.g. for discharge-stage modeling at the scale of a full hydrographic network.

On the other side, from a geomorphological point of view it is unsatisfactory to have to fall back on crude asumptions such as uniform geometry without a longitudinal variability, regardless of the particular shape chosen (be it rectangular or more complex). Indeed, natural rivers do exhibit such variability as one moves downstream, an observation coined downstream hydraulic geometry (DHG ; see e.g. [START_REF] Knighton | River channel adjustment -the downstream dimension[END_REF]). Most noticeably, channel slope decreases and channel width increases for a given, reference discharge quantile along-stream. However, a much shorter-scale (higher spatial frequency) variability is superimposed on this large-scale trend, with alternating facies or morphological units: shallow, fast-flowing stretches (riffles, or crossings between two meanders) alternate with deeper, low-velocity areas (pools, e.g. in the outsides of the bends of meanders). Such sequences typically occur at intervals of 5 to 10 times the bankfull width. For a given discharge, considered a constant at the reach scale, the variability in river cross-sections yields a variability in hydraulic variables (width, depth, hydraulic radius, velocities, etc.) but also in the energy slope. Consequently, any geometric variability "smoothed out" by a uniform approach will come back in the form of a bias in effective friction parameters, thus reducing the robustness of the model.

The parameterization of this variability of rating curves at the scale of a few morphological sequences, coined atmany-stations hydraulic geometry (AMHG), is a new research topic (see e.g. [START_REF] Gleason | Theoretical basis for at-many-stations hydraulic geometry[END_REF]). However, very few approaches based on the simple, 1-D steady-state flow equation (backwater curve) have tried to even qualitatively explain the emergence of AMHG and the longitudinal variability of rating curves at the scale of a few times the bankfull width. It is nonetheless likely that AMHG is the result of a coevolution of geometrical bed parameters [START_REF] Shen | Geomorphological significance of at-many-stations hydraulic geometry[END_REF], and of an emerging covariance structure between these parameters at the reach scale. We can take advantage of this covariance structure to reduce the dimensionality of the bathymetric model, without smoothing out the high frequency variability as in the uniform approach. At the extreme end, it can be represented by strictly periodic functions, allowing at the same time a quasi-analytical resolution of the backwater curve equation (Fourier transform, resolution in the frequency domain, and inverse transform), sending upstream and downstream boundary conditions to ±∞.

Outline of the study

We test this spectral treatment of AMHG on a dataset of 17 reaches of small French rivers extensively studied in [START_REF] Navratil | Débit de pleins bords et géométrie hydraulique : une description synthétique de la morphologie des cours d'eau pour relier le bassin versant et les habitats aquatiques[END_REF]. In our approach, a reach is modeled as a succession of trapezoidal cross-sections, with periodic series for talweg depth anomaly, channel bank and floodplain lateral slope, banfkfull depth and bankfull width. The mean spatial period of the sequences is obtained using wavelet analysis. The approach allows us to reproduce a pattern of variability in rating curves at successive cross-sections, with only a handful of parameters describing the along-stream covariance structure of the main cross-sectional parameters. This parcimonious model is well-suited for inverse problems in which the inflow discharge has to be estimated using the sole measurement of water level, without a knowledge of low-flow bathymetry and friction parameters (e.g. in remotely-sensed discharge estimation).

We propose perspectives to extend the results obtained in 1-D in the case of a 2-D steady-state flow, using a conformal mapping of a periodic half-floodplain into a rectangular, computational domain where a spectral solution is to be calculated.


M Evelt

L Soumah

A B Rinkevich

O Demokritov

A Anane

V Cros

Jamal Ben Youssef

G De Loubens

O Klein

P Bortolotti

S O Demokritov

V E Demidov
email: demidov@uni-muenster.de

Emission of coherent

Recent advances in the studies of the spin-orbit torque [START_REF] Hoffmann | Spin Hall Effects in Metals[END_REF][START_REF] Sinova | Spin Hall effects[END_REF][START_REF] Soumyanarayanan | Emergent phenomena induced by spin-orbit coupling at surfaces and interfaces[END_REF][START_REF] Hellman | Interface-induced phenomena in magnetism[END_REF] opened novel perspectives for nano-spintronics by enabling utilization in spintronic devices of insulating magnetic materials such as Yttrium Iron Garnet (YIG), which is famous for its unprecedented small magnetic damping [START_REF] Cherepanov | The saga of YIG: Spectra, thermodynamics, interaction and relaxation of magnons in a complex magnet[END_REF]. The small natural damping is crucial for the reduction of power consumption of spintronic devices based on the spin-transfer torque phenomena [START_REF] Slonczewski | Current-driven excitation of magnetic multilayers[END_REF][START_REF] Berger | Emission of spin waves by a magnetic multilayer traversed by a current[END_REF]. Additionally, the small damping is highly beneficial to applications utilizing propagating magnons as carriers of spin information, since it results in large magnon propagation length [START_REF] Kajiwara | Transmission of electrical signals by spin-wave interconversion in a magnetic insulator[END_REF][START_REF] Cornelissen | Longdistance transport of magnon spin information in a magnetic insulator at room temperature[END_REF][START_REF] Yu | Approaching soft X-ray wavelengths in nanomagnetbased microwave technology[END_REF][START_REF] Liu | Long-distance propagation of short-wavelength spin waves[END_REF].

It is now well established that spin-orbit torque (SOT) in YIG/Pt systems can be used to enhance the propagation of coherent magnons excited by microwave-frequency field [START_REF] Wang | Control of Spin Waves in a Thin Film Ferromagnetic Insulator through Interfacial Spin Scattering[END_REF][START_REF] Padron-Hernandez | Amplification of spin waves in yttrium iron garnet films through the spin Hall effect[END_REF][START_REF] Evelt | High-efficiency control of spin-wave propagation in ultra-thin yttrium iron garnet by the spin-orbit torque[END_REF] and to achieve the complete compensation of the natural damping resulting in the spontaneous excitation of coherent magnetization auto-oscillations [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF][START_REF] Demidov | Magnetization oscillations and waves driven by pure spin currents[END_REF][START_REF] Safranski | Spin caloritronic nano-oscillator[END_REF]. It is also experimentally and theoretically ascertained, that, in the regime of small driving currents below the damping compensation point, SOT causes excitation of incoherent magnons with a broad spectral distribution [START_REF] Demidov | Control of Magnetic Fluctuations by Spin Current[END_REF][START_REF] Demidov | Chemical potential of quasi-equilibrium magnon gas driven by pure spin current[END_REF]. These magnons are able to propagate in YIG films to significant distances [START_REF] Cornelissen | Longdistance transport of magnon spin information in a magnetic insulator at room temperature[END_REF][20][START_REF] Cornelissen | Magnon spin transport driven by the magnon chemical potential in a magnetic insulator[END_REF][START_REF] Wesenberg | Long-distance spin transport in a disordered magnetic insulator[END_REF][START_REF] Thiery | Nonlinear spin conductance of yttrium iron garnet thin films driven by large spin-orbit torque[END_REF] and can be used for controllable transmission of spin information [START_REF] Cornelissen | Spin-Current-Controlled Modulation of the Magnon Spin Conductance in a Three-Terminal Magnon Transistor[END_REF]. However, the lack of coherence does not allow one to utilize these magnons for information processing based on coherentwave phenomena, where the phase of the waves plays a decisive role [START_REF] Chumak | Magnon spintronics[END_REF].

Transition from the incoherent to the coherent regime is expected to occur when SOT completely compensates the natural damping in the magnetic film resulting in the excitation of coherent auto-oscillations [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF][START_REF] Demidov | Magnetization oscillations and waves driven by pure spin currents[END_REF], which can be treated as an accumulation of a large number of the SOT-excited magnons in a single spectral state [START_REF] Demidov | Chemical potential of quasi-equilibrium magnon gas driven by pure spin current[END_REF]. However, this process is systematically observed in experiments to be accompanied by a nonlinear shift of the magnon frequency, which causes a spatial self-localization of coherent magnons preventing their emission into the surrounding magnetic film [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF][START_REF] Demidov | Magnetization oscillations and waves driven by pure spin currents[END_REF][START_REF] Thiery | Nonlinear spin conductance of yttrium iron garnet thin films driven by large spin-orbit torque[END_REF][START_REF] Demidov | Direct observation of dynamic modes excited in a magnetic insulator by pure spin current[END_REF].

Here we experimentally demonstrate an approach, which allows one to overcome the self-localization phenomena and achieve an efficient emission of SOT-excited coherent magnons into an extended magnetic-insulator film. We show that this can be accomplished by utilizing nanometer-thick Bi-substituted YIG films exhibiting large perpendicular magnetic anisotropy (PMA) [START_REF] Soumah | Ultra-low damping insulating magnetic thin films with large perpendicular magnetic anisotropy get perpendicular[END_REF]. The effects of PMA result in a compensation of the in-plane shape anisotropy of the film. Due to this compensation, the nonlinear frequency shift vanishes leading to the current-independent frequency of the spontaneously formed coherent magnon state. As the result, the created state can emit propagating magnons with the wavelength of about 300 nm. We demonstrate that the Bi doping necessary to achieve sufficiently strong PMA has no significant adverse effect on the density of the electrical current necessary to excite SOT-induced autooscillations and on the magnon propagation length in YIG. This makes the suggested approach very promising for insulator-based spintronic and magnonic applications.

The schematic of our experiment is shown in Fig. 1. The studied devices are based on an extended 20-nm thick film of Bi-doped YIG (BiYIG) (Bi1Y2Fe5O12) with PMA grown by the pulsed laser deposition on substituted Gallium Gadolinium Garnet (sGGG) substrate [START_REF] Soumah | Ultra-low damping insulating magnetic thin films with large perpendicular magnetic anisotropy get perpendicular[END_REF] (see Supplementary information for film characterization). A 6nm thick Pt film is sputter-deposited onto BiYIG and is patterned into a shape of a strip line with the width of 1 μm and the length of 4 μm. Due to the spin-orbit interaction in Pt [START_REF] Hoffmann | Spin Hall Effects in Metals[END_REF][START_REF] Sinova | Spin Hall effects[END_REF] causing the scattering of conduction electrons with opposite magnetic moments toward opposite surfaces of the Pt film, the dc electric current I flowing in the plane of the strip line [28] is converted into the out-of-plane spin current Is (see inset in Fig. 1).

The spin current is injected into the BiYIG film and exerts a torque on its magnetization M. The saturating static magnetic field H0 is applied in the device plane and perpendicular to the direction of the current flow. For the orientation of H0, as shown in Fig. 1, the positive current I results in SOT compensating the natural damping in the BiYIG film [START_REF] Ando | Electric Manipulation of Spin Relaxation Using the Spin Hall Effect[END_REF]. When the damping is completely compensated, large-amplitude coherent magnetization auto-oscillations are excited in the magnetic film above the Pt line. The auto-oscillations can emit coherent magnons into the surrounding BiYIG film, provided that the frequency of the auto-oscillations is not shifted below the spectrum of propagating magnons due to the nonlinear frequency shift.

We detect the SOT-induced magnetization oscillations and the emitted magnons by the space-and time-resolved micro-focus Brillouin light scattering (BLS) spectroscopy [START_REF] Demidov | Magnonic Waveguides Studied by Microfocus Brillouin Light Scattering[END_REF]. We focus the probing laser light with the wavelength of 532 nm and the power of 0.1 mW into a diffraction-limited spot on the surface of the BiYIG film through the sample substrate (Fig. 1) and analyze the spectrum of inelastically scattered light. The resulting BLS signal is proportional to the density of magnons or, alternatively, to the intensity of magnetization oscillations at the position of the probing spot. By rastering the probing-light spot over the surface of the sample, the spatially resolved maps of the magnon density can be obtained.

In Figure 2(a), we present the color-coded BLS intensity in the current-frequency coordinates measured by placing the probing spot in the middle of the Pt line. We find that, at I≥1.8 mA, the BLS spectra exhibit a narrow intense peak (see also inset in Fig. 2(a)) corresponding to the SOT-induced auto-oscillations in the BiYIG film. The onset current IC=1.8 mA corresponds to the current density in the Pt layer of 3×10 11 A/m 2 , which is close to the value obtained for devices based on undoped in-plane magnetized YIG films [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF]. In contrast to the previously demonstrated devices, where the frequency of auto-oscillations strongly decreased with the increase of I due to the nonlinear frequency shift [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF][START_REF] Demidov | Direct observation of dynamic modes excited in a magnetic insulator by pure spin current[END_REF], in BiYIG-based devices, the frequency of the oscillations stays nearly constant over the entire range of I, while the intensity of auto-oscillations monotonously increases with the increase of the driving current (Fig. 2(b)). The frequency of the auto-oscillations fA=5.56 GHz is close to the independently determined frequency of the ferromagnetic resonance (FMR) in the BiYIG/Pt bilayer f0=5.63 GHz (see Supplementary information). Therefore, one can characterize the nonlinear frequency shift by utilizing the Kittel formula for the FMR frequency in an in-plane saturated magnetic film with PMA [START_REF] Kalinikos | The dipoleexchange spin wave spectrum for anisotropic ferromagnetic films with mixed exchange boundary conditions[END_REF]

0 0 0 ( (4 )
)

a e f H H N M γ π = + - . ( 1
)
Here γ is the gyromagnetic ratio, Since the nonlinear frequency shift is known to be responsible for the localization of coherent auto-oscillations, one expects that, in the studied devices, the fulfilment of the magnon propagation conditions will be restored leading to emission of coherent magnons by auto-oscillations excited in the active device area. Indeed, spatially resolved BLS measurements (Fig. 3(a)) confirm the efficient magnon emission. Figure 3(a) shows the color-coded spatial map of the intensity of SOT-induced magnetization oscillations recorded by rastering the probing laser spot over 20 μm by 7 μm area centered at the Pt line (dashed lines in Fig. 3(a)). As seen from these data, the oscillations excited due to SOT are not localized in the active area above the Pt line, but significantly extend into the surrounding BiYIG film suggesting the emission of excited coherent magnons. Figure 3(b) shows the x-dependence of the integrated BLS intensity on the log-linear scale. The dependence exhibits a well-defined exponential decrease enabling determination of the mean propagation length of the emitted magnons ξ=3.7

μm. The spatially resolved measurements also allow us to determine the group velocity of the magnons by comparing the temporal dependences of the BLS intensity at x=0 and 10 μm (Fig. 3(c)). As seen from these data, the pulse of the dynamic magnetization detected at x=10 μm is delayed with respect to that at x=0 by about 58 ns, which corresponds to the group velocity Vg≈0.17 μm/ns.

To get insight into the nature of excited magnons, we calculate the magnon dispersion spectrum in the BiYIG film using the analytical theory developed in Ref.

[31] and independently determined anisotropy field Ha=1.84 kOe (see Supplementary information). Figure 4(a) shows the calculated dispersion curves for magnons propagating parallel (φ=0) and perpendicular (φ=90°) with respect to the direction of the static magnetic field. The dashed horizontal line marks the experimentally determined frequency of the auto-oscillations fA (Fig. 2(a)). As seen from Fig. 4(a), fA is located within the spectrum of propagating magnons allowing their emission from the active device area. We emphasize that magnons with the frequency fA can propagate in all directions with respect to the static field (see the isofrequency curve in the inset in Fig.

4(a)

), which is in agreement with the experimental emission pattern (Fig. 3 Finally, we discuss the dependence of the characteristics of emitted magnons on the magnitude of the static magnetic field H0. Figure 5(a) shows the field dependence of the experimentally determined magnon propagation length ξ and its inverse value " 1/ k ξ = -the imaginary part of the magnon wavevector. As seen from these data, the propagation length significantly increases at small H0, which is associated with the reduction of the magnon relaxation frequency:

0 0 r r H ω ω αγ = +
, where α is the Gilbert damping parameter and 0 r ω is a nonzero offset originating from the influence of the inhomogeneity of the magnetic film. In agreement with the theory [START_REF] Gurevich | Magnetization Oscillations and Waves[END_REF],

0 0 g g " r H k V V ω γ α = +
varies linearly with H0 (line in Fig. 5(a)). Taking into account the experimentally determined group velocity Vg=0.17 μm/ns, which was found to be nearly constant within the range H0=1-2 kOe, we calculate from the slope of the experimental dependence k''(H0) the damping parameter for the BiYIG film α≈8.2×10 -4 . This value is not much larger than that typical for undoped YIG films of the same thickness [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF][START_REF] Collet | Spin-wave propagation in ultra-thin YIG based waveguides[END_REF]. However, it is significantly smaller in comparison with the values typical for metallic magnetic films, making BiYIG a promising material for insulator-based spintronic applications. Note that the found value characterizes the damping in a free BiYIG film, while, in a BiYIG/Pt bilayer, the damping is expected to be significantly increased due to the influence of the spin pumping effect. To characterize the damping in the active device area, we analyze the field dependence of the auto-oscillation onset current (Fig. 5(b)), which is proportional to the relaxation frequency [START_REF] Slavin | Nonlinear Auto-Oscillator Theory of Microwave Generation by Spin-Polarized Current[END_REF]:

C 0 P t / Y I G 0 ( ) r I H β ω α γ = +
, where β is the proportionality coefficient describing the SOT efficiency. Using the results of the previous analysis, we find the Gilbert damping constant αPt/YIG≈1.6×10 -3 . Since in a bilayer constituted by a nanometer-thick lowdamping insulator and Pt, the damping is mostly determined by the spin pumping, the above value is very close to that found for bilayers based on undoped YIG [START_REF] Collet | Generation of coherent spin-wave modes in yttrium iron garnet microdiscs by spin-orbit torque[END_REF]. This fact elucidates the similarity in the auto-oscillation onset current density between the devices studied in our work and those based on undoped YIG films.

In conclusion, we have experimentally demonstrated devices capable of emission of coherent spin-orbit torque excited magnons into magnetic-insulator films. We show that this can be achieved by using films with perpendicular magnetic anisotropy, where the effects of the anisotropy allow one to overcome the nonlinear localization of magnons. The demonstrated devices are flexible with respect to modifications and can be used as coherent sources for magnon-based nanoscale information transmission and processing. We believe that our results will stimulate further developments in the emerging field of insulator-based spintronics. Solid line -linear fit.

  is the magnetization-independent demagnetization factor associated with the anisotropy, a H is the effective anisotropy field, and e M is the effective static magnetization, which decreases with the increase in the current in the Pt layer resulting in the variation of the frequency f0. For BiYIG films used in our study 4 a N π ≈ . Therefore, the e M -dependent term in Eq. (1) vanishes resulting in the vanishing nonlinear frequency shift.

  (a)). Note that, as shown by Fig.4(b), magnons emitted at different angles φ are expected to possess significantly different wavelengths and group velocities. In particular, magnons emitted in the direction perpendicular to the Pt line (φ=0) exhibit the shortest wavelength of about 340 nm and the largest group velocity of 0.21 μm/ns. The latter value is reasonably close to that obtained in the experiment, which confirms the validity of our calculations.

Figure 1 .

 1 Figure 1. Experimental layout. The device consists of an extended 20-nm thick BiYIG film with PMA and 6-nm thick Pt line with the length of 4 μm and the width of 1 μm. The dc current I in Pt is converted into the pure spin current Is (see the inset) inducing local coherent magnetization auto-oscillations, which emit propagating coherent magnons into the surrounding BiYIG film.

Figure 2 .

 2 Figure 2. (a) Color-coded BLS intensity in the current-frequency coordinates measured by placing the probing spot in the middle of the Pt line. fA marks the frequency of SOTinduced auto-oscillations. Inset shows the representative BLS spectrum for I=2.5 mA. (b) Current dependence of the peak BLS intensity. IC marks the current corresponding to the onset of auto-oscillations. The data were obtained at H0=2.0 kOe.

Figure 3 . 6 Figure 4 .

 364 Figure 3. (a) Color-coded spatial map of the BLS intensity recorded by rastering the probing laser spot over 20 μm by 7 μm area. Dashed lines show the contours of the Pt line. (b) x-coordinate dependence of the integrated BLS intensity on the log-linear scale. Symbols -experimental data, solid curves -exponential fit at |x|>2 μm. Vertical dashed lines show the boundaries of the Pt line. (c) Temporal dependences of the BLS intensity recorded at x=0 and 10 μm. Dashed vertical lines mark the beginning of the rear edge of the pulses used to determine the magnon propagation delay. The data were obtained at H0=2.0 kOe and I=2.5 mA.

Figure 5 .

 5 Figure 5. (a) Field dependence of the experimentally determined magnon propagation length and of the imaginary part of the magnon wavevector k''. Symbols -experimental data. Solid line -linear fit of the data for k''. The data were obtained at I=2.5 mA. (b) Field dependence of the auto-oscillation onset current. Symbols -experimental data.

ACKNOWLEDGMENTS

We thank C. Carrétéro and E. Jacquet for their contribution to films growth. We acknowledge support from Deutsche Forschungsgemeinschaft, FASO of Russia (theme "Spin" No. AAAA-A18-118020290104-2), Russian ministry of Education and Science (project No. 14.Z50.31.0025), and that of the ANR Grant ISOLYIG (ref 15-CE08-0030-01). LS is partially supported by G.I.E III-V Lab. France.