Burak Yüksel
email: burak.yueksel@tuebingen.mpg.de

Saber Mahboubi

Cristian Secchi
email: cristian.secchi@unimore.it

Heinrich H Bülthoff

Antonio Franchi
email: afranchi@laas.fr

Design, Identification and Experimental Testing of a Light-weight Flexible-joint Arm for Aerial Physical Interaction

In this paper we introduce the design of a lightweight novel flexible-joint arm for light-weight unmanned aerial vehicles (UAVs), which can be used both for safe physical interaction with the environment and it represents also a preliminary step in the direction of performing quick motions for tasks such as hammering or throwing. The actuator consists of an active pulley driven by a rotational servo motor, a passive pulley which is attached to a rigid link, and the elastic connections (springs) between these two pulleys. We identify the physical parameters of the system, and use an optimal control strategy to maximize its velocity by taking advantage of elastic components. The prototype can be extended to a light-weight variable stiffness actuator. The flexible-joint arm is applied on a quadrotor, to be used in aerial physical interaction tasks, which implies that the elastic components can also be used for stable interaction absorbing the interactive disturbances which might damage the flying system and its hardware. The design is validated through several experiments, and future developments are discussed in the paper.

I. INTRODUCTION

The physical interaction of flying robots, a.k.a. aerial physical interaction, has been a growing interest of the scientists and engineers in the last decade. The control of flying robots during physical interaction has become a challenge in terms of designing new tools or actuators, as well as developing powerful algorithms to allow the exertion of forces and torques on the environment while stabilizing the overall system and protect the expensive hardware.

There are various algorithms presented for aerial physical interaction. In [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF], [START_REF]A nonlinear force observer for quadrotors and application to physical interactive tasks[END_REF] an interconnection and damping assignment passivity-based controller is presented, where the algorithm assigns desired physical properties to the flying robot (quadrotor) depending on the force and torques coming from the environment. In [START_REF] Lee | Mechanics and control of quadrotors for tool operation[END_REF], [START_REF] Nguyen | Hybrid force/motion control and internal dynamics of quadrotors for tool operation[END_REF] passive decomposition methods have been designed for controlling the trajectory of a rigid tool onboard of a quadrotor. Cartesian impedance control for a manipulator placed on an UAV is shown in [START_REF] Lippiello | Exploiting redundancy in cartesian impedance control of UAVs equipped with a robotic arm[END_REF]. An impedance controller for aerial manipulation performed by a ducted-fan vehicle is presented in [START_REF] Forte | Impedance control of an aerial manipulator[END_REF].

In parallel, other studies presented useful tools for improving the performance of physical interaction tasks. In [START_REF] Kondak | Closed-loop behavior of an autonomous helicopter equipped with a robotic arm for aerial manipulation tasks[END_REF], a light-weight industrial arm is attached to a small-size helicopter, under the consortium of [START_REF] Arcas | EU Collab[END_REF]. Smaller scale designs are also presented in the literature. In [START_REF] Kim | Aerial manipulation using a quadrotor with a two dof robotic arm[END_REF] a 2D rigid arm for aerial manipulation is introduced. A tool for surface inspection using a flying robot is developed and presented in [START_REF] Fumagalli | Modeling and control of a flying robot for contact inspection[END_REF]. Besides many different joints and actuated tools, passive ones are also used for physical interaction, as in [START_REF] Gioioso | Turning a near-hovering controlled quadrotor into a 3D force effector[END_REF].

In aerial physical interaction, it is necessary to always guarantee a safe (i.e., non-destabilizing) behavior of the system during any contact with the environment, that can be either desired or unforeseen (e.g., the robot moving in an unknown and/or hostile area). Furthermore, it should be possible to implement explosive motions for useful aerial interaction like aerial repairing or fixing. The use of flexible joints has proven to be successful for the implementation of interactive tasks when using grounded robotic arms or humanoid robots. The elasticity of the joints can in fact be exploited for achieving an intrinsic safe behavior of the system and for amplifying the mechanical performance of a rigid arm by exploiting the energy stored into the elastic element [START_REF] Braun | Optimal variable stiffness control: formulation and application to explosive movement tasks[END_REF]. To the best of the authors' knowledge, the possibility of exploiting the benefits of elastic joints in aerial interaction has not been explored yet. The goal of this paper is to conduct a preliminary work to start filling this gap. We present a novel design of a light weight flexible-joint arm that can be mounted on a small size aerial vehicle for achieving an intrinsically safe aerial physical interaction and that allows to exploit the joint elasticity for implementing aerial explosive tasks (e.g., aerial hammering and aerial throwing). The proposed arm is then mounted on a quadrotor and its benefits during aerial physical interaction are experimentally validated.

The paper is organized as follows: in Section II we present the mechanical design and the mathematical modeling of the flexible-joint arm. The parameters of the linear models for both servo motor and flexible-joint arm dynamics are found by system identification methods and validated using experimental data. Section III is dedicated to the control framework of both flexible-joint arm and quadrotor, separately. The flexible-joint arm is controlled via an unconstrained optimal control strategy to maximize the rigid link velocity by controlling the servo motor velocity and using the potential energy stored on the elastic components. The quadrotor is controlled using near hovering algorithm, which is an effective way of controlling quadrotor state, while it can also be used easily by a human operator for teleoperation tasks. Section IV gives detailed information about the hardware used in flexible-joint arm design, and quadrotor itself. In the same section we explain the experimental setup, and present the results of this study. Finally in Section V we discuss the results of this paper, give useful remarks, and propose future developments on the current design.

II. DESIGN, MODELING, AND IDENTIFICATION

In this section we present the design of the flexible-joint arm in detail, from the CAD model to the actual product. Afterwards we derive the mathematical model of the system, and find the physical parameters in this model by using experimental data with system identification methods. The identified model, together with real measurements, is then exploited for evaluating the evolution of the system.

A. Design of the Flexible-joint Arm

The flexible-joint arm consists of several parts. The rigid parts, with the exception of the actuators, are CAD modeled and 3D printed. As can be seen in Figs. 2 and3, a rigid pulley is attached to the shaft of the servo motor, which is connected to a second rigid pulley via two elastic elements (see also Fig. 1). The first pulley is referred to as active pulley and the second one as passive pulley. The passive pulley is attached to a rigid link, whose objective is to interact with the q πθ environment. The position of the active pulley is measured by the encoder of the servo motor. The measurements for the arm motion are collected with a magnetic encoder attached to the passive pulley. The magnet is placed on a cylindric part attached to the center of the passive pulley, and the encoder is placed on a fixed surface.

A second linear servo motor can let the active pulley slide along the bottom surface in order to regulate the distance between the two pulleys. At present however this feature is not used and the distance has been fixed with a rigid connection. The reason for introducing a second linear motor in the design is to have room for future improvements, such as changing the design from an elastic actuator to a variable stiffness actuator by using, e.g., nonlinear elastic components [START_REF] Ham | Compliant actuator designs[END_REF]. Finally, a rigid box covering the mechanism except the rigid arm is adopted for protection of the hardware components.

The elastic components between two rigid pulleys are chosen as linear springs (see in Figs. 1 and4). The two springs are used as antagonistic pairs, i.e., in case one of the springs contracts, the other one relaxes, and viceversa. This is natural way of designing a flexible arm, similar to biceps and triceps in the human arm; where the muscels cannot push, but only pull [START_REF] Ham | Compliant actuator designs[END_REF]. The advantage of using antagonistic pairs is that, they allow the full elastic behavior in both directions of the rotation.

Besides the mechanics of the system, the electronics and communication with the arm is another part of the design. This is explained in SectionIV-A in detail.

B. Dynamic Model and Identification

The physical parameters of the flexible arm are mostly unknown. In this part of the paper we derive these parameters that will be later used for controlling the flexible-joint arm motion. Additionally, we derive the models of the link and of the motor that will be exploited for evaluating the dynamic behavior of the system during fast motions.

The flexible-joint arm consists of motor-side and link

Preprint version, final version at http://ieeexplore.ieee.org/ 2 2015 IEEE ICRA side dynamics. Both are of second order as shown in [START_REF] Ozparpucu | Optimal control for maximizing link velocity of visco-elastic joints[END_REF]. However, if the motor dynamics is considerably faster than the link dynamics, it is possible to assume that the motor velocity can be controlled instantaneously [START_REF] Ozparpucu | Optimal control for maximizing link velocity of visco-elastic joints[END_REF], [START_REF] Haddadin | On impact decoupling properties of elastic robots and time optimal velocity maximization on joint level[END_REF]. This means that the servo motor behaves as a perfect velocity source. Under this assumption, the linear dynamics of the flexible-joint arm can be described as

M q + D q q + K J q = D θ θ + K J θ (1)
θ = t 0 θd dt + θ 0 (2)
where θ ∈ R and q ∈ R are motor side and link side positions, respectively, M > 0 is the inertia of the rigid link, K J > 0 is the linear spring constant, and D θ > 0 and D q > 0 represent the dissipations (viscous elements) of the motor side and link side, respectively. The desired motor velocity is shown with θd , and θ 0 is the initial motor position. The validity of the assumption of considering motor as a velocity source depends on the load attached to the motor. Since we are introducing a novel design, it is in our interest to see whether the servo motor can be used as velocity source when the elastic and rigid parts are also attached.

Considering as input the desired motor velocity θd , and as output the measured motor velocity θ, the transfer function, in the Laplace domain, is expected to be of the first order, namely

Θ(s) Θd (s) = e -sτ a s + b . ( 3
)
where τ represents the system delay.

