Shiyu Liu
email: shiyu.liu@ls2n.fr

Isabelle Fantoni
email: isabelle.fantoni@ls2n.fr

Abdelhamid Chriette
email: abdelhamid.chriette@ls2n.fr

Decentralized control and state estimation of a flying parallel robot interacting with the environment

Keywords:

Unmanned Aerial Vehicles (UAVs) have great potential to achieve a variety of tasks remotely such as aerial grasping, transporting and manipulating objects. Architectures with multiple UAVs have further enhanced the payload capacity and manipulability of these robots, for instance a Flying Parallel Robot (FPR) where a moving platform is cooperatively supported by multiple quadrotors with passive rigid links. In this paper, we address the vision-based state estimation and decentralized control applied to the multi-UAV parallel robot, taking the FPR as an example. An ArUco marker system is applied to estimate the relative pose of each UAV with respect to the common platform frame, along with the Extended Kalman Filter to reconstruct the robot state without the dependence on any external localization system. The interaction controller is then deployed in a decentralized manner, which is potentially more robust to communication delays or interruptions. The proposed methodology has been validated by real-time experiments demonstrating the teleoperation of the FPR interacting with the environment.

Introduction

Unmanned Aerial Vehicles (UAVs) equipped with additional tools or end-effectors have great potential to accomplish a variety of tasks in remote locations or places considered dangerous to human operators. These aerial robots commonly known as Aerial Manipulators [START_REF] Ollero | Past, Present, and Future of Aerial Robotic Manipulators[END_REF] have been shown capable of performing tasks like grasping, inspection, and transportation, in which the physical interaction with the environment is often involved. Among these aerial manipulators, the ones with multiple UAVs cooperatively supporting a payload or a moving platform have aroused growing interests due to its potential advantages. In fact, the cooperation of multiple UAVs has greatly increased the payload capacity as in cable-suspended aerial robots [START_REF] Bernard | Autonomous Transportation and Deployment with Aerial Robots for Search and Rescue Missions[END_REF][START_REF] Masone | Cooperative transportation of a payload using quadrotors: A reconfigurable cable-driven parallel robot[END_REF][START_REF] Gassner | Dynamic collaboration without communication: Vision-based cable-suspended load transport with two quadrotors[END_REF][START_REF] Loianno | Cooperative transportation using small quadrotors using monocular vision and inertial sensing[END_REF][START_REF] Sanalitro | Fullpose manipulation control of a cable-suspended load with multiple uavs under uncertainties[END_REF] and can achieve the full manipulability in 3dimensional space such as the designs shown in [START_REF] Nguyen | A Novel Robotic Platform for Aerial Manipulation Using Quadrotors as Rotating Thrust Generators[END_REF][START_REF] Park | Pose and posture estimation of aerial skeleton systems for outdoor flying[END_REF]Li, Song, Bégoc, Chriette and Fantoni, 2021b). Similarly to these works, a Flying Parallel Robot (FPR) has been proposed in [START_REF] Six | The Kinematics, Dynamics and Control of a Flying Parallel Robot with Three Quadrotors[END_REF] and validated by experiments in [START_REF] Six | Dynamic modelling and control of flying parallel robots[END_REF], which is composed of a passive parallel architecture supported by a number of quadrotors rigidly connected to the structure by means of spherical joints. By combining the advantages of parallel mechanisms and aerial vehicles, the FPR has shown its potential of accomplishing various tasks by exerting force onto the external objects in the environment [START_REF] Liu | Wrench Estimation and Impedance-based Control applied to a Flying Parallel Robot Interacting with the Environment[END_REF]. A lot of industrial applications can be envisaged such as infrastructure inspection, repairing and pick-and-place operations in remote locations [START_REF] Ruggiero | Aerial manipulation: A literature review[END_REF]. However, the previous works on the FPR and most of other multi-UAV parallel robots still depend on the high-rate exteroceptive measurements of the robot state, which means the necessity of using Motion Capture (MOCAP) systems that are impractical in real-world applications.

The robot localization in multi-UAV parallel robots is an important topic especially for maintaining good inter-UAV positions. Indeed, the access to Global Navigation Satellite Systems (GNSS) in outdoor scenarios has been a solution applied to relatively large-scale systems as in [START_REF] Bernard | Autonomous Transportation and Deployment with Aerial Robots for Search and Rescue Missions[END_REF]; [START_REF] Park | Pose and posture estimation of aerial skeleton systems for outdoor flying[END_REF]. However, the centimetric precision of the GPS is not satisfactory for sufficient calculation of internal configuration avoiding collisions between multiple UAVs, not to mention the degenerated cases where the aerial robots are close to the building or under structures to perform the tasks, or even in GPS-denied environments. Several vision-based solutions have therefore been proposed such as [START_REF] Loianno | Cooperative transportation using small quadrotors using monocular vision and inertial sensing[END_REF]; Li, Ge and Loianno (2021a) using monocular vision and Inertial Measurement Unit (IMU) sensing to localize each individual UAV with respect to the payload. In [START_REF] Loianno | Cooperative transportation using small quadrotors using monocular vision and inertial sensing[END_REF], multiple quadrotors are rigidly connected to a bar-like structure, each estimating its own pose and that of the structure using Visual-Inertial Odometry (VIO) technique. In Li et al. (2021a) for a cablesuspended aerial system, apart from onboard VIO, the information from the embedded monocular camera is further used to estimate the payload and cable states. The vision-based state estimation technique can therefore be applied to other multi-UAV parallel robots as the FPR, where each individual UAV can localize itself relative to the moving platform using intrinsic measurements from the camera and the IMU.

Based on the localization of each UAV with respect to a common body (the payload or the moving platform), a control law can be designed which can potentially be deployed onboard each UAV to reduce the information exchanges and increase the robustness of the whole system. This control methodology is often referred to as a decentralized control scheme extensively applied in swarm robotics and multirobot control community, where each robot is desirably treated as an independent agent running its own perception and control algorithms making the whole system more robust to communication delays or interruptions [START_REF] Coppola | A Survey on Swarming With Micro Air Vehicles: Fundamental Challenges and Constraints[END_REF]. In these systems, the decentralization is often made possible by finding some consensus or maintaining a virtual structure among the agents [START_REF] Chung | A Survey on Aerial Swarm Robotics[END_REF]. The decentralized controller employed in multi-UAV parallel robots is still an open problem, with examples found in the control of a cable-suspended payload by quadrotors using leader-follower [START_REF] Gassner | Dynamic collaboration without communication: Vision-based cable-suspended load transport with two quadrotors[END_REF][START_REF] Gabellieri | A study on force-based collaboration in swarms[END_REF] or force-consensus [START_REF] Thapa | Cooperative Aerial Manipulation with Decentralized Adaptive Force-Consensus Control[END_REF] methods.

In this paper, we propose a state estimation method distributed on each UAV using vision-based technique with ArUco marker system and a decentralized control approach that can also be deployed onboard each UAV with a minimal set of information shared among the UAVs. Our approach allows each individual UAV in the FPR to compute its control based on the onboard measurements and achieve cooperatively an interaction control with sufficiently good estimation of external wrench acting on the robot. This work is therefore taking advantage of the decentralized scheme in control problems of the FPR interacting with the environment based on intrinsic measurements. We have implemented and validated the approach by real-world experiments. It is remarked that in a real-world implementation, a set of onboard sensors is needed so as to estimate certain variables (i.e. attitude and linear velocity) of the quadrotors. In the experiments we conducted, we used the MOCAP system to improve the linear velocity estimates in quadrotor's local frame, which can eventually be replaced by other onboard sensors. The methodology proposed in this paper does not necessarily require any MOCAP system, as the main contribution is to use the vision-based pose estimation technique and decentralized method to control the FPR without dependency on any external localization system.

The rest of the paper is organized as follows. In Section 2, we present the modelling and control of the FPR taking into account the interaction with the environment. Section 3 presents how to estimate the robot state (i.e. pose and velocity) by intrinsic measurements without the need of using any external localization systems. Section 4 addresses the decentralized control scheme for accomplishing the interaction with the environment using only the accessible robot states recovered by the onboard measurements. We present the experimental results of the proposed methodology in Section 5. Finally, we conclude this work and discuss the future works in Section 6.

Dynamic modelling and control

We firstly introduce the modelling and controller for the FPR system, where a moving platform is supported from 𝑛 quadrotors by rigid links (or called legs). The top of each leg is attached to the platform using one Degree-of-Freedom (DoF) revolute joint, while the other end is connected to the quadrotor by means of spherical joint, which is considered attached to the Center of Mass (CoM) of each quadrotor and allows a dynamic decoupling regarding the quadrotor's rotational movements. The kinematic and dynamic models are presented in Section 2.1, based on which external wrench estimation and interaction control methods are formulated in Section 2.2 and Section 2.3.

Kinematic and dynamic models

As shown in Fig. 2, we consider 𝑛 quadrotors that cooperatively support a moving platform (on which an endeffector can be mounted) by rigid legs with passive (i.e. not actuated) joints. The robot configuration can be exclusively represented by the platform's 6D pose and additional DoFs provided by the revolute-joint motions. Therefore, the robot pose is given by a vector where 𝐩 𝑝 ∈ ℝ 3 and 𝔮 𝑝 ∈ ℍ are respectively the 3dimensional position and unit quaternion of the platform with respect to 𝔉 0 , 𝜽 𝑙 = [𝜃 1 , 𝜃 2 , ..., 𝜃 𝑛 ] 𝑇 ∈ ℝ 𝑛 is a vector concatenating the revolute-joint angles of the legs. Notation for relevant variables in this paper is given in Table 1.

𝐪 = [ 𝐩 𝑇 𝑝 𝔮 𝑇 𝑝 𝜽 𝑇 𝑙 ] 𝑇 ∈ ℝ 7+𝑛 (1)
Assuming that the center of the spherical joint is located the CoM of the quadrotor, the geometric relationship between the robot pose and the position of the 𝑖 𝑡ℎ quadrotor can be derived as

𝐩 𝑖 = 𝐩 𝑝 + 𝑟 ⋅ 𝐫 𝑖 + 𝑙 ⋅ 𝐥 𝑖 (2)
where the unit vectors 𝐫 𝑖 and 𝐥 𝑖 are expressed in 𝔉 0 . They can be calculated respectively by 𝐫 𝑖 = 𝐑 𝑝 𝑝 𝐫 𝑖 and 𝐥 𝑖 = 𝐑 𝑝 𝑝 𝐑 𝑙𝑖 𝑙𝑖 𝐥 𝑖 , with 𝑝 𝐫 𝑖 and 𝑙𝑖 𝐥 𝑖 being two constant vectors, 𝑝 𝐑 𝑙𝑖 a rotation matrix between 𝔉 𝑝 and 𝔉 𝑙𝑖 as a function of 𝜃 𝑖 , and 𝐑 𝑝 the rotation matrix of the platform converted from the unit quaternion 𝔮 𝑝 . By concatenating relationships Eq. ( 2) for all the quadrotors, the Inverse Geometric Model (IGM) of the FPR can be expressed by

𝐩 = IGM(𝐪) (3)
with 𝐩 = [𝐩 𝑇 1 , 𝐩 𝑇 2 , ..., 𝐩 𝑇 𝑛 ] 𝑇 ∈ ℝ 3𝑛 a vector concatenating the 3-dimensional positions of all the quadrotors in 𝔉 0 .

Deriving Eq. ( 3), the Inverse Kinematic Model (IKM) can be obtained as

𝐯 = 𝐉(𝐪)𝝂 (4)
where 𝐯 = [𝐯 𝑇 1 , 𝐯 𝑇 2 , ..., 𝐯 𝑇 𝑛 ] 𝑇 ∈ ℝ 3𝑛 is a vector of linear velocities of all the quadrotors in 𝔉 0 , 𝐉(𝐪) ∈ ℝ 3𝑛×(6+𝑛) is a Jacobian matrix, and 𝝂 is the robot velocity given by

𝝂 = [ 𝐯 𝑇 𝑝 𝑝 𝝎 𝑇 𝑝 θ𝑇 𝑙 ] 𝑇 ∈ ℝ 6+𝑛 (5)
with 𝐯 𝑝 , 𝑝 𝝎 𝑝 ∈ ℝ 3 being the linear and angular velocity of the platform expressed respectively in 𝔉 0 and 𝔉 𝑝 , and θ𝑙 ∈ ℝ 𝑛 the derivatives of the leg angles. Knowing the relationship of Eq. ( 3), the robot velocity can be computed from the measurements of quadrotors' linear velocities by the Forward Kinematic Model (FKM) as

𝝂 = 𝐉(𝐪) † 𝐯 (6)
with 𝐉(𝐪) † representing the pseudo-inverse of the Jacobian matrix. Note that to fully actuate the system, the condition of 𝑛 ≥ 3 should be satisfied. When the Jacobian matrix is square (𝑛 = 3), the relation holds true if and only if the Jacobian matrix is not singular (i.e. det(𝐉) ≠ 0).

The dynamic model of the FPR taking into account the external wrench denoted by 𝝉 𝑒 can be described in a matrix form

𝐌(𝐪) ν + 𝐂(𝐪, 𝝂)𝝂 + 𝐠(𝐪) = 𝝉 + 𝝉 𝑒 (7)
where 𝐌(𝐪) ∈ ℝ (6+𝑛)×(6+𝑛) is the generalized inertia matrix, 𝐂(𝐪, 𝝂) ∈ ℝ (6+𝑛)×(6+𝑛) is the Coriolis matrix, 𝐠(𝐪) ∈ ℝ 6+𝑛 is the gravity vector, and 𝝉 is the actuation wrench defined by 𝝉 = [𝐟 𝑇 𝑝 , 𝑝 𝐦 𝑇 𝑝 , 𝐦 𝑇 𝑙 ] 𝑇 ∈ ℝ 6+𝑛 . The actuation wrench can be linked to the thrust forces of the quadrotors by

𝝉 = 𝐉(𝐪) 𝑇 𝐟 (8) with 𝐟 = [𝐟 𝑇 1 , 𝐟 𝑇 2 , ..., 𝐟 𝑇 𝑛 ] 𝑇 ∈ ℝ 3𝑛
a vector concatenating the 3-dimensional thrust forces of all the quadrotors in 𝔉 0 . Note that the dynamic modelling is accomplished using Spatial Vector notation [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], with the Coriolis matrix computed based on a numerical algorithm detailed in [START_REF] Echeandia | Numerical Methods to Compute the Coriolis Matrix and Christoffel Symbols for Rigid-Body Systems[END_REF].

External wrench estimation

Before constructing the interaction control law, a momentum based technique is applied to build an estimate of the external wrench acting on the robot. A momentum-based wrench estimation technique can be adopted according to [START_REF] Liu | Wrench Estimation and Impedance-based Control applied to a Flying Parallel Robot Interacting with the Environment[END_REF], which is summarized as follows.

Let τ𝑒 denote an estimate of the external wrench and  = 𝐌(𝐪)𝝂 be the generalized momentum of the robot. The timedomain expression of the estimated external wrench can be determined based on the computation of the momentum as

τ𝑒 (𝑡) = 𝐊 𝑂 [  (𝑡) -∫ 𝑡 𝑡 0 ( 𝐂(𝐪, 𝝂) 𝑇 𝝂 -𝐠(𝐪) + 𝝉(𝑡) + τ𝑒 (𝑡 -Δ𝑡) ) 𝑑𝑡 - (𝑡 0 ) ] (9)
where

 (𝑡) = 𝐌 ( 𝐪(𝑡) ) 𝝂(𝑡)