In order to validate the model (3) and to estimate the parameters a, b, and τ we performed some experiments and used the non-linear least squares method [START_REF] Ljung | System Identification: Theory for the User[END_REF]. In Fig. 5 the frequency response of the motor is shown for different step inputs as desired velocity profile and with different conditions: loaded and unloaded case. The output is the measured motor velocity that is retrieved using the motor encoder. The different plot colors correspond to pulse trains of different frequencies. Solid lines are corresponding to conditions where no load is attached to the motor. Dashed lines represent the cases when the motor is connected to the arm though the elastic interconnection. As it can be seen, the frequency response of the system stays almost the same in the frequency range of interest, which is around the natural frequency of the total system (this will be identified in a next step). Hence we can consider the motor as velocity source even when it is connected to its load. A good fit for (3) considering different conditions is found for the parameters a = 13.79 and b = 7.175, and τ = 0.2s.

The parameters of the flexible-joint arm given in (1) have to be estimated as well. Let us first write the transfer function of the system dynamics for motor position θ as input and link position q as output. Denoting with ω n and G the natural frequency and the low-frequency gain of the system, respectively, we obtain:

Q(s) Θ(s) = e -sτ K s + µ s 2 + 2ξω n + ω 2 n (4)
where K = GD θ /M , µ = K J /D θ , ω n = K J /M , and ξ = D q /2 √ K J M . The moment of inertia for the rigid link is computed as M = 0.0019 kg•m 2 from the CAD model shown in Fig. 2. Despite the inertia is easy to compute from the geometry of the system, the dissipative parameters such as damping and frictions are hard to retrieve from a simulation. For this reason, we used non-linear least squares [START_REF] Ljung | System Identification: Theory for the User[END_REF] to identify the system parameters. By choosing motor position θ as input and link position q as output we computed the frequency response of the system for different step input profiles. The results are presented in Fig. 6. The best fit for (4) is

Q(s) Θ(s) = e -0.
The physical parameters corresponding to this transfer function are shown in Table I. In Figs. 5 and6 it is clear that the both motor and link can be operated around the natural frequency of the flexible- Motor Velocity Fig. 7: Motor position and velocities. Blue dashed line represents the optimal control input computed using (9) for t f = 3s and ωn = 13.3rad/s. Black dots are the measured motor position and velocities using embedded servo motor encoder in 10 Hz. Red plot are for the corresponding evolutions of the identified transfer functions.

joint arm, ω n , when the servo motor is considered as a perfect velocity source. The natural frequency of the system is especially interesting for us, since we would like to test the new design during fast movement tasks, as will be explained in the next section in detail. Now, by considering the flexible-joint arm as a secondorder LTI system, where desired velocity is input, and motor position and link position are outputs, we have

Θ(s) Θd (s) = e -sτ a s(s + b) . (6)
and

Q(s) Θd (s) = e -sτ (Kas + Kµa) s 4 + (2ξω n + b)s 3 + (ω 2 n + 2ξω n b)s 2 + ω 2 n bs (7)
where both ( 6) and ( 7) can be used to evaluate motor and link motions during explosive movements of the system.

III. CONTROL

A. Control of the Flexible-joint Arm

Due to its compliance, the flexible-joint arm is intrinsically safe and an unforeseen interaction with the environment is absorbed by the arm and it perturbs only slightly the motion of the aerial vehicle it is mounted on. Another interesting interactive feature of flexible-joint arms is the possibility of achieving fast, or even explosive motions. Such movements

Link Velocity

Fig. 8: Link position and velocities. Blue dashed line represents the optimal control input computed using (9) for t f = 3s and ωn = 13.3rad/s. Black dots are the measured link position and velocities using the magnetic encoder in 10 Hz. Red plots are for the corresponding evolutions of the identified transfer functions. can be transformed to hammering or throwing tasks using flying robot, which can be useful especially for construction works in relatively high buildings where a ground robot cannot reach.

By explosive movement, we mean the amplification of link velocity using elastic components of the flexible-joint arm. A rigid arm, directly connected to the actuator without any elastic element, will always have the velocity of the motor. However, a flexible-joint arm can reach higher velocities than actuator can provide, thanks to its elastic components. In the following, we first use an unconstrained optimal control strategy for maximizing the link velocity [START_REF] Haddadin | On impact decoupling properties of elastic robots and time optimal velocity maximization on joint level[END_REF]. Afterwards, the trajectory controller of the flying platform is introduced, which is presented before in [START_REF] Lee | Semi-autonomous haptic teleoperation control architecture of multiple unmanned aerial vehicles[END_REF].

1) Maximizing Velocity of the Flexible-joint Arm: It is shown in the literature [START_REF] Haddadin | On impact decoupling properties of elastic robots and time optimal velocity maximization on joint level[END_REF] that the elastic components, e.g., springs, can be used to amplify the velocity of the actuation source. In this paper we will use unconstrained optimal control strategy to maximize the flexible-joint arm velocity in a specified final time t f . Hence the cost function to be maximized is

J = q(t f ). (8)
Since we found in the previous section the damping factors relatively low, we will neglect them in this paper. Hence, the optimal control policy presented in [START_REF] Ozparpucu | Optimal control for maximizing link velocity of visco-elastic joints[END_REF], [START_REF] Haddadin | On impact decoupling properties of elastic robots and time optimal velocity maximization on joint level[END_REF] for unconstrained undamped mass-spring system is suitable for our purpose. The optimal controller is then

u * = θmax , sin(ω n (t f -t)) > 0 θmin , sin(ω n (t f -t)) < 0 (9)
where the optimal control input is u * = θ * d . Note that we consider the servo motor as perfect velocity source (see Section II) which is however constrained, i.e., u * ∈ [ θmin , θmax ]. This type of controller is called bang-bang controller. In our case no state constraints are considered.

We have applied this controller to the system, on board of a quadrotor, by choosing t f = 3s. The results are shown in Figs. 7 and8. The control input as desired motor velocity is shown in both figures with blue dashed line. The measurements are depicted with black dots, in 10 Hz from both motor and link encoders. Here we also evaluated the evolution of the identified transfer functions in 1 KHz, which is presented with solid red plots. As seen in the lower plot of Fig. 8 the maximum link velocity is reached in around t f = 3s, which is almost four times more than the servo motor velocity.

B. Control of the Aerial Platform

In this part we give an overview of the controller used for quadrotor, to achieve a stable flight while the flexible-joint arm is operating. In addition, we expect from the controller to give a good tracking results when a human operator is controlling its trajectory. In real-life cases, where flying robots are required to be interacting with their environment (such as in out-door scenarios) we believe that controlling the flying system with human input is as important as fully autonomous trajectory control. For this reason, we choose to control the quadrotor using so called near hovering control, as deeply explained in [START_REF] Lee | Semi-autonomous haptic teleoperation control architecture of multiple unmanned aerial vehicles[END_REF].

Briefly, the controller is developed for teleoperation tasks using quadrotors, where it allows operator more focusing on high-level tasks, while low-level ones are hidden from the operator such as underactuation of the translational dynamics. The goal of the controller is to make quadrotor follow a smooth trajectory in translational motion, using the inputs given by the human operator, while trying to keep hovering configuration, i.e., small pitch and roll angles, as much as possible. In this paper we used this controller for controlling the translational trajectory of the quadrotor, where the novel flexible-joint arm is attached to its center of gravity as shown in Fig. 9.

IV. HARDWARE AND EXPERIMENTS

The flexible-joint arm is a combination of 3D printed rigid parts, elastic components, actuator, and measurement units with their electronics. The 3D printed parts and elastic elements are introduced previously in detail. In this section we will present the actuator, i.e., the servo motor, the measurement units, and the electronics. Moreover, we will briefly give information about the flying platform and the setup for the experiments. Finally we present the experimental results.

A. Hardware of the Flexible-joint Arm

The main limitation of the design is the weight and size of the flexible-joint arm, since it is planned to be mounted on a small-scale quadrotor. For this reason we chose Dynamixel AX-12A servo motor, which is both velocity and position controllable. The motor provides velocity, position, load, and temperature measurements. The communication with servo motor is done using serial communication block of Simulink-Matlab, which is running in 10 Hz [START_REF] Mensink | Characterization and modeling of a dynamixel servo[END_REF]. The rest of the simulink file is running in 1 KHz, which includes the controller of the flexible-joint arm. The position and velocity measurements of the passive pulley (so the rigid link) are acquired using 10-bit AEAT-6010 magnetic encoder 1 . The encoder is fixed to the body of the actuator, whereas the rotating part (magnet) is directly connected to the passive pulley. The encoder readings are acquired using serial channel of Arduino-Uno2 , and transferred to another serial communication block in the same simulink file where the servo motor is controlled. The encoder measurements are also done in 10 Hz. Simulink file was running under Windows-8 machine using Matlab-2013b.

B. Experimental Setup