is the momentum at time instant 𝑡 with  (𝑡 0 ) initialized by zero, τ𝑒 (𝑡) is the estimate at current timestamp which depends on the previous estimate τ𝑒 (𝑡 -Δ𝑡), and 𝐊 𝑂 is a positive-definite diagonal matrix for the estimation gains. By approximating any two consecutive time steps (supposing Δ𝑡 small enough), the estimation of external wrench by Eq. ( 9) is equivalent to a first-order filtering towards the real values [START_REF] Liu | Wrench Estimation and Impedance-based Control applied to a Flying Parallel Robot Interacting with the Environment[END_REF].

Interaction controller

Based on the estimated external wrench, a reaction behavior of the FPR can be designed during the interaction with the environment. An impedance-based control law with desired wrench tracking ability is implemented in [START_REF] Liu | Wrench Estimation and Impedance-based Control applied to a Flying Parallel Robot Interacting with the Environment[END_REF], in which a generalized form of the desired impedance behavior of the robot can be chosen as (10) with 𝐌 𝑑 , 𝐁 𝑑 and 𝐊 𝑑 ∈ ℝ (6+𝑛)×(6+𝑛) being the positivedefinite diagonal matrices respectively for the desired mass, damping and stiffness coefficients, 𝜺 𝜏 = -τ𝑒 -𝝉 𝑑 𝑒 representing the tracking error of the interaction wrench, and 𝜺 𝑞 (𝐪 𝑑 , 𝐪) the tracking error of the robot pose determined in [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF].

𝐌 𝑑 ( ν𝑑 -ν) + 𝐁 𝑑 (𝝂 𝑑 -𝝂) + 𝐊 𝑑 𝜺 𝑞 (𝐪 𝑑 , 𝐪) = 𝜺 𝜏
Based on the desired impedance in Eq. ( 10), an auxiliary input corresponding to the acceleration of the system can be determined by

𝐮 = ν𝑑 + (𝐌 𝑑 ) -1 ( 𝐁 𝑑 (𝝂 𝑑 -𝝂) + 𝐊 𝑑 𝜺 𝑞 (𝐪 𝑑 , 𝐪) -𝜺 𝜏 ) (11)
which is then applied to the dynamic model in Eq. ( 7) to obtain the thrust force commands as

𝐟 = [ 𝐉(𝐪) 𝑇 ] † ( 𝐌(𝐪)𝐮 + 𝐂(𝐪, 𝝂)𝝂 + 𝐠(𝐪) -τ𝑒 ) (12)
with [ 𝐉(𝐪) 𝑇 ] † being the pseudo-inverse of the transpose of the Jacobian matrix. When 𝑛 = 3, the Jacobian matrix is square and the pseudo-inverse is replaced by a simple matrix inversion, which necessitates det(𝐉) ≠ 0.

We remark that the achievement of thrust force for each quadrotor is handled by the onboard flight controller tracking the attitude and thrust magnitude commands (𝔮 𝑑 𝑖 , 𝑓 𝑑 𝑖 ) computed from each 3-dimensional thrust force vector 𝐟 𝑖 using the method presented in [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF]. It is additionally remarked that this control method requires high-frequency and stable measurements of the robot pose and velocity relative to the global reference frame 𝔉 0 , which is usually known by external localization systems, meaning for indoor applications the dependence on Motion Capture (MOCAP) system.

Vision-based state estimation

In this section, we present the robot state estimation based on the vision-based method using ArUco marker system that can be deployed onboard each quadrotor to avoid using any external localization systems. We hereby demonstrate how to construct the state estimation using real-world data from the onboard visual and inertial measurements.

Relative pose measurement

Supposing that each quadrotor can estimate its pose relative to the platform frame 𝔉 𝑝 , the robot pose can be partially reconstructed, consisting in roll and pitch of the platform and the internal configuration (i.e. leg angles). We therefore implement a real-time pose estimation process using ArUco marker detection algorithm (Romero-Ramirez, Muñoz-Salinas and Medina-Carnicer, 2021). This visionbased estimation method requires the visual marker to be always visible to the camera, an assumption satisfied by carefully mounting the camera on the quadrotor pointing towards the marker.

Consider that an ArUco marker is attached below the platform, which has a marker frame denoted by 𝔉 𝑀 (see Fig. 3(a)). The marker can be detected and tracked by the ArUco detection algorithm using the monocular camera attached onboard each quadrotor, for which the relative pose of the marker expressed in the camera frame (given by 𝔉 𝐶𝑖 ) is estimated and can be represented by a homogeneous matrix written as

𝐶𝑖 𝐓 𝑀 = [ 𝐶𝑖 𝐑 𝑀 𝐶𝑖 𝐩 𝑀 𝟎 𝑇 3 1 ] ∈ ℝ 4×4 (13)
where 𝐶𝑖 𝐑 𝑀 ∈ 𝑆𝑂(3) and 𝐶𝑖 𝐩 𝑀 ∈ ℝ 3 are respectively the rotation matrix and the translation vector for the relative pose of the marker expressed in quadrotor 𝑖's camera frame. Note that 𝟎 3 is a 3-dimensional column vector of zeros. Then, by further multiplying two fixed transformations 𝑝 𝐓 𝑀 and 𝐶𝑖 𝐓 𝑖 calibrated beforehand, which respectively represent the marker pose relative to 𝔉 𝑝 and the camera pose relative to 𝔉 𝑖 , we can obtain an estimate of the platform pose relative to the body-fixed frame 𝔉 𝑖 of quadrotor 𝑖 by

𝑖 𝐓 𝑝 = ( 𝐶𝑖 𝐓 -1 𝑖 ) 𝐶𝑖 𝐓 𝑀 ( 𝑝 𝐓 -1 𝑀 ) (14)
with 𝐶𝑖 𝐓 -1 𝑖 and 𝑝 𝐓 -1 𝑀 being the inverse transformations, knowing that for a homogeneous transformation 𝐴 𝐓 𝐵 between two arbitrary frames 𝔉 𝐴 and 𝔉 𝐵 , its inverse is given by

𝐴 𝐓 -1 𝐵 = 𝐵 𝐓 𝐴 = [ 𝐴 𝐑 𝑇 𝐵 -𝐴 𝐑 𝑇 𝐵 𝐴 𝐩 𝐵 𝟎 𝑇 3 1 ] ∈ ℝ 4×4 (15)
The inverse transformation of Eq. ( 14), 𝑝 𝐓 𝑖 = 𝑖 𝐓 -1 𝑝 , allows each quadrotor to locate itself commonly in the platform frame.

Before reconstructing the robot pose from knowledge of 𝑝 𝐓 𝑖 , the estimated relative pose needs to be further filtered to be able to output at a constant frequency with reduced estimation noises and errors that might occur due to the issues commonly found in vision systems such as motion blur, unfavourable lightening conditions and loss of objects in the camera's angle of view. To achieve this, an Extended Kalman Filter (EKF)-based algorithm is adopted and detailed as follows.

EKF-based pose estimation

The relative pose initially given by 𝑝 𝐓 𝑖 is rearranged by a vector defined as the system state of the EKF

𝐱 = [ 𝑝 𝐩 𝑇 𝑖 𝑝 𝔮 𝑇 𝑖 ] 𝑇 ∈ ℝ 7 (16)
with 𝑝 𝔮 𝑖 ∈ ℍ a unit quaternion converted from the rotation matrix 𝑝 𝐑 𝑖 . To develop the prediction model in the EKF, two reasonable hypotheses are made:

• 𝑝 ṗ𝑖 is supposed to be zero at each instant, since each quadrotor is rigidly attached to the passive structure, and its relative position is constant with non-changing (or slow movements of) leg angles.

• 𝑝 q𝑖 is only affected by the angular velocity of each quadrotor, as the rotational movements are fully decoupled with the passive architecture and the angular velocity of the platform is negligible compared to the fast rotational dynamic response of quadrotors.

Consider now the system equations for the prediction and correction models of the EKF written as [START_REF] Sabatini | Quaternion-based extended kalman filter for determining orientation by inertial and magnetic sensing[END_REF])

𝐱 𝑘 = 𝐟 (𝐱 𝑘-1 , 𝐮 𝑘 ) + 𝐰 𝑘-1 𝐳 𝑘 = 𝐡(𝐱 𝑘 ) + 𝐯 𝑘 (17)
where 𝐮 𝑘 is the system input being the angular velocity 𝑖 𝝎 𝑖 of the quadrotor 𝑖 for each time instant 𝑘, and 𝐰 𝑘-1 , 𝐯 𝑘 are vectors respectively for the prediction and measurement noises modeled by Gaussian noises as 𝐰 𝑘-1 ∼  (𝟎, 𝐐 𝑘-1 ) and 𝐯 𝑘 ∼  (𝟎, 𝐑 𝑘 ) with 𝐐 𝑘-1 , 𝐑 𝑘 the covariance matrices.

Based on the hypotheses, the prediction model and the state transition matrix can be determined by

𝐟(𝐱 𝑘-1 , 𝐮 𝑘 ) = 𝐱 𝑘-1 + ⎡ ⎢ ⎢ ⎣ 𝟎 1 2 𝑝 𝔮 𝑖 • | | | | 0 𝑖 𝝎 𝑖 | | | | ⎤ ⎥ ⎥ ⎦ 𝑇 𝑠 ( 18
)
𝐅 𝑘 = 𝜕𝐟 𝜕𝐱 | | |𝐱 𝑘-1 ,𝐮 𝑘 = [ 𝟏 3×3 𝟎 3×4 𝟎 𝑇 3×4 𝟏 4×4 + 𝛀( 𝑖 𝝎 𝑖 )𝑇 𝑠 ] ( 19
)
where 𝑇 𝑠 is the time step of the prediction process, and 𝛀(.) is a (4 × 4) skew-symmetric matrix associated to a vector of angular velocity 𝝎 commonly found in quaternion-based EKF algorithms [START_REF] Sabatini | Quaternion-based extended kalman filter for determining orientation by inertial and magnetic sensing[END_REF], given by

𝛀(𝝎) = 1 2 [ 0 -𝝎 𝑇 𝝎 -[𝝎] × ] ∈ ℝ 4×4 (20)
with [.] × ∈ ℝ 3×3 being a conventional skew-symmetric matrix of a 3-dimensional vector. Remark that the evolution of quaternion between two iterations in Eq. ( 18) and Eq. ( 19) is written by the first-order and second-order of Taylor series expansion of exp ( 𝛀( 𝑖 𝝎 𝑖 )𝑇 𝑠 ) , as 𝛀( 𝑖 𝝎 𝑖 ) is assumed to be constant within the time interval 𝑇 𝑠 , when 𝑇 𝑠 is small enough.

The correction model and the observation matrix is by the way easy to derive as the full vector of the system state is considered to be directly measured, which gives

𝐡(𝐱 𝑘 ) = 𝐱 𝑘 , 𝐇 𝑘 = 𝜕𝐡 𝜕𝐱 | | |𝐱 𝑘 = 𝟏 7×7 (21)

Robot state reconstruction

After having estimated and filtered the relative pose of each quadrotor expressed in 𝔉 𝑝 , we now focus on how to recover some variables of the robot state that could be sufficient for the control.

The internal configuration (i.e. the leg angle) for each quadrotor was previously computed using the geometric model knowing the global pose of the platform and global positions of quadrotors. It can now be computed with only knowledge of relative position 𝑝 𝐩 𝑖 by

𝜃 𝑖 = arccos ( 𝑝 𝐫 𝐴 𝑖 𝐵 𝑖 ⋅ 𝑝 𝐩 𝐴 𝑖 ‖ 𝑝 𝐫 𝐴 𝑖 𝐵 𝑖 ‖ ⋅ ‖ 𝑝 𝐩 𝐴 𝑖 ‖ ) (22)
with 𝑝 𝐩 𝐴 𝑖 representing the position of the revolute joint's center expressed in 𝔉 𝑝 , and 𝑝 𝐫 𝐴 𝑖 𝐵 𝑖 = 𝑝 𝐩 𝑖 -𝑝 𝐩 𝐴 𝑖 being a vector describing the leg's direction expressed in 𝔉 𝑝 .

For representing the orientation of a rigid body without any global reference frame, we introduce a notion of flat frames similarly as in [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF]. It is a frame (denoted by 𝔉 𝐹 𝑗 attached to an arbitrary body 𝑗) defined with zero roll and pitch relative to a global reference frame 𝔉 0 , remaining an unknown yaw with respect to 𝔉 0 (as shown in Fig. 3(b)).

Therefore, the roll and pitch of a body can be expressed with respect to its flat frame, whose values are the same as the ones represented in any global reference frame 𝔉 0 . We remark that the roll and pitch of each quadrotor can be sufficiently well estimated by the onboard IMU sensor and the filtering technique, which are denoted by 𝐹 𝑖 𝜙 𝑖 and 𝐹 𝑖 𝜗 𝑖 expressed in its flat frame 𝔉 𝐹 𝑖 . Along with the relative orientation between 𝔉 𝑝 and 𝔉 𝑖 estimated and filtered onboard, the platform's orientation is known by the relationship

𝐹 𝑖 𝐑 𝑝 = 𝐹 𝑖 𝐑 𝑖 𝑖 𝐑 𝑝 (23)
where 𝑖 𝐑 𝑝 is the rotation matrix extracted from 𝑖 𝐓 𝑝 and 𝐹 𝑖 𝐑 𝑖 = 𝐑 𝑦 ( 𝐹 𝑖 𝜗 𝑖 )𝐑 𝑥 ( 𝐹 𝑖 𝜙 𝑖 ) is the rotation matrix computed from unit rotations (in 𝑥 and 𝑦 axes) for roll and pitch angles of quadrotor 𝑖.

Finally, the roll and pitch of the platform expressed in its flat frame can be reconstructed by

𝐹 𝑝 𝜙 𝑝 = atan2 ( 𝐹 𝑖 𝐑 𝑝 (3, 2), 𝐹 𝑖 𝐑 𝑝 (3, 3) ) 𝐹 𝑝 𝜗 𝑝 = -asin ( 𝐹 𝑖 𝐑 𝑝 (3, 1) ) (24)
where 𝐑(𝑟, 𝑐) refers to the element at the 𝑟-th row and 𝑐-th column of the matrix.

We notice that the flat frame attached to the platform can be served as a common frame to express all the other states of each quadrotor 𝑖, for which the transformation matrix from 𝔉 𝑖 to 𝔉 𝐹 𝑝 is computed by

𝐹 𝑝 𝐓 𝑖 = 𝐹 𝑝 𝐓 𝑝 𝑝 𝐓 𝑖 (25)
where the translation vector of which are then used to compute the velocity vector of the FPR by the FKM given in Eq. ( 6), in which the linear velocity of the platform is now expressed in 𝔉 𝐹 𝑝 .

Decentralized control

In this section, we focus on how to deploy the estimation and control methods in a decentralized manner allowing each quadrotor to perform the state estimation and its control based on the onboard and intrinsic measurements.

Control problem formulation

From the previous section, we have reconstructed a partial set of robot pose (i.e. the roll and pitch of the platform and the leg angles), with the position and yaw of the platform relative to any global reference frame remained unknown. We remark however that these unknown variables are equivalent to the flat outputs of a single quadrotor [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF], which are commonly controllable by their derivatives (such as linear velocity and yaw rate) for a single quadrotor. We further remark that the first-order derivatives of the unknown variables of the platform can be computed using the kinematic model as presented in Section 3.3. Therefore, we choose the variables on which we have the direct control as a vector