The experiments are conducted on the flying platform, which is a quadrotor UAV 3 . The maximum payload of the quadrotor is about 2 kg, hence we designed the flexiblejoint arm as light as possible (total mass of the flexiblejoint arm is 360 gr). The maximum reachable rotation angle of the flexible-joint arm is limited to q ∈ [-1, 0.55] rad. The experiments are done in-door, and we use a motion capture system (MoCap to have good estimation of the quadrotor position. On board of the quadrotor there are Inertial Measurement Unit (IMU), flight controller, brushless motor controller, battery and markers for the MoCap. The total system weight is 1360 gr. The communication with onboard electronics of the quadrotor and the MoCap is done on an Ubuntu 12.04 machine, using ROS-fuerte [START_REF] Quigley | ROS: an open-source Robot Operating System[END_REF] and Telekyb software, which is developed for controlling the aerial robots [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor control[END_REF].

The quadrotor is controlled by a human operator, using a standard joypad connected to the Ubuntu machine. Besides the operator, another person was holding the security stick (see the videos attached), which is connected to the quadrotor from top. The connection is done with ropes for security reasons, which have no tension the during flight.

C. Experiments and Results

In this paper we designed two experiments; fast movement of the arm, and stable aerial physical interaction task (see Fig. 9).

For the first experiment, the quadrotor is actively keeping a hovering condition, and the optimal controller presented in ( 9) is applied to the flexible arm. The motion of the arm is represented in Fig. 10, which is a snapshot of the real experiment. The results are given in Fig. 8, where the link velocity is clearly amplified compared to desired and measured servo motor velocities given in Fig. 7. Thanks to the storage of elastic energy due to the flexible joint, the motion of the link is amplified and this is the biggest step to take towards the implementation of explosive motions. Once the motion amplification is achieved, it is only necessary to endow the arm with a proper tool for achieving the desired explosive task (e.g., aerial hammering or throwing).

The second experiment is for the aerial physical interaction task using flexible-joint arm on the quadrotor. A platform with a slopped surface (blue colored surface in Fig. 12) is designed for this experiment. The goal is to have a stable flight, while quadrotor equipped with the arm is sliding on the blue surface in safe. For this experiment, we have set the proportional gain of the position controller to zero for near-hovering controller along the y-axis [START_REF] Lee | Semi-autonomous haptic teleoperation control architecture of multiple unmanned aerial vehicles[END_REF]. The reason is to let the quadrotor drift along the direction of motion, while the flexible-joint arm is interacting with the environment (see Fig. 12). The flexible-joint arm is position-controlled in this case, meaning that we send the desired motor position and velocity values to a set value at will. The results are given in Fig. 11. Left side of the figure shows the positions of the quadrotor in three Cartesian axes, with blue solid lines. The gray areas are representing the hovering condition of the quadrotor, meaning that its position is fixed. After the first grayed area, quadrotor starts to descend and approaches to the blue platform so that the flexible arm touches the surface. Later the quadrotor flies along they-axis while trying to keep contact with the sloped surface (see Fig. 12). After holding in hovering position in the second grayed area, it flies back along the + y direction by keeping the contact. Finally it arrives to initial position, ascends, and stays in hovering condition as in the beginning. The pitch angles (rotation along y axis) are given as black solid lines in Fig. 11. The near hovering controller is keeping the system and the interaction stable during the whole flight. The red plots in the same figure are showing the motor positions (θ), and the link positions (q) of the joint arm. Depending on the direction of the flight, different motor positions are set to the servo motor. The changes in the link position is clearly seen in the figure, which is following the profile of the interaction surface. The motor positions, on the other hand, are fixed unless it is not commanded. Thanks to the flexibility of the joint the arm can safely interact with an unknown environment. Furthermore, since the elasticity of the joint absorbs the impact with the environment, the motion of the UAV is left almost unperturbed by the unforeseen impact with the environment.

We encourage the reader to see the video attached to this paper, where we present these two experiments in detail. In addition, two preliminary hammering tests have been shown, using flexible-joint arm attached on a flying quadrotor. In the first hammering experiment, flexible-joint arm is controlled by a standart joystick, while in the second one we use leapmotion sensor 4 to operate the arm based on human hand motion, using ROS-groovy interface.

V. CONCLUSIONS

In this paper we presented a novel light weight flexiblejoint arm for light-weight flying robots, to be used in aerial physical interaction tasks. By using system identification methods, the physical parameters of the new design are estimated and presented. An optimal controller to maximize the elastic link velocity is used to perform fast motion tasks. Such design is not only limited to explosive interaction tasks, but can also be used to achieve stable contact with the environment thanks to its elastic design. The design and the controller are tested on a flying quadrotor and validated with experimental results.

A next step is to improve this design to variable stiffness actuator (VSA), and perform various aerial physical interaction tasks together with a flying system using specific controllers for aerial physical interaction. In fact, the flexiblejoint arm proposed in the paper can be extended to a variable stiffness joint arm, by adding a linear motor and changing elastic components from linear springs to nonlinear ones, such as using rubbers.

Fig. 1 :

 1 Fig. 1: The light-weight flexible-joint arm.

Fig. 2 :

 2 Fig. 2: CAD Model of the flexible-joint arm.

Fig. 3 :

 3 Fig. 3: CAD Sketch of the flexible-joint arm.

Fig. 4 :

 4 Fig. 4: Schematic figure of the flexible-joint arm with notations used in the paper. The linear motor is grayed out and inactive for the current design. A future development includes activating this part too, for upgrading the current elastic actuator design to a variable stiffness actuator one.

Fig. 5 :

 5 Fig.5: Frequency response of motor for different step inputs. The input is desired motor velocity, and output is measured motor velocity. The dashed lines represent the conditions 'with load' (wL), and the solid lines the ones with 'no load' (nL). Different colors correspond to pulse trains of different frequencies (T1, T2, T3, and T4).

Fig. 6 :

 6 Fig. 6: Frequency response of the flexible-joint arm system for different step inputs. The input is the motor position and the output is the measured link position. Different colors present the trials with different step input profiles in motor velocities.

Fig. 9 :

 9 Fig.9: Flexible-joint arm is attached to the quadrotor system. Two experimets are designed; explosive arm movement (on the left) and aerial physical interaction (on the right). The second one is a task, where flexible-joint arm is sliding on a sloped surface, while quadrotor is performing a stable flight. The rigid link mechanism consists of a mass mc = 44gr at the tip of the link, in a distance of l = 18cm to the link center of gravity, where m l = 44gr.

Fig. 10 :

 10 Fig. 10: Explosive movement of the flexible arm during flight on the quadrotor. On the left picture arm swings right, and on the right picture it swings left. More detail is given in the video attached to this paper.

Fig. 11 :

 11 Quadrotor Positions

Fig. 12 :

 12 Fig. 12: Experiments for aerial physical interaction. Three figures represent the evolution of the motion along the -ŷ-axis, as mentioned in Fig. 11, where the results are shown. Three axes of Cartesian coordinates are shown in the middle figure in three different colors. Full experiment is available in the video attached to this paper.

TABLE I :

 I Estimated parameters of the flexible-joint arm

Preprint version, final version at http://ieeexplore.ieee.org/ 3 2015 IEEE ICRA

http://www.avagotech.com/docs/AV02-0188EN

http://arduino.cc/en/Main/arduinoBoardUno

https://www.mikrocontroller.com Preprint version, final version at http://ieeexplore.ieee.org/ 5 2015 IEEE ICRA

Preprint version, final version at http://ieeexplore.ieee.org/ 6 2015 IEEE ICRA

https://www.leapmotion.com/


Nicolas Staub

Antonio Franchi

Battery-aware Dynamical Modeling and Identification for the Total Thrust in Multi-rotor UAVs using only an Onboard Accelerometer

We propose and experimentally validate a new class of models for the total thrust generation in multi-rotor UAVs which is suitable for low-and middle-end platforms. Differently from typical models assuming to instantaneously control the rotor spinning velocity, in the proposed class we consider that the total thrust has its own dynamics and its final value explicitly depends both on the pseudo-setpoint commands given to the motor driver and the measurement of the battery terminal voltage. We compare the different model instances within the class using a principled experimental setup in which the total thrust is precisely measured using a motion capture system as ground truth, instead of relying on a setup based or noise-prone force sensors. We then show that the use of a dynamical model that includes also the battery terminal voltage significantly improves the prediction ability of the model in terms of accuracy. Finally we show how the proposed model can be identified using on-board only acceleration measurements, achieving a surprisingly good accuracy when compared with the ground truth case. We expect that the use of the proposed model will be important both in case of precise flight control and in the case of aerial physical interactive tasks.

I. INTRODUCTION

Unmanned Aerial Vehicles (UAVs) have been studied for decades, among those vehicles multi-rotor platforms raised recently the interest of the robotics research community for their agility, affordability and robust mechanical design. Typical multi-rotor platforms are actuated by a set of propellers with parallel rotational axes, whose number and spatial dispersion can vary among the platforms. This design is efficient in hovering, however it leads to underactuation, which must be addressed in the flight control.

Nowadays the aerial robotics community aims at enabling physical interaction for aerial systems like, e.g., for inspection [START_REF] Franchi | EU Collab[END_REF] or assembly [START_REF] Arcas | EU Collab[END_REF] tasks. To do so, control strategies have been proposed, based on different approaches. Nguyen et. Al. propose in [START_REF] Nguyen | Hybrid force/motion control and internal dynamics of quadrotors for tool operation[END_REF] a hybrid force/position control, with a tangential motion control and a normal force control.Yüksel et. Al. [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] propose an Interconnection and Damping Assignment Passivity Based Control (IDA-PBC) scheme, where the total energy of the system is reshaped to mimic a target dynamics. Some other methods using more classical admittance [START_REF] Augugliaro | Admittance control for physical human-quadrocopter interaction[END_REF] and impedance schemes [START_REF] Gioioso | Turning a near-hovering controlled quadrotor into a 3D force effector[END_REF] are also proposed in the literature. In this context the knowledge of the force exerted on the environment by the aerial robot is of paramount importance to achieve a precise interaction force control. Furthermore, it can be used also to retrieve external disturbance forces like, e.g., due to the wind. 1 CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France 2 Univ de Toulouse, UPS, LAAS, F-31400 Toulouse, France nicolas.staub@laas.fr,antonio.franchi@laas.fr

To know the interaction force, one could either measure it directly at the point of contact or estimate it. The measurement approach relies on force sensors, whose ratio performance over weight and compactness is not suitable to be embedded in multi-rotors in the low/mid-cost range. Moreover to measure any force acting on the system by this mean, the sensor should be located at the interaction point, which is highly restricting for interaction applications and not possible for wind estimation. This leads to our choice of investigating an estimator based approach.

This observer should require low computation, and must be easy to reproduce, those two requirements are set to enforce usability of this method. The models used in the estimation process should remain as simple as possible to limit computational cost but still capture all the essential dynamics of the system in order to produce a trustworthy estimate. Once the models are set, their parameters estimation should be highly reproducible to be performed quickly after any modifications of the systems.

Estimation algorithms have been already studied in the literature. They mainly fall within two categories: force mappings and model based observers. A static mapping between commanded thrust and actual force has been proposed by [START_REF] Bellens | A hybrid pose/wrench control framework for quadrotor helicopters[END_REF] achieving a ±1 Newton accuracy. A more accurate, yet static, mapping of the force produced by a typical brushless motor for aerial vehicle has been proposed by [START_REF] Spica | An open-source hardware/software architecture for quadrotor UAVs[END_REF], based on discretized force measurements for desired commands. In [START_REF] Ruggiero | Impedance control of VToL UAVs with a momentum-based external generalized forces estimator[END_REF] a momentum based external generalized force observer is presented, which requires to measure the whole dynamical system state, the control torques and the thrust to produce an estimate of the external forces and torques. The methods has been proven to work indoor by using a precise off-board motion capture tracking. Yüksel et. al [START_REF] Yüksel | A nonlinear force observer for quadrotors and application to physical interactive tasks[END_REF] proposed a non linear force observer, which requires to compute online both the position and attitude up to their 2nd derivatives. In [START_REF] Augugliaro | Admittance control for physical human-quadrocopter interaction[END_REF] a classical Unscented Kalman Filter is presented to estimate the external force and torque acting on a quadrotor. This approach has also been tested with the use of an off-board motion capture system.

The main contribution of our work is to present a new kind of approach to estimate the total thrust produced by the propellers which is fast, computationally inexpensive and can be implemented only using onboard measurements (a singleaxis accelerometer), and without the knowledge of the multirotor state and/or it's derivatives. Contrarily to other existing static approaches [START_REF] Bellens | A hybrid pose/wrench control framework for quadrotor helicopters[END_REF], [START_REF] Spica | An open-source hardware/software architecture for quadrotor UAVs[END_REF], [START_REF] Sa | System identification, estimation and control for a cost effective open-source quadcopter[END_REF] our approach is dynamical, since it takes into account the dynamics of the propeller in changing its rotational speed. In particular we propose a class of prediction models that use only the pseudo-setpoint of the propeller motor controller and the battery terminal voltage information, hereafter simply referred to as battery voltage. We conduct an experimental investigation to find the best and simplest model in this class using the acceleration retrieved from a motion capture system as ground-truth to validate our accelerometer-only-based approach. By comparing the prediction error for a set of possible inputs, we show that the addition of the battery voltage information in the model provides a manifest better force prediction. We alsot show experimentally that the usage of only onboard acceleration measurement for identification does not result in a significant degradation of the prediction when compared to the use of ground truth for the same purpose.

The rest of the paper is organized as follows. Section II presents our proposed class of models for the force generation. Those models are proposed based on observations made in the literature or on our test bed. In Section III the method to estimate the parameters of the aforementioned models are derived. Section IV describes the experimental set up that we used to acquire the data and Sec. V presents the results of the experiment analysis with respect to our claims. Finally Section VI presents a summary of the presented work and ongoing developments.

II. THE PROPOSED MODEL CLASS

In this section we extend the classical model of an underactuated multi rotor vehicle in order to include the presence of an unknown motor-driver dynamics for each propeller. We first start by recalling the standard model.

A. Multi-rotor Vehicle Dynamics and Standard Motor Model

We denote with {A} the inertial frame whose origin is O A and unit vectors along the axes are denoted by { a 1 , a 2 , a 3 }. The vector p ∈ R 3 represents the position of the Center of Mass (CoM) of the aerial vehicle expressed in {A}. Denote with {B} the body fixed frame whose origin coincides with the CoM of the vehicle, and axes w.r.t. {A} are represented by the unit vectors { b 1 , b 2 , b 3 }. We assume, as standard in the aeronautic field, that {B} is oriented with the North-East-Down (NED) convention, i.e., b 3 is pointing downwards when the aerial vehicle is hovering. The orientation of the vehicle in {A} is provided by the rotation matrix

A R B = R = ( b 1 b 2 b 3 ) [12].
Denote with N the number of rotors of the vehicle, which are all spinning on the plane defined by b 1 and b 2 , thus about b 3 . The equation of motion of the vehicle's CoM can be written as:

m p = mg a 3 + Rf, (1)
where m ∈ R + is the mass of the vehicle, g a 3 is the gravity acceleration, and f ∈ R 3 is the total force acting on the vehicle, expressed at the center of gravity in {B}. The dominant aerodynamics of each rotor i, for i = 1, . . . , N, produces a force (thrust) F i b 3 (see, e.g., [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]). In contactfree flight, the total force, expressed in body frame, is then

f = N ∑ i=1 F i b 3 + δ = F T b 3 + δ , (2)
where

F T = N ∑ i=1 F i (3)
is the total thrust, and δ comprises second order aerodynamic forces mainly due to flapping and drag effects, that are typically neglected in nominal working conditions [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF]. It can be assumed [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF] that the thrust produced by the rotor i is instantaneously related to its spinning velocity ϖ i by the following relation

F i = c F ϖ 2 i , (4)
where c F > 0 is an aerodynamic constant that depends on the specific properties of the propeller used.

B. Possible Drawbacks of the Standard Model

Many of the previous works, see e.g. [START_REF] Mahony | Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor[END_REF], assume that ϖ i is the control input of the system and that c F is known. Under these two assumptions the F i for i = 1 . . . N can be considered as control inputs, which makes the control problem simpler. However, these assumptions are hard to be met in reality for the reasons that are explained in the following. First of c F is in many cases unknown, depends on the propeller type, and even vary within the same type because of normal large-scale production variations. Second of all the motor has its own dynamics. Therefore ϖ i cannot be changed instantaneously acting on the torque applied by the motor driver.

In fact, in the best case, the motor control input is usually a setpoint for ϖ i , which is tracked by the motor driver with a certain dynamics and accuracy. The control accuracy of F i depends then on the accuracy in the knowledge of the parameter c F (needed in order to generate the setpoint for the motor driver) and on the accuracy of the driver in tracking the velocity setpoint. For typical multi-rotor platforms, especially in the low/mid-cost range it is not possible to send a setpoint for the spinning velocity to the motor drivers, but rather a pseudo-setpoint, e.g., a PWM signal, which is monotically related to the steady state spinning velocity, if all the other flight conditions are constant.

In both cases (for setpoint and for pseudo-setpoint) the battery voltage level has an impact on the behavior of the motor controller, as noted, e.g., in [START_REF] Sa | System identification, estimation and control for a cost effective open-source quadcopter[END_REF]. In particular, the battery voltage decreases as long as the battery discharges due to the increase of the internal resistance. This measurement can be found among all the type of platforms. During contactfree flight, the dependence to the battery voltage can be overcome by using an adaptive term in the control law which compensates the discharge of the battery. Such a method can take the form of a mass estimator or an integral term [START_REF] Spica | An open-source hardware/software architecture for quadrotor UAVs[END_REF], [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor control[END_REF], which increases the average of the commands sent to the motors, thus compensating the voltage drop. This approach is working because the only force acting on the system are the weight and the thrust generated by the propellers. For physical interaction this technique can not be used since a force sensor is unavailable in our assumptions and the interaction force involved in the balance is typically unknow.

In order to have a good control of the force exerted by the system it is instead important to have control on the total thrust F T exerted by the rotors. Therefore another viable approach, considered in this paper, is to take into account the battery voltage influence and the motor dynamics directly in the input-output nonlinear model. Since we are interested in controlling the force we try to identify the direct relation between the (pseudo)-setpoints and the battery voltage as inputs and the force as output.

F (u i , β) u 1 , .., u N β 1 τ 1 s + -

C. Model Based on (Pseudo)-Setpoint and Battery Level

Denote with u i ∈ [u min , u max ], with u min > 0 the control input of the motor driver, which may represent either a setpoint or a pseudo-setpoint for the driver of the i motor of the vehicle. Denote with β the battery voltage level, that we assume to be measured. Our goal is then to propose a modeling of the relation between u i and β and F i that is as simple as possible but captures all the relevant dynamics.

Based on the analysis of experimental data, see e.g., [START_REF] Spica | An open-source hardware/software architecture for quadrotor UAVs[END_REF], [START_REF] Sa | System identification, estimation and control for a cost effective open-source quadcopter[END_REF], we have made the following observations: O 1. if u i is kept constant for a time window of the order of magnitude of ≈ 1 second then F i reaches a constant value that is is monotonically increasing w.r.t. u i and monotonically decreasing w.r.t. β .

Motivated by the previous observation we propose the following simple model for the force

Ḟi = 1

τ i F i (u i , β ) -F i (5)
where F(u i , β ) is an unknown nonlinear map, and τ i is an unknown time constant, this model can be represented as a nonlinear block and a linear first order system, see Fig. 1.

In order to provide a simple expression of F T let us consider the following reasonable assumptions: A 1. F i (u i , β ) and τ i are the same across all motor controllers, i.e., for i = 1 . . . N A 2. F i (u i , β ) is a smooth function that can be well approximated by a finite polynomial of a suitable order in the region of interest u i ∈ [u min , u max ] and β ∈ [β min , β max ] Given Assumptions 1 and 2 we write

τ i = τ, ∀ i = 1 . . . N (6) F i (u i , β ) = F(u i , β ) ≈ n u ∑ j=0 n β ∑ k=0 (α jk u j i β k ), ∀ i = 1 . . . N (7)
where n u , n β ∈ N have to be chosen taking into account the desired model complexity, and α jk , with j = 1 . . . n u , k = 1 . . . n β , are n u • n β parameters to be estimated. Using (3),( 5),( 6), [START_REF] Bellens | A hybrid pose/wrench control framework for quadrotor helicopters[END_REF] we can write the dynamics of the total thrust F T as

ḞT ≈ N ∑ i=1 1 τ i n u ∑ j=0 n β ∑ k=0 (α jk u j i β k ) -F i (8) = 1 τ n u ∑ j=0 n β ∑ k=0 (α jk β k N ∑ i=1 u j i ) -F T (9) = 1 τ n u ∑ j=0 n β ∑ k=0 α jk v jk -F T (10)
where we compactly wrote v jk = β k ∑ N i=1 u j i . The equation ( 10) is describing a class of sufficiently simple nonlinear dynamical models for the total thrust of a multi-rotor platform using only the (pseudo)-setpoint and the battery voltage, both information that can bee found on nearly every, if not all, platforms. This class has been introduced in order to illustrate two of our assertions which constitute part of our contribution: Claim 1. Using only (pseudo)-setpoint and battery voltage in the model, it is possible to predict the actual total thrust with a 'good' confidence.

Claim 2. Furthermore the use of the battery voltage information improves 'substantially' the quality of the prediction.

The model class size is only limited by the choice of n u and n β , among the models present in the class some have better performance than others. Our goal is to find the simplest model (i.e., the one with the lowest number of parameters) that provides a sufficient prediction performance compared to more complex models.

III. IDENTIFICATION OF MODEL PARAMETERS

Model [START_REF] Yüksel | A nonlinear force observer for quadrotors and application to physical interactive tasks[END_REF] represents a class of models, depending on the values of the parameters n u , n β . In this class of models, three quantities are needed in order to estimate the parameters τ and the α jk 's: the model inputs, i.e., the (pseudo)-setpoint u i and the battery voltage β , and the model's output, i.e., the total thrust F T . The information on F T can be replaced by acceleration information, given that the mass m of the system is known. In practical situations two cases, detailed in the following, can happen.

The case in which p can be precisely measured (e.g., using a motion capture system). In this case the thrust force in the inertial frame, A f = Rf, can be computed by directly employing [START_REF] Franchi | EU Collab[END_REF]. The total thrust F T is then computed as A f . We shall use this case as ground truth in the experiments.

The second case is when the vehicle is equipped with a calibrated IMU which measures the proper acceleration of the vehicle expressed in the body frame {B}, i.e.,

ā = R T ( p -g z) + η IMU , (11)
where η IMU is some additive noise with zero mean. This case will be used in the experiments to show that the model can be effectively identified using only onboard measurements. Using (1),( 2) in [START_REF] Sa | System identification, estimation and control for a cost effective open-source quadcopter[END_REF] we obtain In typical conditions (i.e., at low speed and at a certain distance from the ground ) F T is much larger than the third component of δ . Therefore, we can write

ā = F T m z + 1 m δ + η IMU . (12)
F T = m z T ā + η F ≈ m āz . ( 13
)
Where η f represents a negligible contributions of both the aerodynamic effects and the IMU noise. Assuming that the mass m of the vehicle can be measured before the flight, we shall use [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor control[END_REF] as an onboard measurement of the total thrust. In particular, we make the following claim that will be tested in the experimental section:

Claim 3. The readings of calibrated low/middle-cost accelerometers are sufficiently accurate to perform a reliable identification of the parameters of the proposed model, leading to practically acceptable onboard total thrust prediction.

A. Use of Estimated Model in Contact Tasks

Notice that ( 13) is valid only in contact-free flight, i.e., when the non-gravitational forces acting on the vehicle can be expressed by [START_REF] Arcas | EU Collab[END_REF]. During contact, (2) includes also the interaction forces f e , thus it can be rewritten as

m p = mg a 3 + Rf + f e . ( 14
)
The estimation of F T can be performed in the same way as before if f e is known.

If instead the estimation of τ and F is performed using measurements taken during contact-free flight phases, the identified model can then be used to estimate f e in contact phases since both τ and F are constant over time. We leave the validation of this claim as a future work. while in the following we provide the method to estimate the parameters of the force model.

IV. EXPERIMENT DESIGN

In order to estimate the parameters in [START_REF] Yüksel | A nonlinear force observer for quadrotors and application to physical interactive tasks[END_REF] we designed an experimental set up based on the Telekyb framework [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor control[END_REF] and composed by a motion capture system (mocap), a control computer and a quadrotor, see Fig. 2. The mocap provides the tracked object position at 100Hz with millimeter-scale accuracy. A very accurate estimation of the velocity and acceleration is obtained post-processing the measured position with a non-causal Savitzky-Golay filter [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF].

The quadrotor mechatronics is based on the Mikrokopter platform, whereas the flight control software has been replaced by the Telekyb [START_REF] Grabe | The TeleKyb framework for a modular and extendible ROS-based quadrotor control[END_REF] one. Two Xbee radio transmitters are used to send commands to the quadrotor and to record the telemetry data, respectively. The data recorded are: the battery voltage, a timer for synchronization, the four motor pseudo-setpoints (PWM signals) sent by the flight controller to the brushless controllers and the onboard acceleration measurements provided by the threee 3D LIS344alh accelerometers (0.0057g 0 m/s 2 resolution and ±2g 0 m/s 2 range). The accelerometers are calibrated following the procedure detailed in [START_REF] Spica | An open-source hardware/software architecture for quadrotor UAVs[END_REF]. The control computer runs the Telekyb control framework under the middleware ROS. The flight controller used for the quadrotor is a near-hovering scheme which allows to follow trajectories where the roll and pitch of the quadrotor remain within ±30 • .

The thrust is then computed twice using (i) the acceleration measurement from the motion capture (ground truth) and (ii) the accelerometer reading, as explained in Sec. III.

Any trajectory that spans the the battery range and the pseudo-setpoint range of interest can be used for the identification. The former requirement asks for a flight duration that discharges the battery enough and the latter one asks for a sufficiently rich acceleration content of the tracked trajectory. To meet the arena-size constraint as well, a vertical trajectory has been chosen providing enough space to reach high accelerations while remaining near the hovering attitude. To ensure the spectral richness of the acceleration the vertical trajectory is composed of five sinusoids with different pulsations:

z(t) = ρ + 5 ∑ i=1 a i sin(ω i t + φ i ), (15)
where ρ = 2 is a position offset and

[ω 1 , . . . , ω 5 ] = 2π • [0.3, 0.2, 0.1, 0.5, 0.6] [a 1 , . . . , a 5 ] = 0.3 • [ω -2 1 , 0.95ω -2 2 , 0.9ω -2 3 , 0.8ω -2 4 , 0.7ω -2 5 ] φ = π • 1 2 , 1 3 , 1 4 , 1 5 , 1 6
are the pulsations, the amplitudes, and the phase shifts of the sinusoids, respectively. With this choice the input is then persistently exciting with order 10, making us able to identify of up to 10 parameters in [START_REF] Yüksel | A nonlinear force observer for quadrotors and application to physical interactive tasks[END_REF], see, e.g., [START_REF] Ljung | System Identification -Theory For the User[END_REF].

V. EXPERIMENTAL RESULTS

To estimated the parameters of (10) we used a predictive error method coupled with a grey-box model. This method is a numerical optimization with a cost function based on the norm of the prediction [START_REF] Soderstrom | System Identification, ser. Prentice Hall International Series In Systems And Control Engineering[END_REF]. Three different analysis have been conducted on the data and are described hereafter.

In order to illustrate the experimental setup, the video attached to the paper shows a preliminary experiment with a simpler trajectory.

ID of the model 1 2 3 4 5 6 7 8 9 10

n u 0 0 1 1 0 2 1 2 2 2 n β 0 1 0 1 2 0 2 1 2 2
Numb. of param. prediction error (N) Fig. 3: Boxplot of the prediction error for the proposed models (see Tab. I), the parameters are estimated with the ground truth and the prediction is compared against the ground truth.

A. Model Order Choice and Mocap-based Identification

The number of parameters in ( 10) is (n u + 1) • (n β + 1) + 1, thus it depends on the chosen value of n u and n β . Our goal is to find n u and n β such that the system dynamics is well described by the estimated model but with the smallest number of parameters.

As explained in Sec. IV the maximum number of parameters that we can identify with our experimental setup is 10. Therefore we can consider n u , n β ∈ {0, 1, 2}. The 10 resulting models are summarized in Table I, where models 9 and 10 have the same structure but the parameters are estimated with a different initial guess as explained in the following. We discard model 1 since the system dynamics can not be described just by a constant. Models 2 and 5 are also discarded as they do not contain information on the (pseudo)-setpoint. The initial guess for the estimation of the parameters of more complex models is provided by the solution of the antecedent in the model class, thus creating the following orders: 3 → 4 → 7 → 9 and 3 → 6 → 8 → 10.

We first estimated the parameters using a set of data from the mocap (we shall present the estimation using onboardonly sensor in Sec. V-C). The estimated parameters are then validated against another set of data from the mocap. To better validate the ability of prediction of our models w.r.t. the battery effect, the battery voltages in the 2 sets span completely different values. Boxplots graphs of the absolute value from the prediction error are presented in Fig. 3, the predictions are quite good overall, with a prediction error median around or under 0.1 N, which is remarkable considering that the range of the recorded force during the experiment is [7.5 N, 12.2 N]. This demonstrates that the proposed class of models is able to predict the total thrust accurately thus validating our first claim. The complete set of estimated parameters is presented in Table II. Notice that the time constant is of the same magnitude for all the models, with model 9 having a slightly higher one.

B. Discussion on Battery Influence

From the previous model perfomances one can investigate the importance of the battery. Models 3 and 6, where battery voltage information is not used, are both outperformed (higher median, wider dispersion) by the models containing both pseudo-setpoint and battery information, i.e., 4, 7, 8, 9, and 10. As models 3 and 4 have the same number of inputs the use of battery data is the only changing factor that can explain the better fitting of 4 with respect to 3. Moreover, we can notice that despite of the fact that model 6 has more inputs than model 4, its prediction ability is worse. In general we can conclude that the models in which n β = 0 can better predict the system behavior by taking into account the battery voltage drop along the flight. This fact removes, e.g., the need for an adaptation term in the flight controller. This result validate our second claim, i.e., that the use of the battery voltage information improves 'substantially' the quality of the prediction. A direct comparison between model 4 and 6 on a chunk of the validation set is shown in Fig. 4.

C. Accelerometer-based Identification

An analysis was conducted to determine the validity of parameters estimation based only on the onboard accelerometer measurements. Then the quality of the prediction of those models has been compared to ground truth (i.e., the validation set used in the mocap case). The prediction error against the mocap is presented in Fig. 5.

The comparison of Fig. 3 and Fig. 5 shows that the parameter estimation process can also be conducted without a mocap at the cost of minimal variation of the prediction error, allowing the proposed force prediction framework to be deployed without a mocap. This result validates our third claim, i.e., that the use of a calibrated low/middle-cost accelerometer is sufficient to perform a reliable identification of the parameters of the proposed model, leading to practically acceptable onboard total thrust prediction. The complete set of estimated parameters can be found in Table III.

D. Choice of the Best Model

From our analysis, model 4 results to be the best compromise between accuracy and complexity. In fact, this model resulted able to describe the dynamics of the system with the same level of prediction ability of models with more parameters. Furthermore, model 4 is computationally lightweight and therefore can be easily implemented, e.g., on a microcontroller. Furthermore, this model keeps substantially the same the prediction ability if the estimation of the parameters is conducted using onboard accelerometer data.

VI. CONCLUSIONS

In this work we proposed a class of dynamic models to predict the force (total thrust) generated by a underactuated multi-rotor system. The major contributions have been to consider the motor dynamics, include battery voltage information in the prediction model, and to identify the model parameters using only onboard accelerometer measurements. We conducted an experimental investigation to find the best model among the simplest ones in this class using the acceleration retrieved from mocap measurements as groundtruth.

By comparing the prediction error for a set of possible inputs, we have shown that the addition of the battery voltage information in the model provides a manifest better force prediction. From our experimental result it was also clear that the usage of only onboard acceleration measurement for identification does not result in a significant degradation of the prediction when compared to the use of ground truth for the same purpose.

Nevertheless the assumption is made that the only force acting on the system during the record of the data for parameter estimation are the weight force and the total thrust, thus an outdoor estimation of the parameters in windy conditions is not possible. However, the proposed models can be used for an outdoor estimation of external forces acting on the system, like e.g., wind.

In the near future we are working toward the integration of this work in a physical interaction framework, where we want to use the force prediction to close the control loop on external forces.

Fig. 1 :

 1 Fig. 1: Block diagram of the proposed model

Fig. 2 :

 2 Fig. 2: The multi-rotor platform used at LAAS

2 3 3 4 4 5 7 7 10 10

 10

Fig. 4 :

 4 Fig. 4: Comparison of two estimated model one not including battery info (model 6) and one including battery info (model 4). Both models have the same number of parameters, the one containing battery information is clearly more performant.

Fig. 5 :

 5 Fig.5: Boxplot of the prediction error for the proposed models (see Tab. I), the parameters are estimated with onboard acceleration and the prediction is compared against the ground truth

TABLE I :

 I Different models considered in the identification.

	0.5
	0.4
	0.3
	0.2
	0.1
	0	3	4	6	7	8	9	10
							org/	4	2015 IEEE ICRA

Preprint version, final version at http://ieeexplore.ieee.

TABLE II :

 II Recapitulative table of the estimated parameters of the models when the estimation is conducted with mocap data.

TABLE III :

 III Recapitulative table of the estimated parameters of the models when the estimation is conducted with accelerometers data.

Preprint version, final version at http://ieeexplore.ieee.org/ 2

Preprint version, final version at http://ieeexplore.ieee.org/ 3

Preprint version, final version at http://ieeexplore.ieee.org/ 5

Preprint version, final version at http://ieeexplore.ieee.org/ 6 2015 IEEE ICRA


Matteo Gagliardi
email: matteo.gagliardi@tuebingen.mpg.de.

Giuseppe Oriolo
email: oriolo@dis.uniroma1.it

Heinrich H Bülthoff

Antonio Franchi
email: antonio.franchi@laas.fr

Image-based road network clearing without localization and without maps using a team of UAVs

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The use of coordinated teams of robots allows to require less from each hardware device and at the same time achieve higher performance and robustness. For such reasons the design of controllers for multiple robots gained great relevance in the literature, see, e.g., [START_REF] Howard | Experiments with a large heterogeneous mobile robot team: Exploration, mapping, deployment and detection[END_REF], [START_REF] Pasqualetti | On cooperative patrolling: Optimal trajectories, complexity analysis, and approximation algorithms[END_REF] and references therein. Particular interest has been paid to the exploration of a given region [START_REF] Howard | Experiments with a large heterogeneous mobile robot team: Exploration, mapping, deployment and detection[END_REF], [START_REF] Franchi | The sensor-based random graph method for cooperative robot exploration[END_REF], [START_REF] Durham | Distributed pursuit-evasion without global localization via local frontiers[END_REF], allowing the execution of several subtasks: guaranteed clearing for pursuit-evasion [START_REF] Durham | Distributed pursuit-evasion without global localization via local frontiers[END_REF], [START_REF] Gerkey | Visibility-based pursuitevasion with limited field of view[END_REF], models for forest area coverage [START_REF] Zheng | Multi-robot forest coverage[END_REF], autonomous or semiautonomous control for rescue purposes [START_REF] Doroodgar | A hierarchical reinforcement learning based control architecture for semi-autonomous rescue robots in cluttered environments[END_REF], mine detection [START_REF] Kumar | Cognitive maps in swarm robots for the mine detection application[END_REF], exploration models for underwater regions [START_REF] Rathnam | Distributed communicative exploration under underwater communication constraints[END_REF], just to mention a few. Along this line, in this work we present an algorithm for the guaranteed clearing of a network of roads using a team of UAVs.

The clearing problem is well established in the research on robotics: so far, however, there have been still few results suitable for real-world applications. In fact, a large fraction of such research work has developed models using assumptions far from those suitable for realistic applications, e.g., sensors with infinite range, a priori knowledge of the environment or target model, or perfect localization. Furthermore, From a theoretical point of view, important results have been achieved in this field using three main approaches: graphbased, visibility-based and sweeping research. In graphbased approaches the team is typically supllied with a complete or partial knowledge of the environment to be explored [START_REF] Rocha | Decision methods for distributed multirobot patrol[END_REF]- [START_REF]Multi-robot surveillance: an improved algorithm for the GRAPH-CLEAR problem[END_REF]. In [START_REF] Kolling | Pursuit-evasion on trees by robot teams[END_REF] the authors improve their previous works by proposing a clearing strategy that does not require a map of the environment. Using visibilitybased approaches one can cope with the lack of a priori information about the environment, as, e.g., in [START_REF] De Hoog | Role-based autonomous multi-robot exploration[END_REF] where a kind of communication chain through limited-range sensors with a perfect localization in the environment is created in order to report all the data acquired to a command center. An important contribution to the visibility-based approach is proposed in [START_REF] Gerkey | Visibility-based pursuitevasion with limited field of view[END_REF] in which the clearing is performed by terrestrial agents using sensors with limited range capabilities. Solutions using Simultaneous Localization and Mapping (SLAM) make possible to avoid a priori assumptions about the environment, to the detriment of the memory space and computational load. Another relevant work is [START_REF] Carpin | Variable resolution search with quadrotors: Theory and practice[END_REF], where the searching problem is addressed using a variable resolution data structure and field experiments are also presented. In [START_REF] Durham | Distributed pursuit-evasion without global localization via local frontiers[END_REF], [START_REF] Durham | Distributed pursuit-evasion with limited-visibility sensor via frontier-based exploration[END_REF] all the previously adopted assumptions are removed at once: each agent uses limited-range planar sensors and shares information about its own local exploration frontiers with its neighbors. Almost all the aforementioned works developed algorithms that are suited for ground robots. However, using UAVs instead of ground robots makes possible the inspection, monitoring and search of large areas that may be inaccessible by terrestrial vehicles. Multi-rotor UAVs have, in particular, gained a lot of attention in the community [START_REF] Michael | Opportunities and challenges with autonomous micro aerial vehicles[END_REF], [START_REF] Franchi | Shared control: Balancing autonomy and human assistance with a group of quadrotor UAVs[END_REF]. Along this line, [START_REF] Kleiner | Guaranteed search with large teams of unmanned aerial vehicles[END_REF] presents an algorithm to search for ground targets using multiple UAVs employing the technique of sweep lines.

Here, we propose a novel vision-based outdoor clearing algorithm for a group of UAVs, equipped with limited-fieldof-view down-facing cameras, that relies neither on global localization nor on a mapping of the area to be cleared. The method call for each UAV to store only the current image and the previous one in order to work properly, thus requiring limited computational capabilities and memory space, a feature that is mandatory, e.g., in the so called nano-UAVs.

The paper is organized as follows. In Sec. II a formal definition of the scenario is provided. A theoretical explanation of the clearing strategy is exposed in Sec. III. Sec. IV describes the image processing part of the algorithm and Sec. V presents the control architecture adopted for the coordination of the whole team. Finally, Sec. VI shows simulation results with a group of 8 quadrotor UAVs and Sec. VII concludes the paper with some final remarks. We assume that the aerial occlusions caused by the buildings can be neglected, i.e., the ground environment, as seen from the UAVs, is 'sufficiently flat'.

Each UAV is equipped with a rigidly attached downfacing monocular camera having finite resolution L × L and both the vertical and horizontal field of view equal to α ∈ (0, π). An NED (North-East-Down) frame, denoted with R i : O ix i , y i , z i , is attached to the camera of each UAV U i . The origin O i is located at the focal center of the camera and the z i axis is directed along the focal axis. The position of O i in W is denoted by p i = (p x i , p y i , p z i ) ∈ R 3 and R W i denotes the rotation matrix representing the orientation of R i with respect to W . We assume that each UAV can store only its current and previous image (see the discussion in Sec. I).

The UAV system is also endowed with a decentralized communication infrastructure that provides two different channels. A High-Bandwidth Short-Range (HBSR) channel is used to transmit the acquired image to every robot that is closer than D HBSR max . A Low-Bandwidth Long-Range (LBLR) channel is used to transmit brief messages to every robot that is closer than D LBLR max D HBSR max . Each UAV is equipped with a flight controller that guarantees the following conditions during the whole mission:

• altitude of the camera is kept at a given constant desired value pz , i.e., p z i = pz ; • the focal axis of the camera z i is kept parallel to z W ;

• the x i axis is kept parallel to x W ;

• the horizontal velocity ( ṗx i , ṗy i ) is regulated to a value u i ∈ R 2 that is assigned by the clearing algorithm described in the following section. The four requirements on the flight controller can be easily achieved using, e.g., a barometer for the regulation of the altitude, the IMU measurements and a magnetometer for the regulation of z i and x i , respectively, and a velocity sensor (e.g., based on the optical flow plus the barometer to recover the scale) for the regulation of the horizontal speed.

Being the image acquisition a discrete process, we consider a discretized time, denoting with k a generic time step and with I k i the image acquired by the i-th UAV camera at k. We assume that a specific algorithm is available to extract the road in I k i that is designed and tuned for the environment at hand. Such algorithm is able to extract the roads in I k i provided that the road width, expressed in pixels, is always contained in [W L min ,W L max ], where W L min > 0 and W L max < L. In order to obtain a well-posed problem setting, it is needed that:

• the altitude, the image resolution, and the road width allow a reliable segmentation, i.e.:

W min 2 pz tan (α/2) ≥ W L min L , W max 2 pz tan (α/2) ≤ W L max L , (1)
where W min and W max are maximum and minimum width of the roads in A; • D LBLR max is larger than the diameter of A, i.e., the maximum distance between any two points of A. The first assumption guarantees a successful road extraction while the second one guarantees that two UAVs can always communicate using the LBLR channel. Finally we assume that the lines where the road network crosses ∂ A apart from the starting line S are constantly patrolled by some guards. In the case that S is the only access to A this assumption is not needed.

Notice the following facts:

• the set R is unknown (no map);

• each UAV can use only the current and previous image;

• each UAV ignores both its latitude and its longitude (no localization). We can now formally define our objective.

Clearing Problem: An unknown number of moving agents (e.g., vehicles) travels on the roads of R. Their motion is continuous but their initial position is unknown and their speed is unpredictable and unbounded. An agent is assumed to be detected only if it falls in the field of view of a camera or if it is spotted by a guard (i.e., if it crosses ∂ A).

The goal for the UAV team is to move and deploy over the area A at constant altitude pz so as to guarantee that any moving agent that starts moving in R S at t 0 is, in a finite time, detected. The UAVs enter the scene from the same starting line, and in a way such that S is covered by the camera of at least one of them.

III. CLEARING STRATEGY A. Preliminary Results

In order to precisely describe our clearing strategy we start by formalizing the problem from a set theory point of view. Define as U k C ⊂ U the subset of UAVs whose images at time step k are actually used for the clearing task in a way to be described below. We call these UAVs cleaners. V spans both the north and the east sides. However, by prediction it is clear that such a movement will not ensure F k V ⊂ Īk+1 , and it is therefore not executed. Fig. 2 introducing U k C is that at time step k not all the UAVs provide an image which is used to clean the environment. Denote with Īk the union of all the points of the plane < x w , y w > that are visible by the cleaners at time step k, i.e., Īk = ∪ i∈U k C Īk i . Consider a ground agent, and denote with T k its trail, i.e., the area swept by the agent during its motion from time step k -1 to time step k. The ground agent is detected at k if

T k ∩ ∂ A ∪ Īk = / 0, (2)
i.e., if T k either crosses ∂ A or is visible (either entirely or partially) by one of the cleaners at time step k.

A region R S ⊂ R S is cleared at time step k if each agent in that region has been detected, i.e., if for every agent

T k ∩ R S = / 0 ⇒ ∃l ≤ k : T l ∩ ∂ A ∪ Īl = / 0,
The clearing problem can then be formulated as clearing the whole set R S .

Instrumental to the description of the clearing strategy is the introduction of V k , the connected visited subset of R S , defined as

V k = Conn S, R S ∩ ∪ l=1...k Īl (3)
i.e., the maximal connected subset of R S that includes S and has been covered by the cleaners up to time k. We define the frontier of V k as

F k V = ∂V k \ ∂ R S = ∂V k ∩ R • S , (4)
i.e., the points of the boundary of V k that are not contained in the boundary of R S (i.e., are part of the interior of R S ). Figure 2a depicts the sets defined so far in a typical case. Notice that the symbol • denotes the interior of a set. The entire F k V is in general the union of several connected components, i.e.,

F k V = ∪ j∈F k V F k V j , where each component F k V j
is referred to as a sub-frontier of F k V . The following result (based on the "pushing back frontiers" concept [START_REF] Durham | Distributed pursuit-evasion without global localization via local frontiers[END_REF], [START_REF] Durham | Distributed pursuit-evasion with limited-visibility sensor via frontier-based exploration[END_REF]) is the main principle around which our strategy is built.

Proposition 1. If V k is cleared and F k V ⊂ Īk+1 (5) then V k+1 is also cleared and F k+1 V ⊂ ∂ Īk+1 . If, additionally, F k+1 V = / 0, then R S is cleared at time k + 1.
The formal proof of the previous proposition is omitted here for brevity. Proposition 1 implies the following results: Corollary 1. Any strategy ensuring that 1) condition (5) holds at every k ≥ 1, and 2) there exists a finite K such that V K = R S , solves the clearing problem.

The second property of Corollary 1 implies that if

F 1 K = / 0 then F k V ⊂ Īk+1 • (6)
holds a number of times m > 1.

B. Image-based Clearing Strategy

Assume that by suitable image processing we are able to identify the portion of V k in each Īk and its frontier F k V . We shall give more details on how to implement such image processing in Section IV. Furthermore, in order to ease the description of the clearing algorithm, we temporarily assume that the commanded motion of the UAVs is perfectly executed. In Section IV we will show how this last assumption can be relaxed.

Our clearing method is based on Corollary 1. Assume that each cleaner can either stay still or move, at constant speed, along the 8 cardinal and ordinal directions, thus obtaining 9 motion primitives on the horizontal plane. Denote with D the set of 9 possible horizontal displacements that are obtained staying still or moving in one of the possible directions between two consecutive instants, i.e., D = {∆ 0 , ∆ N , ∆ NE , ∆ E , ∆ SE , ∆ S , ∆ SW , ∆ W , ∆ NW } where ∆ 0 denotes the null displacement.

At time step k = 1 the group starts with U 1 C = {U 1 } and U 1 located on top of the starting line S. Let D 1 1 ⊂ D be the set of all the displacements of U 1 such that (5) holds at time step k = 1. At the same time, it is easy to see that the following properties hold

S ∈ Ī1 1 , V 1 ⊂ Ī1 1 , F 1 V ⊂ Ī1 1 . Proposition 2. If U 1 moves in D 1 \ {∆ 0 } then F 1 V ⊂ Ī2 • (7)
A motion primitive candidate ∆ ∈ D is selected taking into account the location of F 1

V on the 4 sides composing ∂ Ī1 as follows:

• if F 1
V is entirely contained in one of the sides, the motion candidate is the displacement in the cardinal direction of the side itself, e.g., if F 1

V is entirely contained in the North side of the image then the candidate will be ∆ N (see Fig. 2a); direction, e.g., if F 1 V is contained in the West, North, and East sides then the candidate will be ∆ N ;

• if F 1 V spans
• if F 1
V spans two contiguous sides, the motion candidate is the displacement toward the ordinal direction between the two sides; e.g,, if the F 1

V is contained in the North and East sides the candidate will be ∆ NE (see Fig. 2a);

• if F 1
V spans only two non-contiguous sides or all four sides of the image, the candidate is ∆ 0 . In order to check the feasibility of the motion candidate U 1 performs a prediction of the location of F 1 V in Ī p(1) after the execution of the motion candidate. The motion candidate is feasible if it does not violate [START_REF] Gerkey | Visibility-based pursuitevasion with limited field of view[END_REF]. In case of feasibility U 1 actually executes the candidate primitive.

If the motion candidate is not feasible (see Fig. 2b) then the null motion primitive ∆ 0 is selected. Such situation implies that one or more UAVs have to be included in U 2 C to continue the clearing. To this end, the frontier F 1

V is shrunk and partitioned in sub-frontiers. Then U 1 assigns some of the sub-frontiers to other UAVs until its motion candidate results feasible. The UAVs who receive an assignment are then promoted to cleaners and included in U 2 C . This event is referred to as splitting. A splitting decision can be taken at each step k by any cleaner.

In the case, referred to as an encounter, where two or more cleaners get so close to each other that their images overlap, each cleaner deletes from its frontier the portions that are contained in the others' images. Figure 3 illlustrates this case.

Finally, if at a certain step the frontier of a UAV becomes empty, its clearing phase is temporarily over and it becomes a helper. The behavior of the helpers is described in Sec. V.

IV. IMAGE PROCESSING

In Section III the basic clearing strategy has been described under the assumption that an image processing scheme is available for reconstructing the set F k V at each step k. In this section we want to justify this premise by outlining a suitable algorithm that works also if a certain level of noise is present in the segmentation process.

Consider the segmented image Īk i and denote with ∂ Īk i the boundary of Īk i formed by the four lines of pixels on its four sides. Denote with s k 1 the segmented points in Ī1 1 and define ∂

F V k F V k F V k I k - I k -- I k - ∂ -ε I k -- i k - i s k i s ∂I k -- i ∂ i i i (a)
s k i = s k 1 ∩ ∂ Ī1 1 .
Consider now the case k = 1, where only the current image Ī1

1 is available. However, since by assumption S belongs to Ī1

1 we have that F 1 V can be computed as

F 1 V = Conn S, s k 1 ∩ ∂ s k 1 . (8)
Consider now the case k > 1. After a registration of the two images Īk i and Īk-1 i , it is possible to obtain a compound image that includes both, see Fig. 4a. Then consider a 'shrunk' version of ∂ Īk-1 i that is at distance ε from ∂ Īk-1 i and denote this internal stripe with ∂ -ε Īk-1 i . The new frontier can then be computed as

F k V = Conn F k-1 V , (s k i ∪ s k-1 i ) ∩ ∂ -ε Īk-1 i ∩ ∂ s k i . ( 9
)
Loosely speaking, ∂ -ε Īk-1 i represents a virtual barrier that avoids a backward propagation of the frontier computation.

Regarding the validation of candidate primitives, this can be easily done by checking if F k V lies in the interior of a frame displaced with the nominal commands. Furthermore, in order to deal with noise on the motion velocity input, which would prevent the aforementioned algorithm to work, the frontier is evaluated on a shrunk virtual boundary of the image, see Fig. 4b. In this way, even if the actual motion is off by the shrinking amount the frontier will still lie in the next image, thus allowing a correct execution of the algorithm. Such a solution requires a trade-off between two opposite demands:

• amount of shrinking as small as possible to make the clearing process faster; • amount of shrinking as large as possible to achieve robustness to imprecise motions of the camera (i.e., of the UAV).

V. FORMATION CONTROL

Dealing with a network of many roads, the usage of a single robot would not be able in general to clear all the paths belonging to the network. As already mentioned, during the splitting event, each cleaner can entrust a subset of F V to another agent that automatically joins the set U C . The objective of the formation control is to provide, whenever possible, each cleaner with a sufficient number of nearby agents, not belonging to U C , with which is possible to Fig. 5: Block scheme of the main operative steps of our algorithm. Starting from the left side: UAVs (input: velocities and data coming from other element of the team; output: current images); image-processing block (output: selected displacement for the cleaners and relative positions for the followers); formation control (A= splitting event; B=end of clearing; C=integration of helpers; output -desired velocities); the motion controllers represent the final step of the system. share the frontier if necessary. This feature requires the implementation of a collaboration strategy among the agents.

Our approach is to organize the team into a variable number of formations whose agents perform specific roles. In general each formation is composed by: one cleaner, up to two wing agents denoted with U W (called left and right wing) and up to n -3 tail agents, denoted with U T . The formation configuration chosen is shown in Fig. 6. An additional role is the helper: the number of helpers can be up to n and they are not part of a formation. Such roles are not fixed for each agent but change during the clearing process. The above formation choice ensures that all communication required is feasible by only using visual odometry and the HBSR telecommunication system.

During the clearing process each agent, except the cleaner, has an agent, among the members of its own formation, called target, with respect to which it has to maintain a certain relative position. In particular:

• Cleaner: it does not have any target and it is controlled directly by the motion strategy of Sec. III-B, that uses the image processing of Sec. IV; • Wings: their target is the cleaner (dashed horizontal red connections in Fig. 6). During the normal clearing phase they have to maintain a position such that their images always overlap two opposite halves of the image of their own target. The wings are the only agents that can be called to share, if necessary, a portion of F V of their cleaner during a splitting event.

• Tails: the tail agents are arranged as shown in Fig. 6.

The target of each tail agent is the previous tail agent (dashed green connections), except for the first whose target is the cleaner (dashed vertical red connection). After a splitting event the tail agents must fill all the empty positions left by the wing agents promoted to cleaners, and become wing agents themselves.

At k = 1 all the team is organized in an unique formation containing all the n agents1 ; the formation follows U 1 executing its task until a splitting event occurs. In this case, always communicating through the HBSR channel, all the tail agents with odd positions abandon the formation in order to follow the new cleaner composing a new, independent formation as shown in Fig. 7. The same mechanism can be extended to each formation present at the generic time step k. In the case that the tail of the formation is empty the possible vacant wing roles cannot be replaced. If a cleaner concludes its clearing process, then all the agents belonging to the associated formation become helper UAVs. When a UAV becomes a helper it sends a signal of availability on the LBLR channel; all the cleaners then reply on the same channel providing a suitable signal in order to let the helper understand the direction where to move to approach them. Following a suitable selection strategy, the helpers elect a target among the cleaners and move toward it following the direction of origin of the LBLR signal. As soon as the the distance between the helpers and their target gets sufficiently small ( =i.e., HBSR communication is available) a merging event occurs: all the helpers become tails of the selected cleaner. Note that the whole process is completely decentralized.

A summary block-scheme representing the main operative steps of the proposed algorithm is in Fig. 5.

VI. SIMULATIONS

In order to validate the proposed algorithm we have performed a wide range of simulations. Many different kind of road networks have been considered by using real satellite images of urbanized environments, taken by the Google Maps service2 , and imported on the ground of our simulated world in the V-REP simulator 3 . A group of realistically simulated quadrotors has been used as UAV platform. Each UAV has a rigidly attached camera which is therefore subject to the same tilting rotations that the quadrotor performs in order to change its velocity. Image processing has been implemented using the OpenCV4 library, the formation control through Matlab 5 and communication is handled in ROS 6 .

We report the results of a representative simulation which illustrates the main features of the algorithm, see Fig. 8. A video of the simulation can be seen at https://www.youtube.com/watch?v=ESlVS40h1XE. For this simulation we selected a map of an area which is about 500 m 2 wide, located near Lancaster (California). All the roads are approximately 6 m wide and visible from an aerial point of view. All the most frequent possible characteristics of a network of roads are represented: dead-ends, loops and junctions. In this simulation we used a team of 8 UAVs, each of them is equipped with a camera having an field of view of α = 80 • and a resolution L = 256. The desired altitude is pz = 70 m; the clearing algorithm runs at a frequency of f = 8 Hz and the absolute value of the nonzero displacement is |∆| = 0.125 m. The flight control loop that regulates the speed of the quadrotor is based on the controller described in [START_REF] Lee | Semi-autonomous haptic teleoperation control architecture of multiple unmanned aerial vehicles[END_REF] and runs at 100 Hz.

Segmentation of the street has been performed incrementing the contrast of the roads with respect to the background and then using the OpenCV library. As it can be noticed from the video, in spite of a certain amount of noise and outliers in the segmentation the clearing process is executed correctly by the team. The trajectories executed by the cleaner UAVs are represented in Fig. 8 as continuous lines with a different color for each UAV; the colored squares represent a stroboscopic view of the portions of the area visible by the cleaners along their trajectories every 10 s. The dashed lines represent the approximate trajectories executed by each cleaner after it becomes a helper. Notice that in this particular case one of the 8 UAVs never assumes the role of cleaner.

Note that the UAVs do not the street profile precisely; this is because each cleaner tries to push back all its current sub-frontiers in order to clear as much area of the street as possible. Note also how the whole network contained in A is covered by the union of squares representing the fields of view, thus demonstrating that the clearing is successfully achieved.

Figure 9 shows some magnified views of the salient phases of the clearing process summarized in Fig. 8.

In Fig. 10a we report the cumulative percentage of the road cleared by each cleaner at each step. At the first step a high increment is present because in the first frame acquired by the first cleaner all the visible road area contributes to the clearing. During the following steps only the new areas actually matters. Whenever a UAV is promoted to cleaner a new steep increment is present. The cleared area for each robot is obviously constant during the time intervals in which the robot is not a cleaner. The plot shows a remarkably uniform distribution of the workload among all the UAVs used as cleaners in the group.

In Fig. 10b the distances covered by each agent over time are shown: the small peaks along the curves represent the accelerations performed by the helpers going toward a formation to merge with it. Note how velocities are approximately constant, showing that the UAVs do not stop (f) After magenta group received helpers from the yellow group, it splits creating the gold group which then encounters the red group.

VII. CONCLUSION

We have presented an algorithm for aerial clearing of a completely unknown network of roads. Under a limited number of reasonable assumptions the algorithm allows a team of UAVs to complete the task by organizing themselves into formations, without any help from mapping or localization systems. Thanks to these features the robots involved in such a process use a very small amount of memory space. Through a series of simulations we determined the effectiveness of the proposed algorithm.

An interesting prospect for future works is the experimental implementation of the proposed clearing strategy with real aerial robots using, as a preliminary experiment, a reconstructed urban scenario in an indoor environment.

Fig. 1 :

 1 Fig. 1: Main formalism for our scenario.

  The set U k C might be the whole U or just one UAV. The reason for Selection of a motion candidate: at step k the primitive ∆ N is selected because the frontier F k V is located on the north side of the image Īk ; at step k + 1 the motion ∆ NE is instead selected since the frontier F Motion candidate validation: ∆ NE has been selected as possible motion at step k since F k

Fig. 3 :

 3 Fig. 3: Encounter between two cleaners U 1 and U 2 : each cleaner deletes from its frontier the potions that fall into the other's image. U 1 remains without frontiers and the frontier of U 2 is reduced to the only northern part.

  Main notation of the image-based process. ∂M ∂s (b) Image-processing for a real image using the margin box ∂ M.

Fig. 4

 4 Fig. 4

Fig. 6 :

 6 Fig. 6: Red: cleaner; orange: right and left wings; green: tails. All lines represent HBSR communication needed during the clearing and splitting phases; the dashed represent in particular data coming from the targets.

Fig. 7 :

 7 Fig. 7: Splitting strategy.

Fig. 8 :

 8 Fig. 8: Trajectories of the cleaners and stroboscopic view of the covered areas. The dashed orange line represents the boundary of the region A; the white line represents the starting line S. A video of the simulation is available at https://www.youtube.com/watch?v=ESlVS40h1XE

  (a) At start red group is full. (b) First split between red and blue groups. (c) Red and blue groups proceed independently. (d) The helpers created from the blue group merge with the yellow group. (e) Second split of red group generates the magenta one.

Fig. 9 :

 9 Fig. 9: Salient phases of the simulation with quadrotor UAVs performing the clearing task on a real network of roads.

  Distances covered by each UAV.

  Fig. 10

For reasons of clarity we assume n > 3, so that at least a tail is present.

https://maps.google.com/

http://www.coppeliarobotics.com/ Preprint version 5

http://opencv.org/

http://www.mathworks.com/

http://wiki.ros.org/

Preprint version 6

Preprint version

13th ECC