𝝌 = [ 𝐹 𝑝 𝐯 𝑇 𝑝 𝐹 𝑝 ψ𝑝 𝐹 𝑝 𝜗 𝑝 𝐹 𝑝 𝜙 𝑝 𝜽 𝑇 𝑙 ] 𝑇 (27)
where the yaw rate of the platform 𝐹 𝑝 ψ𝑝 can be computed from the angular velocity 𝝎 𝑝 for which the reconstruction is detailed in Section 3.

We remark that the control problem related to driving the controllable variables of Eq. ( 27) to some desired values can be referred to as an easily teleoperable system, since the platform roll and pitch and the leg angles can be regulated by desired values, while a human operator may focus on control of the 3-dimensional position and heading of the platform by sending the velocity commands in 𝔉 𝐹 𝑝 independent of changes in platform angles. To further ensure the smoothness of the evolution of the robot states, the teleoperation commands (denoted by 𝝌 𝑑 ) sent by the operator are further processed by a simple low-pass filter, which outputs the desired derivative commands χ𝑑 to the controller as well. Once received the teleoperation commands composed of 𝝌 𝑑 and χ𝑑 on the controller side, these commands can be parsed to extract the desired trajectory data (i.e. 𝐪 𝑑 , 𝝂 𝑑 and ν𝑑 ). Meanwhile, it is also expected that the operator provides the desired interaction wrench 𝝉 𝑑 𝑒 to the robot.

Communication-based decentralized control

On each quadrotor, the control law itself is based on the interaction controller presented in Section 2.3. However, the wrench estimation and impedance controller require the knowledge of some states of the other quadrotors. One may remark that this information can be known from the desired values and setpoints for these variables, as discussed in Non-Communicating controller in [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF]. However, the sharing of information is obligated to build an unbiased estimation of the external wrench, which is necessary to ensure the stability of the overall control system. Therefore, we introduce for each quadrotor 𝑖 an information vector composed of the minimum set of variables to share with the other quadrotors

𝜻 𝑖 = [ 𝜃 𝑖 𝐹 𝑝 𝐯 𝑇 𝑖 𝐹 𝑝 𝔮 𝑇 𝑖 𝑓 𝑑 𝑖 ] 𝑇 (28)
where 𝜃 𝑖 is the estimated leg angle attached with the quadrotor 𝑖, 𝐹 𝑝 𝐯 𝑖 and 𝐹 𝑝 𝔮 𝑖 are measurements of quadrotor 𝑖's linear velocity and attitude expressed in 𝔉 𝐹 𝑝 , and 𝑓 𝑑 𝑖 is the thrust magnitude command that has been most recently sent to the low-level flight controller.

From 𝑓 𝑑 𝑖 and 𝐹 𝑝 𝔮 𝑖 shared among the quadrotors, the system input 𝐟 computed in the latest control iteration can be recovered and used to estimate the external wrench in the momentum-based observer. It is noted that to reduce the estimation error, we share the actual attitude of the quadrotor instead of the desired values. The actual thrust magnitude produced by each quadrotor might not perfectly corresponds to the desired command, resulting that any mismatch between the actual and desired thrust will be considered as a disturbance and appear in the external wrench estimates. Regarding the potential communication delays when sharing the information among the UAVs, we should make sure that the communication frequency must be higher enough (at least equal to or greater than the control frequency) such that the synchronization of measurements can be negligible.

Once the 3-dimensional thrust force 𝐹 𝑝 𝐟 𝑖 (now expressed in 𝔉 𝐹 𝑝 ) for each quadrotor is computed from the high-level control law in Eq. ( 12), the desired commands for the lowlevel flight controller can be determined [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF], which is given by a vector

𝐮 𝑖 = [ 𝑓 𝑑 𝑖 𝐹 𝑝 𝔮 𝑑 𝑖 𝑇 ] 𝑇 (29)
and achieved by the onboard thrust/attitude controller such as the one presented in [START_REF] Brescianini | Nonlinear Quadrocopter Attitude Control[END_REF].

Adjustment on attitude commands

As all the quadrotors do not share a common reference frame now, the attitude command for each quadrotor 𝐹 𝑝 𝔮 𝑑 𝑖 expressed in 𝔉 𝐹 𝑝 should be further adjusted to be coherent with its onboard attitude. We consider that the onboard attitude of each quadrotor is expressed in an unknown inertial frame denoted by 𝔉 0𝑖 , which shares common roll and pitch angles with those expressed in 𝔉 𝐹 𝑝 (angles along axes perpendicular to the gravity vector) but a different yaw. Therefore, the correction of attitude commands is simply to correct the yaw setpoints which should be coherent with those of onboard measurements.

We can know from 𝐹 𝑝 𝐑 𝑖 a relative yaw between quadrotor 𝑖's frame and 𝔉 𝐹 𝑝 , denoted by 𝐹 𝑝 𝜓 𝑖 , and a measurement of quadrotor 𝑖's yaw denoted by 0𝑖 𝜓 𝑖 expressed in its unknown reference frame. The corrected yaw setpoint expressed in 𝔉 0𝑖 can therefore be computed by

0𝑖 𝜓 𝑑 𝑖 = 0𝑖 𝜓 𝑖 + ( 𝐹 𝑝 𝜓 𝑑 𝑖 -𝐹 𝑝 𝜓 𝑖 ) (30)
with 𝐹 𝑝 𝜓 𝑑 𝑖 the desired yaw originally obtained from the desired attitude 𝐹 𝑝 𝔮 𝑑 𝑖 .

Discussion on the overall control

The overall control system is distributed onboard each quadrotor, with the decentralized control scheme depicted in Fig. 4. The stability of the overall control system to the external disturbances can be ensured as long as the estimation errors remain bounded [START_REF] Liu | Wrench Estimation and Impedance-based Control applied to a Flying Parallel Robot Interacting with the Environment[END_REF]. This assumption can be considered true as all the UAVs share their most recent attitudes and thrust magnitude commands, making the estimation error performed by each individual quadrotor being bounded. The estimation errors due to measurement noises or communication delays on the robot states might further affect the stability of the system, which is however not studied in this paper.

Experimental results

In this section, we present the implementation and experimental results of the proposed methods applied to the FPR for interacting with the environment. A video of experiments can be found at https://youtu.be/eTaze1PvIP0.

Experimental setup

The prototype of the FPR studied within this work is composed of a triangle-form platform and three legs attached with quadrotors, as shown in Fig. 1. The aluminium platform weighs 250 g and the rigid legs are each 1.043 m long with a mass of 66 g. The custom quadrotors used in the FPR are based on a 34-cm Lynxmotion Crazy2Fly frame, controlled by an open-source flight autopilot which is based on a custom-built PX4 firmware (software) and Pixhawk 5 (hardware) [START_REF] Meier | PX4: A node-based multithreaded open source robotics framework for deeply embedded platforms[END_REF]. The onboard flight controller has been experimentally tuned with aggressive attitude and angular rate gains to ensure rapid convergence to the attitude setpoints. On each quadrotor, a Raspberry Pi 4B is mounted as a companion computer to handle the communication with the Pixhawk 5 via Ethernet by UDP protocol. The communication between quadrotors is handled by ROS2 Galactic, through a 5-GHz Wifi network.

The ArUco marker with ID 0 and the size of 0.109 m is selected from the dictionary "ARUCO_MIP_36h12". Images are captured by the onboard Raspberry Pi V2 cameras at 20 Hz with the relative pose estimated by ArUco detection algorithm and filtered by EKF at the same rate. The interaction controller with external wrench estimator is run onboard each quadrotor at 50 Hz, with information shared among the quadrotors at the same rate. The onboard flight controller is run at 250 Hz.

During the experiments, the FPR was flown by teleoperation using a 4-axis gamepad controller inside a 4×6×3.5 m flight arena equipped with an 8-camera Qualisys MOCAP system. The MOCAP provides the ground truth for postflight analysis. We remark that MOCAP data is used for several technical reasons that can be totally avoided in future works

• For safety reason: emulating the relative pose 𝑝 𝐓 𝑖 with representative noises in [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF] if the detection is lost for more than 0.5 s (corresponding to 10 control iterations) or a false detection is identified (with translation error >10 cm or rotation error >10°).

It may be handled by redundant detection sources in the future.

• For practical reason: emulating high-level perception of the environment (replacing human operator) to obtain the platform velocity commands in 𝑥, 𝑦 and yaw axes in Experiment II. It can be achieved by adding a perception system with an additional camera on the platform.

• Sensor fusion onboard the quadrotors to improve the linear velocity estimates of each quadrotor with realistic noises added. It could thus be replaced by optical flow (estimating 𝑥, 𝑦-axis velocities) and distance (for correcting 𝑧-axis velocity) sensors. An additional analysis of the impact of velocity estimation noises and errors is given in the appendix.

It should be noted that MOCAP data is used in a way that only the local frame is needed (such as for improving the body-frame linear velocity estimates, and for emulating camera data), which can be replaced by other onboard sensors and cameras, making the proposed methodology still valid for controlling the robot without knowing any global frame associated to an external localization system.

Experiment I: precise positioning and pick-up of load

In the first experiment, the platform is teleoperated by a human operator for the precise positioning towards a target (tennis ball) hanging in the air on which a load is suspended corresponding to total weight of about 400 g, as shown in Fig. 5(a). The robot configuration is maintained with a flat platform orientation (zero roll and pitch) and leg angles of 50°, a configuration beneficial to ensure the visibility towards the marker from the onboard cameras of all the quadrotors. The operator thus attempts to pick up the load within 45 s after the take-off.

We first evaluate the real-time pose estimation by ArUco detection algorithm and EKF-based filtering. The results of the estimated, filtered values and the ground truth known by MOCAP are recorded and plotted in Fig. 6, with the rootmean-square error (RMSE) summarized in Table 2. Note that each axis of the relative position and orientation are separately analyzed. In Fig. 6, it can be noticed that the outliers of the estimated values emphasized within the red dots have been well filtered from the EKF outputs. The RMSE of the filtered values relative to the ground truth are all bounded within 3 cm for relative position and 3°f or relative orientation, which has been improved compared to the initial estimation results from the ArUco detection algorithm.

The precise positioning of the platform and the pickup of load are then evaluated. The operator was asked to perform the same positioning task 6 times to reduce the human factor affecting the results. However, from the pilot's point of view, the system is quite easy to manipulate, similar to the manual pilot of a single quadrotor without much difficulties. The results on the evolution of positioning errors for the 6 tasks (5 success and 1 failure) are shown in Fig. 7, with successful tasks well identified by the positioning error being stabilized below 15 cm. The 𝑧-axis external force of the platform estimated by one quadrotor during the successful task 1 ○ is plotted in Fig. 8, from which the pick-up of the load is illustrated with the estimated values converged to the ground truth of around -4 N. The overall success rate (5/6) of the pick-up tasks by teleoperation seems very encouraging.

In terms of attitude tracking of quadrotors, the results shown in Table 3 and Fig. 9 demonstrate that the lowlevel attitude control has been well performed. Note that the attitude tracking results are represented by ZYX (yawpitch-roll) Euler angles. The yaw angles of quadrotors are estimated from the outputs of the onboard pose estimation process. The magnitude of yaw tracking errors is around 10°d ue to the camera calibration and estimation errors, which is within acceptable range and doesn't affect the performance of the system due to the characteristic of quadrotors that yaw motion is decoupled with the thrust force the quadrotor can generate.

Experiment II: contact-based interaction with a board

In the second experiment, the FPR is controlled to interact with a 30 × 30 cm wooden board in the environment and exert desired contact force on the board. A region of dimension 10 × 10 cm is pointed out on the surface of the board, with its centre considered as the target for contactbased interactions. Fig. 5(b) illustrates the interaction experiment performed by the FPR. The teleoperation still involves the human in the loop, which is undesired for validating the methodology independent of the proficiency of human operators. Therefore, an emulated pose estimation between the target and the platform is set up in this experiment running at the highest level to generate the desired velocity on 𝑥 and 𝑦 axes 𝐹 𝑝 𝐯 𝑝,𝑥∕𝑦 and yaw rate 𝐹 𝑝 ψ𝑝 using the MOCAP data, leaving the 𝑧-axis linear velocity of the platform, the desired contact force and robot configuration (i.e. platform orientation and leg angles) controlled by the human operator. This however may be done by adding an additional camera on the platform perceiving the contact surface and estimating the relative pose using model-based or learningbased algorithms.

Fig. 10 shows the evolution of the positioning error between the platform and the target, and Fig. 11 are the desired and estimated values of the contact force between the platform and the board. It can be seen that during the 7: Evolution of 3-dimensional positioning of the platform with respect to the target in 6 performed tasks in Experiment I. Note that the black dashed line is the criterion for a successful task (≈ 15 cm). The blue and purple curves (Success 1 ○ and 3 ○) are two successful cases where the load was released before the end of 45 s.

steady contact (the time interval of [32, 55] s) the positioning error is stabilized at a value around 5 cm. From Fig. 11, the estimated contact force well tracks the desired values of 5 N during the interaction, except the peak between [27, 32] s probably because an impact occurs at the beginning of the UAV These results have therefore shown that the proposed decentralization of the state estimation and interaction control is capable of positioning the platform with good precision and accomplishing interactions with the environment by the platform. In addition, we have also observed in the experiment that the decentralized controller can better deal with the communication delays or interruptions compared to the centralized one especially in the emergency case where the communications among the quadrotors were lost. While a crash is expected and often experienced with a centralized controller, each quadrotor in a decentralized scheme is still capable of making use of desired states as investigated in [START_REF] Liu | Decentralized Control and Teleoperation of a Multi-UAV Parallel Robot Based on Intrinsic Measurements[END_REF] to maintain its own control and land properly without crashing.

Analysis of failure cases

In real-world applications of our proposed method, there exists failure cases where the vision-based detection and pose estimation technique may be defeated. This is principally due to the fact that the ArUco marker detection system necessitates direct line-of-sight to the marker. In some specific leg configurations, such as when angle between the platform's plane and the leg's direction is large, the marker visibility from the camera attached on the quadrotor is worsen, which results in an increase of the marker detection losses by the ArUco detection algorithm. Fig. 12 demonstrates examples of captured images by the onboard camera of one quadrotor during the experiment, in which Fig. 12a is a successful detection case in normal robot configuration and Fig. 12b shows a failure case with a wide leg configuration.

The marker visibility issue may potentially be resolved by carefully regulating the leg angles within a safe interval (such as between 35°and 65°), which was validated during the experiments conducted in this work. This is a reasonable assumption since the potential application of the FPR consists in exerting forces onto the external environment with a configuration where the platform is above the quadrotors. An opposite configuration (platform underneath quadrotors) has less interests and the transitions between two sets of configuration require the crossing of singularity, which is totally out of the scope of this paper. However, the failure cases related to environmental conditions such as brightness changes still present the frustrating aspect of the visionbased method. In the experiments of this work, we must have favorable lightening conditions from the side way of the robot to ensure the vision-based detection algorithm to be working correctly.

Conclusion

In this paper, we present the distributed state estimation technique based on intrinsic measurements from the monocular camera and the IMU to recover the robot states without any external localization system. An impedancebased interaction control method enhanced with the external wrench estimation and desired wrench tracking capability has been deployed in a decentralized manner, allowing the whole system to interact with the environment. The experiments conducted on a flying parallel robot have shown the capability of this robot to accomplish real manipulation tasks, such as contact-based inspection or repair on the surface of industrial infrastructures.

Future works can continue on the implementation of the fully MOCAP-independent FPR equipped with a complete set of sensors to provide enough onboard information. A redundant detection algorithm might also be investigated, using multiple sources of information to recover the robot pose such as the IMU attached to the platform, and/or distancebased measurement units. The perception of environment may be done as well providing high-level commands to alleviate the difficulties related to the teleoperation by human operators.

Figure 1 :

 1 Figure 1: A flying parallel robot composed of a moving platform supported by three quadrotors with rigid legs.

Figure 2 :

 2 Figure 2: System description of a generic flying parallel robot. The frames are defined by axes in red (𝑥), green (𝑦), blue (𝑧).

Figure 3 :

 3 Figure 3: Definition of frames for the pose estimation on each quadrotor. (a) frames for the detected relative pose and the fixed transformations; (b) flat frames attached to the platform and each quadrotor 𝑖. Note that the transparent bodies represent the poses in their respective flat frames.

Figure 4 :

 4 Figure 4: Decentralized control diagram. The blocks within the red dashed rectangle are implemented on each UAV. The interaction control block includes the external wrench estimation and impedance-based controller.

Figure 5 :

 5 Figure 5: Experimenting scenarios of proposed state estimation and decentralized control methods. (a) precise positioning and pick-up of load; (b) contact-based interaction with a board.

Figure 6 :

 6 Figure 6: Relative pose of UAV 1 in the platform frame over the course of the first successful test in Experiment I, with ArUco estimates (blue curves), filtered outputs (green curves) and ground truth (red dashed curves). (a)-(c): 𝑥, 𝑦, 𝑧 elements of the relative position 𝑝 𝐩 1 , (d)-(g): 𝑤, 𝑥, 𝑦, 𝑧 elements of the unit quaternion 𝑝 𝔮 1 representing the relative orientation.

Figure 8 :

 8 Figure 8: Evolution of the 𝑧-axis external force estimates of the platform during the successful experiment (Success 1 ○). Note that the curves of the force estimates of all the quadrotors are overlapped.

Figure 9 :

 9 Figure 9: Attitude tracking results of quadrotor 3 during one instance of successful pick-up experiments.

Figure 10 :

 10 Figure 10: Evolution of 3-dimensional positioning error during the contact-based interaction. (a) positioning error over the whole flight; (b) positioning error during the steady contact.

Figure 11 :

 11 Figure 11: Tracking of desired contact force during the interaction experiment.

Figure 12 :

 12 Figure 12: Images captured by the onboard camera of one quadrotor. (a) successful detection instance; (b) failure instance.

  𝑆𝑂(3) rotation matrix of the platform relative to 𝔉 0 𝔉 𝑝 platform frame, with CoM located at point 𝑃 𝐑 𝑖 ∈ 𝑆𝑂(3) rotation matrix of quadrotor 𝑖 relative to 𝔉 0 𝔉 𝑙𝑖 , 𝔉 𝑖 leg 𝑖's frame and quadrotor 𝑖's frame, with 𝐯 𝑝 , 𝐯 𝑖 ∈ ℝ 3 linear velocity of the platform, quadrotor 𝑖 in 𝔉 0 CoM located respectively at point 𝐴 𝑖 and 𝐵 𝑖 𝑝 𝝎 𝑝 ∈ ℝ 3

	Table 1
	Notation and symbols
	𝔉 0	global reference frame	𝐑 𝑝 ∈ body-frame angular velocity of the platform
	𝔉 𝐹 *	flat frames attached to an arbitrary body	𝑖 𝝎 𝑖 ∈ ℝ 3	body-frame angular velocity of quadrotor 𝑖
	𝐩 𝑝 , 𝐩 𝑖 ∈ ℝ 3 position of the platform, quadrotor 𝑖 in 𝔉 0	𝐟 𝑝 ∈ ℝ 3	force exerted on the platform expressed in 𝔉 0
	𝔮 𝑝 , 𝔮 𝑖 ∈ ℍ	unit quaternion of the platform, quadrotor 𝑖	𝑝 𝐦 𝑝 ∈ ℝ 3	moment exerted on the platform expressed in 𝔉 𝑝
		relative to 𝔉 0	𝐦 𝑙 ∈ ℝ 𝑛	a vector concatenating moments around the legs
	𝜽 𝑙 ∈ ℝ 𝑛 θ𝑙 ∈ ℝ 𝑛	a vector concatenating the leg angles 𝜃 𝑖 a vector concatenating the rates of leg angles	𝐟 𝑖 ∈ ℝ 3 𝑓 𝑖 ∈ ℝ	thrust force of quadrotor 𝑖 expressed in 𝔉 0 thrust force magnitude of quadrotor 𝑖
	𝐫 𝑖 , 𝐥 𝑖 ∈ ℝ 3	unit vectors of directions ⃖⃖⃖⃖⃖⃖ ⃗ 𝑃 𝐴 𝑖 , ⃖⃖⃖⃖⃖⃖⃖ ⃗ 𝐴 𝑖 𝐵 𝑖	𝟎 𝑘	a 𝑘-dimensional vector of zeros
	𝑟, 𝑙 ∈ ℝ	platform radius (length of ⃖⃖⃖⃖⃖⃖ ⃗ 𝑃 𝐴 𝑖 ) and leg length 𝟏 𝑘×𝑘	an identity matrix of dimension (𝑘 × 𝑘)

  𝐹 𝑝 𝐓 𝑝 is zero and the rotation matrix of 𝐹 𝑝 𝐓 𝑝 is similarly computed by 𝐹 𝑝 𝐑 𝑝 = 𝐑 𝑦 ( 𝐹 𝑝 𝜗 𝑝 )𝐑 𝑥 ( 𝐹 𝑝 𝜙 𝑝 ). From the rotation matrix 𝐹 𝑝 𝐑 𝑖 , the quadrotor 𝑖's attitude expressed in 𝔉 𝐹 𝑝 represented by the unit quaternion 𝐹 𝑝 𝔮 𝑖 can be also computed. This further allows to express the linear velocities of all the quadrotors commonly in 𝔉 𝐹 𝑝 by

𝐹 𝑝 𝐯 𝑖 = 𝐹 𝑝 𝐑 𝑖 𝑖 𝐯 𝑖 (26)

Table 2

 2 Root-mean-square error of the ArUco detection and EKF-based filtering results of the relative pose between each quadrotor and the platform. Results are calculated based on 6 repetitive tasks performed in Experiment I.

Table 3

 3 Root-mean-square error on the attitude tracking of the quadrotors during one successful test in Experiment I.

Table 4

 4

	, which demonstrates a good maintenance of robot
	configuration during the interaction experiments. Although
	the detection failures and interruptions have been notified
	in the experiments, the percentage of successful detection
	presented in Table 4 (from 88.3% to 95.3% for all the
	quadrotors) has been an encouraging progress applying the
	vision-based state estimation technique in the real-world
	scenario.

Table 4

 4 Root-mean-square error on the tracking of the platform orientation and leg angles (left) and successful detection rate (right) in two experiments. Note that results of experiment 1 are computed by the average values of 6 performed tasks.

Acknowledgment

The authors acknowledge Julian Erskine and Damien Six for their valuable suggestions and support on this work. This work was carried out within the framework of ROBOTEX 2.0, the French Infrastructure in Robotics under the grants ROBOTEX (EQUIPEX ANR-10-EQPX-44-01) and TIRREX (EQUIPEX+ grant ANR-21-ESRE-0015).

Conflict of interest

The authors declared no potential conflict of interest with respect to the research, authorship and publication of this work.

Appendix Analysis of velocity estimation errors

The experimental validation of the proposed method in the paper was still dependent on the Motion Capture (MO-CAP) system. One of the major reason was to improve the linear velocity estimates locally measured by the quadrotors IMU sensors which inevitably present the drifts in velocity measurements. Although it is more a technical concern to replace the MOCAP by onboard sensors such as optical flow and distance sensor to correct the body-frame velocity estimates by sensor fusion techniques, it is still important to prove that our method is valid with acceptable range of estimation errors and noises of the quadrotor linear velocities.

To achieve this, we conducted ROS-based real-time simulations under Gazebo environment (as shown in Fig. 13), with Gaussian noises added both on the pose estimation and velocity measurements of quadrotors. We considered the relative pose estimation errors as 4 cm in position and 4°in orientation, while the linear velocity estimates are disturbed by Gaussian noises with the standard derivation as 0.05m/s on each axis.

In the simulation, we controlled the FPR by sending the velocity commands to the decentralized controller, equivalently to the real-world experiments we analyzed in the paper. Fig. 14 shows the disturbed linear velocities of three quadrotors in the simulation. The platform velocity and its desired values are presented in Fig. 15, and the tracking of the platform orientation is demonstrated in Fig. 16, which together show that the robot is still well controlled by teleoperation even with the disturbed measurements. In reality, the measurements noises and drifts in the quadrotor linear velocities can be improved by equipping more complete set of onboard sensors, which allow the controller to maintain its performance that has been proven in this paper.


Baptis Delaporte

Thomas Arsou

Gunnar Kepp

Swen Jullien

Christophe Menkes

Thomas Ibanez

StormR: An R package to quantify and map the tropical storms and cyclones' winds characteristics

. This database provides a fairly comprehensive record of tropical storms and cyclones with a 3-hours temporal resolution since 1841. However any storm track data can be used as long as the mandatory fields are provided.

Storm track data can be extracted using a specified point location, a user defined spatial polygon shapefile, a country or a cyclone basin name. The main functions of the StormR package allow us to generate wind speed and direction fields as reconstructed from storm track data and a parametric cyclone model. Different models and models combinations can be chosen by the user. By default the spatial resolution is set to 2.5 min (~4.5 km at the equator), but a finer spatial resolution of 30 s (~1 km at the equator) and coarser spatial resolutions of 5 min (~9 km at the equator) or 10 min (~18.6 km at the equator) can be set. The temporal resolution is set to 1 hour by default but finer resolutions of 45 min, 30 min, or 15 min can be set.

Once wind speed is generated for each cell or specific location and each time step, StormR functions can compute summary statistics on wind speed over the lifespan of a storm. Summary statistics encompass the maximum sustained wind speed, the power dissipation index or total power dissipated by a tropical storm [START_REF] Emanuel | The power of a hurricane: An example of reckless driving on the information superhighway[END_REF][START_REF] Emanuel | Increasing destructiveness of tropical cyclones over the past 30 years[END_REF], and the duration of exposure to winds reaching defined speed thresholds. By default the duration of exposure is computed for each Saffir-Simpson Hurricane Scale threshold values for tropical cyclone categories, i.e., 33,43,50,58, and 70 𝑚.𝑠 -1 [START_REF] Simpson | The Hurricane Disaster-Potential Scale[END_REF], but can be defined by the user.

Statement of need

Globally, an average of 86 tropical cyclones per year occurred over the past four decades [START_REF] Murakami | Detected climatic change in global distribution of tropical cyclones[END_REF]. These disturbances cause severe damages to ecosystems, infrastructures, and peoples [START_REF] Costanza | The Value of Coastal Wetlands for Hurricane Protection[END_REF][START_REF] Parks | Tropical cyclone exposure is associated with increased hospitalization rates in older adults[END_REF][START_REF] Uriarte | Hurricane María tripled stem breaks and doubled tree mortality relative to other major storms[END_REF]. Climate change has likely increased the proportion of category 3-5 tropical cyclones on the Saffir-Simpson hurricane wind scale over the past four decades and climate scientists are predicting with high confidence that the proportion of the most intense and potentially the most destructive tropical cyclones (category 4-5) would increase by +10% even if warming is limited to 1.5°C (Intergovernmental Panel on Climate Change, 2023). StormR R package responds to the need for an easy to use tool that helps to better understand winds generated by storms and tropical cyclones in relation to potential damages.

To our knowledge two R packages are available on R CRAN: HurreconR [START_REF] Boose | HurreconR: Models hurricane wind speed, wind direction, and wind damage[END_REF] and hurricaneexposure [START_REF] Anderson | hurricaneexposure: Explore and map county-level hurricane exposure in the United States[END_REF]. As highlighted by their names, these two packages have a strong focus on the North American basin where tropical cyclones are named hurricanes. Another important limitation of those packages is that they both rely on a single model to reconstruct wind speed. The HurreconR package relies on the HURRECON model [START_REF] Boose | Landscape and Regional Impacts of Hurricanes in New England[END_REF][START_REF] Boose | Landscape and Regional Impacts of Hurricanes in Puerto Rico[END_REF], a modification of the [START_REF] Holland | An Analytic Model of the Wind and Pressure Profiles in Hurricanes[END_REF] model, and the hurricaneexposure package relies on the Willoughby's model [START_REF] Willoughby | Parametric Representation of the Primary Hurricane Vortex. Part II: A New Family of Sectionally Continuous Profiles[END_REF]. Many model exists and none of these is the best for all tropical cyclones and storms [START_REF] Yan | Research progress on tropical cyclone parametric wind field models and their application[END_REF]. Compared to those packages, the StormR R package allows to reconstruct wind behaviour for tropical storms and cyclones anywhere, anytime, and this with a set of models including those used in HurreconR and hurricaneexposure packages.

Package workflow

The workflow of the stormR R package consists of four main consecutive steps:

1-The defStormsDataset() function is used to create a stormsDataset object from a NetCDF file, e.g., exported from a storm track database such as IBTrACS database (https://www.ncei.noaa.gov/products/international-best-track-archive) [START_REF] Knapp | The International Best Track Archive for Climate Stewardship (IBTrACS): Unifying Tropical Cyclone Data[END_REF][START_REF] Knapp | International Best Track Archive for Climate Stewardship (IBTrACS) Project[END_REF].

2-

The defStormsList() function is used to perform queries on the stormsDataset object and extract storms of interest based on their location and date of occurrence. The track data for these storms are then stored into a stormsList object. The plotStorms() function can be used to see the extracted data.

3-

The temporalBehaviour() and spatialBehaviour() functions are used to compute wind speed and direction and three summary statistics: the maximum sustained wind speed, the power dissipation index, and the duration of exposure to winds reaching a given wind speed along the life span of the cyclones. The temporalBehaviour() function returns time series at specific point locations (list of data frames) while the spatialBehaviour() function returns spatial data (rasters).

4-

The writeRast() and plotBehaviour() functions are used to export or see the final results.

Figure 1 :

 1 Figure 1: Simplified workflow and main functions of the stormR R package. External storm track database in grey, main functions in blue, and R objects created by the stormR functions in yellow.

Figure 2 :

 2 Figure 2: Example of outputs for the tropical cyclone Pam (2015) nearby Vanuatu with default settings. Storm track data (SSHS=Saffir-Simpson Hurricane Scale, TD=tropical depression, TS=tropical storm), maximum sustained wind (MSW), power dissipation index (PDI), and duration of exposure to wind stronger than 58 𝑚.𝑠 -1 . Plots were generated with the plotStorms() and plotBehaviour() functions as shown in the spatialBehaviour article.

Delaporte et al. (2024). StormR: An R package to quantify and map the tropical storms and cyclones' winds characteristics. Journal of Open Source Software, 9(93), 5766. https://doi.org/10.21105/joss.05766.

Acknowledgements

Work was supported by a Hermon Slade Foundation grant grant HSF 19105. GK was partially supported by the Montpellier Advanced Knowledge Institute on Transitions (MAK'IT) as part of a Visiting Scientist Fellowship while working on parts of the manuscript.