Gianluca Antonelli

Filippo Arrichiello

Stefano Chiaverini

Paolo Robuffo

Adaptive trajectory tracking for quadrotor MAVs in presence of parameter uncertainties and external disturbances

The paper presents an adaptive trajectory tracking control strategy for quadrotor Micro Aerial Vehicles. The proposed approach, while keeping the typical assumption of an orientation dynamics faster than the translational one, removes that of absence of external disturbances and of perfect symmetry of the vehicle. In particular, the trajectory tracking control law is made adaptive with respect to the presence of external forces and moments, and to the uncertainty of dynamic parameters as the position of the center of mass of the vehicle. A stability analysis as well as numerical simulations are provided to support the control design.

I. INTRODUCTION

Over the last years, the use of Micro Aerial Vehicles (MAVs) as robotic platforms has witnessed the growing attention of the robotics community. Typical envisioned tasks for MAVs involve traffic surveillance, monitoring of air pollution, area mapping, inspection of collapsed building or dangerous sites, and agricultural applications such as pesticide spraying. Among the different advantages offered by MAVs, their capabilities to take off and land in limited spaces and to hover above targets are of crucial importance.

Several research groups have recently focused their activities on the development of high-performance flight control and motion planning algorithms for single and/or multiple MAVs; for the experimental validation, they usually exploit suitable indoor facilities equipped with an external visual tracking system for fast and accurate state estimation [START_REF] Michael | The GRASP multiple micro-UAV testbed[END_REF], [START_REF] Franchi | Bilateral shared control of multiple quadrotors[END_REF]. At the same time, other research groups have addressed the problem of estimating the MAV position/orientation online by only relying on onboard sensing (usually, cameras and IMU) and (limited) computation capabilities, see [START_REF] Weiss | Monocular-SLAM-based navigation for autonomous micro helicopters in GPS-denied environments[END_REF], [START_REF] Shen | Autonomous indoor 3D exploration with a micro-aerial vehicle[END_REF], [START_REF] Bry | State estimation for aggressive flight in GPS-denied environments using onboard sensing[END_REF], [START_REF] Grabe | On-board velocity estimation and closed-loop control of a quadrotor UAV based on optical flow[END_REF]. Furthermore, the potentiality of MAVs has been demonstrated in cooperative transportation tasks [START_REF] Fink | Planning and control for cooperative manipulation and transportation with aerial robots[END_REF], aerial grasping of moving objects [START_REF] Spica | Aerial grasping of a moving target with a quadrotor UAV[END_REF], and in complex missions involving distributed sensing and coordination among several MAVs [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with time-varying topology[END_REF], [START_REF] Franchi | Shared Control: Balancing Autonomy and Human Assistance with a Group of Quadrotor UAVs[END_REF], [START_REF] Giordano | A Passivity-Based Decentralized Strategy for Generalized Connectivity Maintenance[END_REF].

For all the applications involving MAVs, ensuring a good and robust performance of the underlying flight controller represents a fundamental requirement. The works in [START_REF] Castillo | Stabilization of a mini rotorcraft with four rotors[END_REF], [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to VTOL drones[END_REF] offer an interesting introduction to the theory involved, as well as some experimental results concerning quadrotors vehicles, i.e., MAVs equipped with four thrusters aligned in the same direction -a nowadays widespread solution in research labs. Quadrotor MAVs are, as well-known, underactuated mechanical systems (6 Degrees-Of-Freedom (DOFs) vs. 4 control inputs). A common strategy to handle the quadrotor motion is to only control its 3D position and yaw angle, since these are flat outputs for the system [START_REF] Fliess | Flatness and defect of nonlinear systems: Introductory theory and examples[END_REF]. To counteract external disturbances such as, e.g., the wind, integral-based actions may be implemented; however, without a proper knowledge of the dynamics under investigation, the use of adaptive or integral control actions may eventually feed the closed loop system with an additional disturbance rather than reducing it, see, e.g., [START_REF] Antonelli | On the use of adaptive/integral actions for 6-degrees-offreedom control of autonomous underwater vehicles[END_REF]. At the best of our knowledge, the only adaptive control for MAVs is proposed in [START_REF] Palunko | Agile load transportation. Safe and efficient load manipulation with aerial robots[END_REF]; in that solution, the mathematical model is rewritten in order to make clear its linear dependency with respect to the sole center of gravity position, which is then successively used in a feedback linearization approach.

This paper presents an adaptive control technique for trajectory tracking of quadrotors that removes the assumption of known dynamic parameters (in particular the position of the center of mass) and considers the presence of constant external forces/moments. With respect to [START_REF] Palunko | Agile load transportation. Safe and efficient load manipulation with aerial robots[END_REF], this paper proposes a different adaptive approach, it explicitly considers the external disturbances in the stability analysis, it isolates a minimal set of parameters to be identified, and it finally takes into account the effect of a wrong projection from desired forces to thrusters velocities due to a poor barycenter estimation. A detailed mathematical analysis and numerical simulations show the efficacy of the proposed approach. The proposed adaptive controller would consistently improve the reliability/robustness of MAVs operations, especially in outdoor environments (wind) or in applications involving picking and releasing of small loads (barycenter shift).

II. MODELING

A. Kinematics

A rigid body is completely described by its position and orientation with respect to a reference frame Σ i , Oxyz assumed earth-fixed and inertial. Let η 1 ∈ R 3 η 1 = x y z T be the vector of the body position coordinates in an earthfixed reference frame. The vector η1 is the corresponding time derivative (expressed in the earth-fixed frame). If one defines ν 1 = u v w T as the linear velocity of the origin of the body-fixed frame Σ b , O bx b y b z b with respect to the origin of the earth-fixed frame, expressed in the bodyfixed frame (from now on: body-fixed linear velocity), then the following relation between the defined linear velocities holds:

ν 1 = R B I η1 , (1)
where R B I is the rotation matrix expressing the transformation from the inertial frame to the body-fixed frame.

Define η 2 ∈ R 3 as

η 2 = φ θ ψ T ,
the set of body Euler-angle coordinates in a earth-fixed reference frame. Those are commonly named roll, pitch and yaw angles and correspond to the succesive elementary rotations around x, y and z in the fixed frame [START_REF] Siciliano | Robotics: modelling, planning and control[END_REF]. The vector η2 is the corresponding time derivative (expressed in the inertial frame). Let define ν 2 = p q r T as the angular velocity of the body-fixed frame with respect to the earth-fixed frame expressed in the body-fixed frame (from now on: body-fixed angular velocity). The vector η2 does not have a physical interpretation and it is related to the body-fixed angular velocity by a proper Jacobian matrix:

ν 2 = J k,o (η 2 ) η2 . (2)
Figure 1 shows the defined frames and the elementary motions.

B. Dynamics

The rigid body dynamics of a quadrotor, in matrix form, is given by:

M RB ν + C RB (ν)ν + τ v,W + g RB (R B I ) = τ v , (3)
where

τ v = τ T 1 τ T 2 T , the vector τ 1 = X Y Z T
collects the linear forces acting on the rigid body expressed in a body-fixed frame, while

τ 2 = K M N T collects
the moments acting on the rigid body expressed in a bodyfixed frame.

The matrix M RB is constant, symmetric and positive definite, i.e.,

Ṁ RB = O, M RB = M T RB > O.
Its unique parametrization takes the form:

M RB = mI 3 -mS(r b C ) mS(r b C ) I O b ,
where

I 3 is the (3 × 3) identity matrix, r b C ∈ R 3
is the center-of-mass position expressed in body-fixed frame and I O b is the inertia tensor expressed in the body-fixed frame.

Notice that (3) can be simplified if the origin of the bodyfixed frame is chosen coincident with the central frame, i.e., r b C = 0 and I O b is a diagonal matrix.

The gravity force, acting in the center of mass r b C , is represented in body-fixed frame by:

g RB (R B I ) = mR B I g I 0 3×1
where g I = 0 0 9.81

T .

The term τ v,W ∈ R 6 represents external disturbances such as the wind; its effect on the vehicle is modeled as a constant disturbance in the earth-fixed frame that is further projected onto the vehicle-fixed frame. To this purpose, let define as γ v,W ∈ R 6 a vector of constant parameters; the current disturbance in the vehicle-fixed frame can be modelled as

τ v,W = Φ v,W (R I B )γ v,W = R B I O 3×3 O 3×3 R B I γ v,W (4)
where the (6 × 6) regressor matrix Φ v,W simply expresses the force/moment coordinate transformation between the two frames.

It is possible to rewrite eq. ( 3) by exploiting the linearity in the parameters as:

Φ v ( ν, ν, R B I )γ v = τ v (5)
where γ v ∈ R 16 is the vector of the dynamic parameters collecting the mass (1 parameter), the first moment of inertia (3 parameters), the inertia tensor (6 parameters) and the 6 elements of the disturbance γ v,W . The same equation may be easily rewritten with respect to the variables expressed in the inertial frame η, η, η following the guidelines of, e.g., [START_REF] Siciliano | Robotics: modelling, planning and control[END_REF].

In the following, the term Φ xy ∈ R 2×16 will denote the first two rows while φ z ∈ R 1×16 the third one of the regressor matrix expressed in the inertial frame. Following the guidelines of [START_REF] Gautier | Direct calculation of minimum set of inertial parameters of serial robots[END_REF], well established in robotics [START_REF] Antonelli | A systematic procedure for the identification of dynamic parameters of robot manipulators[END_REF], it is possible to further elaborate the regressor and classify the parameters among the sets: unidentifiable, identifiable alone and identifiable in linear combination. As an example, the body is affected by a vertical force caused by both the gravity and the wind; those effects can not be separated and the corresponding parameters will be identifiable in linear combination only. For sake of space, in this work the details are omitted; the controller tested in the following, in fact, will only consider the parameters that affect the steady state error.

C. Thrust

Quadrotors are equipped with 4 thrusters aligned along the body-fixed z axis with position p b t,i ∈ R 3 , and each of them provide each a force and a moment according to

f i = bω 2 t,i τ t,i = dω 2 t,i for i = 1, . . . , 4
where ω t,i is the angular velocity of the i th rotor, b and d are the thrust and drag coefficients. Figure 2 reports the common motor position with relevant variables. Notice that the body-fixed frame is positioned in the geometric center of the vehicle, i.e., in the intersection between the two thrusters axes.

x y z O x b y b z b O b f 1 f 2 f 3 f 4 1 2 3 4 ω t,1 ω t,2 ω t,3 ω t,4 τ t,1 τ t,2 τ t,3 τ t,4
l Fig. 2. Thrusters related variables It holds

τ 1 = 0 0 4 i=1 f i T .
By initially assuming that the center of gravity is coincident with the geometric center of the quadrotor and with a proper numeration of the thrusters (the first is on the positive x b axis and the numeration follows counter-clockwise looking from above), it holds:

τ 2 =   l(f 2 -f 4 ) l(f 1 -f 3 ) -τ t,1 + τ t,2 -τ t,3 + τ t,
 p b t,i -   0 0 r C,z     ×   0 0 f i   = p b t,i ×   0 0 f i   ∀ i
and thus the displacement along z b of the center of gravity does not modify the moment contribution. A value r C,x = 0 and r C,y = 0 introduces distortion effects along the pitch and roll directions respectively. After few computations it holds

M = l(f 1 -f 3 ) + r C,x (f 1 + f 3 ), K = l(f 2 -f 4 ) + r C,y (f 2 + f 4 ).
It is finally possible to write the mapping from the angular velocities to the force-torque at the vehicle:

    Z K M N     = B v     ω 2 t,1 ω 2 t,2 ω 2 t,3 ω 2 t,4     (6)
with

Bv =    b b b b 0 b(l + rC,y) 0 -b(l -rC,y) b(l + rC,x) 0 -b(l -rC,x) 0 -d d -d d    .

III. QUADROTOR ADAPTIVE CONTROL

The thrusters velocities may be assumed as the control input for the quadrotor control problem. The dynamics of the low level motor controller, in fact, can be typically neglected with respect to the vehicle dynamics; thus, we can assume u = ω 2 t,1 ω 2 t,2 ω 2 t,3 ω 2 t,4

T . Since the controller will output the desired force Z c K c M c N c T at the vehicle, the control input is computed by

u = B -1 v     Z c K c M c N c     (7)
where B -1 v ∈ R 4×4 is the inverse of ( 6)

B -1 v =          l -r C,x 4bl 0 1 2bl - l -r C,x 4dl l -r C,y 4bl
1 2bl 0 l -r C,y 4dl l + r C,x 4bl 0 - 1 2bl - l + r C,x 4dl l + r C,y 4bl - 1 2bl 0 l + r C,y 4dl          .
It is interesting to evaluate what happens if the mapping from the desired forces to the thrusters velocities is computed with the estimated mapping (r C ), while the effective mapping is physically related to r C :

    Z K M N     = B v | rC B -1 v rC     Z c K c M c N c     (8)
i.e.,

    Z K M N     =        1 0 0 0 rC,y 2 1 0 br C,y 2d rC,x 2 0 1 - br C,x 2d 0 0 0 1            Z c K c M c N c     (9)
where the terms non belonging to the unitary matrix represent a coupling effect that may arise if the center of mass is wrongly estimated or neglected. Also, force along z b and moment around z b are not affected by a wrong estimation of the center of mass and thus Z = Z c and N = N c . In the following, an adaptive control law for quadrotor position and yaw regulation will be developed by taking into account persistent external disturbances and unknown center of mass position. The assumption that the roll and pitch dynamics are faster than the position one is made. Figure 3 sketches the control loop.

A. Altitude control

Let define z = z dz ∈ R, s z = ż + λ z z ∈ R with λ z > 0 and γv = γ v -γv with the hat symbol denoting the estimate of the corresponding variable.

We consider a scalar Lyapunov candidate function V > 0:

V (s z , γv ) = m 2 s 2 z + 1 2 γT v K γ,z γv η 1d ψ d φ d , θ d Z c pos or   Kc Mc Nc   B -1 v u motors w 2 t,i B v    Z K M N    τ v,W η plant
V (s z , γv ) = s z mz d -mz + mλ z ż -γT v K γ,z γv = s z (φ z γ v -cos φ cos θZ) -γT v K γ,z γv in which φ z dependency is φ z (z d + λ z ż, η, R B I ). V is made negative semidefinite by selecting Z = 1 cos φ cos θ (φ z γv + k pz s z ) γv = K -1 γ,z φ T z s z yielding V (s z , γv ) = -k p,z s 2 z . (10)
We can now prove the system stability in a Lyapunovlike sense using Barbalats Lemma [START_REF] Khalil | Nonlinear Systems[END_REF]. Since V (s z , γv ) is lower bounded, V (s z , γv ) ≤ 0 and V (s z , γv ) is uniformly continuous, then V (s z , γv ) → 0 as t → ∞ and thus s z → 0.

As usual in adaptive control technique, we cannot prove asymptotic stability of the whole state, since γv is only guaranteed to be bounded. It is interesting to implement a simpler version of the controller aimed at compensating the sole persistent dynamic terms, i.e., those terms preventing a null steady state error, yielding:

Z = 1 cos φ cos θ (γ z + k pz s z ) (11) γz = k -1 γ,z s z ( 12
)
in which γ z ∈ R embeds the joint contribution of the gravity and the vertical wind effects.

B. Horizontal position control

Let us recall the elementary rotation around z as

R z = cos(ψ) sin(ψ) -sin(ψ) cos(ψ)
and define as ηxy = x dx y dy T ∈ R 2 and s xy = ηxy + λ xy ηxy ∈ R 2 with λ xy > 0.

We consider a Lyapunov candidate function V > 0:

V (s xy , γv ) = 1 2 ms T xy s xy + 1 2 γT v K γ,xy γv
whose time derivative is given by

V = ms T xy ṡxy -γT v K γ,xy γv = s T xy mη d,xy -mη xy + mλ xy ηxy -γT v K γ,xy γv = s T xy Φ xy γ v -ZR T z c φ s θ -s φ -γT v K γ,xy γv with Φ xy (η d,xy + λ xy ηxy , η, R B I ) ∈ R 2×16 .
This can be made negative semidefinite by selecting the virtual inputs φ and θ to solve

c φ s θ -s φ = 1 Z R z (Φ xy γv + k p,xy s xy ) , γv = K -1 γ,xy Φ T xy s xy .
Notice that the necessity to integrate the dynamic parameters with two different adaptive laws for the horizontal and vertical controllers is consistent. Each estimate, in fact, is used for its respective controller. In both cases, Lyapunov theory implies boundedness, but not convergence to zero, of the errors. In case of persistent excitation, the two estimates should converge to the same values. Due to the lack of space, further discussion is out of the scope of this paper.

In this case too, by considering only the persistent dynamic terms, the controller reduces to a simple

c φ s θ -s φ = 1 Z R z γxy + k p,xy s xy (13)
γxy = k -1 γ,xy s xy (14)
where γxy ∈ R 2 represents the sole wind effect supposed constant in the inertial frame, and φ and θ can be easily computed yielding to φ d and θ d to be sent to the orientation controller.

C. Orientation control

The orientation control receives as input the desired roll, pitch and yaw; the formers are obtained by the position control equations. Notice that, in this case, it is necessary to explicitly consider the presence of a coupling effect among the desired and obtained forces as shown in eq. ( 8):

K = K c + rC,y 2 Z c + br C,y 2d N c , M = M c + rC,x 2 Z c - br C,x 2d N c , N = N c .
It is worth noticing that neither the altitude nor the yaw control loop are affected by rC . The convergence to a steady state value for Z c and N c can thus be assumed, in addition, in absence of external moment disturbance along z b , at steady state N c = 0. In any case, roll and pitch control can be

K = K c + 1 2 Z c + b d N c rC,y M = M c + 1 2 Z c - b d N c rC,x .
The disturbance value is unknown and its effect may be compensated by resorting to several adaptive control laws well known in the literature.

D. Center of mass estimation

In case a simple PD control is used for pitch and roll control the steady state error is experienced. This effect can be counteracted linking the roll-pitch error to a proper integral estimate of the center of mass position according to:

ṙC,x ṙC,y = -k rC θ d -θ φ d -φ (15)
to be used in eq. [START_REF] Fink | Planning and control for cooperative manipulation and transportation with aerial robots[END_REF].

IV. SIMULATIONS

Numerical simulations in a wide range of operative conditions have been run in order to verify the effectiveness of the proposed adaptive controller. Due to lack of space the sole reduced controller will be reported here with the aim to prove that the persistent terms have been effectively compensated resorting to a minimal set of dynamic parameters whose number is 5. The sampling time of sensors and controller has been set to T = 1 ms; Tables I and II report the parameters used in simulation and, when applicable, the initial values used by the controller.

The simulation is run with a constant disturbance:

γ v,W = 0.5 0.6 0 0 0 0 T [N,Nm]
and by requiring a regulation displacement of 1 m in the 3 directions x, y and z and 20 deg in yaw.

In the following, the controller will be compared with its non adaptive version; in the latter, despite the name, the sole altitude part of the controller will benefit from the adaptive action. Comparison with [START_REF] Palunko | Agile load transportation. Safe and efficient load manipulation with aerial robots[END_REF] is not possible since the latter does not consider external disturbances, other than displacement of the center of mass that is not an external disturbance, and thruster mapping correction.

Figure 4 reports the norm of the position (top) and yaw (bottom) errors by applying the proposed control law (bluesolid line) and its non adaptive version (red-dashed line). As expected, the position error of the adaptive law goes to zero while the non adaptive version suffers from the presence of the external disturbance. Both the yaw errors goes to zero due to the mathematical properties of the quadrotor. Figure 7, finally, reports the parameters used in the simulation (gray) together with the estimated one (blue). The steady state values is always the true one. This is due to the tautological persistent excitation arising when compensating for the sole persistent terms. V. CONCLUSIONS Adaptive control for quadrotors is addressed in this paper. Convergence of the position and yaw error to zero is guaranteed in presence of constant external force/moment disturbance and unknown dynamic parameters, in particular for what concerns the center-of-mass position. A Lyapunovbased stability analysis has been used to design the controller, and results of numerical simulations are reported to validate it.

Fig. 1 .

 1 Fig. 1. Frames and elementary vehicle's motion

4 

 4 Let us now consider a center of gravity not coincident with the origin of the body-fixed frame, i.e. r b C = r C,x r C,y r C,z T = 0 T . Due to the distributive property over the addition of the cross product, it is possible to individually consider the displacement components. It is easy to verify that 

Fig. 3 .

 3 Fig. 3. Sketch of the control loop

Fig. 4 .Figure 5

 45 Fig. 4. Norm of the position (top) and yaw (bottom) errors by applying the proposed control law (blue-solid line) and its non adaptive version (reddashed line) Figure 5 reports the roll and pitch angles. On the top the desired (gray) and real (blue) values for the adaptive version while on the bottom the desired (gray dashed) and real (red dashed) values for the non adaptive simulation. Both controllers require a certain inclination at steady state to counteract for the wind.

Fig. 5 .

 5 Fig. 5. Roll and pitch angle; on the top the desired (gray) and real (blue) values for the adaptive version while on the bottom the desired (gray dashed) and real (red dashed) values for the non adaptive simulation Figure 6 reports the control effort in terms of force along z b (top) and moments (bottom) for the first seconds with

Fig. 6 .

 6 Fig. 6. Force along z b (top) and moments (bottom) by applying the proposed control law (blue-solid line) and its non adaptive version (reddashed line)

Fig. 7 .

 7 Fig. 7. Simulated parameters (gray) and estimated ones (blue)


Member, IEEE Dongjun J Lee
email: djlee@snu.ac.kr.

Member, IEEE, Hyoung Antonio Franchi
email: anto-nio.franchi@tuebingen.mpg.de.

Member, IEEE Il Son

Student Member, IEEE Ha Changsu

Member, IEEE Heinrich H Bülthoff

Member, IEEE Paolo Robuffo Giordano

C Ha

Semi-Autonomous Haptic Teleoperation Control Architecture of Multiple Unmanned Aerial Vehicles

Keywords: haptic feedback, multiagent control, passivity, teleoperation, unmanned aerial vehicles

We propose a novel semi-autonomous haptic teleoperation control architecture for multiple unmanned aerial vehicles (UAVs), consisting of three control layers: 1) UAV control layer, where each UAV is abstracted by, and is controlled to follow the trajectory of, its own kinematic Cartesian virtual point (VP); 2) VP control layer, which modulates each VP's motion according to the teleoperation commands and local artificial potentials (for VP-VP/VP-obstacle collision avoidance and VP-VP connectivity preservation); and 3) teleoperation layer, through which a single remote human user can command all (or some) of the VPs' velocity while haptically perceiving the state of all (or some) of the UAVs and obstacles. Master-passivity/slave-stability and some asymptotic performance measures are proved. Experimental results using four custom-built quadrotor-type UAVs are also presented to illustrate the theory.

I. INTRODUCTION

Due to the absence of human pilots on-board, unmanned aerial vehicles (UAVs) can realize many powerful aerospace applications with reduced cost/danger and possibly higher performance than the conventional pilot-driven aerial vehicles: surveillance and reconnaissance, fire-fighting and rescue, remote sensing and exploration, pesticide spraying and geophysical survey, logistics and payload transport, and ad-hoc communication gateway, to name just few. See [START_REF] Valavanis | Advances in Unmanned Aerial Vehicles: State of the Art and the Road to Autonomy[END_REF], [2]. Deploying multiple UAVs will further enhance these applications by infusing them with the benefits of multi-agent systems (e.g., better performance via cooperation such as higher payload transport and faster domain coverage; better affordability than a single/bulky system; robustness against single point failures, etc.). On the other hand, many real UAV applications take place in environments, which are unstructured, uncertain and not precisely known a priori. For such cases, fully-autonomous control of the UAVs is typically infeasible/impossible, and, instead, to impose human's intelligence on the task to cope with such uncertainty, some teleoperation of their behaviors is desired, if not absolutely necessary 1 . Now, suppose that a large number of UAVs is presented to a human user and s/he is required to teleoperate their motions all at once. This would define a daunting task for the human user, since we humans can control well only a small number of degrees-of-freedom (DOFs) simultaneously, yet, such many UAVs are characterized by a large number of DOFs. However, if we examine many practical aerospace applications deploying multiple UAVs, we can also see that, very often, the task for each UAV can be split into a component, that is rather simple and mathematically well-defined (e.g., maintaining relative distance, avoid collision/obstacles, etc.) and another component, that is mathematically obscure and requires complex intelligent/cognitive information processing (e.g., how to drive UAVs in the presence of uncertainty; whether to proceed/stop when obstacles appear, etc.).

With this distinction in mind, in this paper, we propose a novel semi-autonomous haptic teleoperation control framework for multiple UAVs, which enables a single remote human user to stably teleoperate the overall (abstracted) motion of the multiple UAVs with some useful haptic feedback, while the UAVs are reacting autonomously among themselves and against local obstacles so as to render themselves collectively as a flying teleoperated deformable object. More specifically, our semi-autonomous teleoperation control architecture consists of the following three control layers (see Fig. 1):

1) UAV control layer, which enforces each UAV to (unilaterally) track the trajectory of its own first-order kinematic Cartesian virtual point (VP). For this, we assume availability of some reasonably-good trajectory tracking control law for each UAV (e.g., [START_REF] Koo | Output tracking control design of a helicopter model based on approximate linearization[END_REF], [START_REF] Mahony | Robust trajectory tracking for a scale model autonomous helicopter[END_REF], [START_REF] Aguiar | Trajectory-tracking and path-following of underactuated autonomous vehicles with parametric modeling uncertainty[END_REF], [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF], [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF], [START_REF] Pflimlin | Position control of a ducted fan vtol uav in crosswind[END_REF], [START_REF] Oishi | Switched nonlinear control of a vstol aircraft[END_REF], or that presented in Sec. II-B), which then allow us to abstract each UAV by their kinematic VP for the purpose of their semi-autonomous teleoperation control design, while bypassing the low-level control issues of UAVs (e.g., control under-actuation [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF], [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF]);

2) VP control layer, which modulates the motion of the multiple VPs in such a way that, as a whole, in a distributed manner (i.e., each VP is sensing/communicating only with their own neighboring VPs on a certain timeinvariant connectivity (or information) graph G), they collectively behave as a multi-nodal flying deformable object, whose shape autonomously deforms according to local artificial potentials (designed for VP-VP/VPobstacle collision avoidance and VP-VP connectivity maintenance), while whose bulky motion is driven by the teleoperation (velocity) command received from the master side; 3) teleoperation layer, which enables a remote human user to tele-drive some (or all) of the VPs (i.e., control set N c ), while haptically perceiving the state of some (or all) of the real UAVs (i.e., sensing set N s ). For this, passive set-position modulation (PSPM [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF]) is adopted due to its implementation flexibility (e.g., can accommodate master-slave kinematic/dynamic dissimilarity and various forms of haptic feedback signals), guarantee of passivity (i.e., interaction stability with human users), and less conservative passifying action (thus, better performance). Although there are numerous results for the single UAV motion control (e.g. [START_REF] Koo | Output tracking control design of a helicopter model based on approximate linearization[END_REF]- [START_REF] Oishi | Switched nonlinear control of a vstol aircraft[END_REF]) and some results for the single UAV haptic teleoperation (e.g. [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF], [START_REF] Lam | Haptic feedback in uninhabited aerial vehicle teleoperation with time delay[END_REF], [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF], [START_REF] Masone | Interactive planning of persistent trajectories for human-assisted navigation of mobile robots[END_REF]), it would be fair to say that the problem of haptic teleoperation of multiple UAVs starts being considered fairly recently.

The framework of [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], which was later applied to UAVs in [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], utilizes passive decomposition [START_REF] Lee | Passive decomposition of mechanical systems with coordination requirement[END_REF], [START_REF] Lee | Passive decomposition approach to formation and maneuver control of multiple rigid bodies[END_REF] to precisely maintain the formation shape (i.e., shape system) among the agents, while their overall motion (i.e., locked system) is teleoperated. These results [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], however, demand allto-all communication among the agents and also are limited to master-slave position-position teleoperation, which is not so suitable for the UAV teleoperation, since the master workspace is bounded, yet, that of the slave UAVs is not (i.e., kinematic dissimilarity [START_REF] Lee | Feedback r-passivity of lagrangian systems for mobile robot teleoperation[END_REF]). Passive decomposition, extended to nonholonomic systems [START_REF] Lee | Passive decomposition and control of nonholonomic mechanical systems[END_REF], was also used in [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF] for the teleoperation of multiple wheeled mobile robots. This work [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF] has some similarity with the current work (e.g., master-position/slave-velocity teleoperation via PSPM), yet, still requires all-to-all communication and does not address collision-avoidance/connectivity-preservation among the slave robots.

To our knowledge, the works of [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF] are the very first results on the haptic teleoperation of multiple UAVs, that do not require all-to-all communication among the UAVs (i.e., distributed) and also consider the issues of collision and connectivity among them. Here, [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF] is the conference version of the current manuscript, while [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], later expanded in [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], considers the switching leader-follower information topology among the second-order dynamic VPs with constant time-delay between the master and the leader UAV. Subsequent developments/improvements ensuing these [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF] are as follows: 1) in [START_REF] Franchi | Modeling and control of UAV bearing-formations with bilateral high-level steering[END_REF], a partially-decentralized shared control system for multi-UAVs is proposed, that is based only on bearing (angular) measurements locally obtainable from camera-like sensors; 2) in [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF], a decentralized approach, that can enforce global connectivity (e.g., for steady information flow among the UAVs) in the presence of graph switching and teleoperation, is presented; and 3) in [START_REF] Son | Human-centered design and evaluation of haptic cueing for teleoperation of multiple mobile robots[END_REF], by using the control framework proposed in [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF] and in this paper, the impact and effectiveness of different haptic feedback for multi-UAV haptic teleoperation is studied from a perceptual point of view.

Differently from the works mentioned above [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF]- [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF], our semi-autonomous haptic teleoperation control architecture, proposed first in [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF] and detailed/completed in this paper, possesses the following properties: 1) the information flow (i.e., connectivity graph) among the UAVs is distributed (cf. [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF], [START_REF] Franchi | Modeling and control of UAV bearing-formations with bilateral high-level steering[END_REF]) and their collective shape can reactively deform according to the external environment (cf. [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF]); 2) master-passivity/slave-stability (with first-order kinematic VPs) is enforced, which is likely less conservative than master-passivity/slave-passivity of other results using the second-order dynamic VPs (cf. [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF]); 3) any forms of haptic feedback signal can be adopted without jeopardizing master-passivity/slavestability even in the presence of communication unreliability2 (cf. [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF]); and 4) the human user can freely choose any "control set" N t and any "sensing set" N s from multiple UAVs for tele-control/sensing depending on task objectives and conditions (cf. [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF]). Our semi-autonomous control architecture has also served as the foundation for some of those subsequent results (e.g., kinematic VPs and flexible PSPM with masterpassivity/slave-stability: [START_REF] Franchi | Modeling and control of UAV bearing-formations with bilateral high-level steering[END_REF], [START_REF] Son | Human-centered design and evaluation of haptic cueing for teleoperation of multiple mobile robots[END_REF], [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF]).

A portion of this paper was presented in [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF]. The current version has been substantially revised from [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF], particularly with: 1) full experiment (i.e., with real UAVs) as compared to the only semi-experiment (i.e., with simulated UAVs) in [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF] (Sec. III); 2) complete explanation on implementing the UAV control layer, which was only alluded in [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF] (Sec. II-B); and 3) whole new introduction, improved organization and significantly expanded explanations of technical results/details. Some high-level description of our semi-autonomous architecture was also reported in [START_REF] Franchi | Shared control: Balancing autonomy and human assistance with a group of quadrotor uavs[END_REF], yet, without technical details, which are fully provided in this paper.

The rest of this paper is organized as follows. Sec. II introduces some preliminary materials and details the three control layers: UAV control layer in Sec. II-B; VP control layer in Sec. II-C; and teleoperation layer in Sec. II-D. Sec. III presents experimental results using four custom-built quadrotor-type UAVs with hardware/software details. Sec. IV summarizes the paper with some comments on future research directions.

II. SEMI-AUTONOMOUS TELEOPERATION CONTROL ARCHITECTURE

A. Slave UAVs and Master Haptic Device

Let us consider N UAVs, whose 3-DOF Cartesian positions are denoted by x i ∈ ℜ 3 , i = 1, 2, ..., N . Here, we are interested in the case where a single human user teleoperates the Cartesian positions x := [x 1 ; x 2 ; ...; x N ] ∈ ℜ 3N of the N UAVs simultaneously. For this, we do not require the UAVs to be of a specific type. We rather allow them to be of any types (e.g., swarm of heterogeneous UAVs) as long as a reasonablyperforming trajectory tracking control exists for them so that we can drive each x i to faithfully track a smooth reference trajectory. See Sec. II-B.

One class of such UAVs, that possesses a well-performing trajectory tracking control, is so called vectored-thrust (or thrust-propelled) UAVs [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF], whose 6-DOF (under-actuated) dynamics in SE(3) is given by

m i ẍi = -ρ i R i e 3 + m i ge 3 + δ i (1)
J i ẇi + S(w i )J i w i = γ i + ζ i , Ṙi = R i S(w i ) (2)
where m i > 0 is the mass, x i ∈ ℜ 3 is the Cartesian center-ofmass position w.r.t. the NED (north-east-down) inertial frame (with e 1 , e 2 , e 3 representing N, E and D directions), ρ i ∈ ℜ is the thrust control input along the body-frame e 3 , R i ∈ SO(3) is the rotational matrix describing the body NED frame of UAV w.r.t. to the inertial NED frame, w i ∈ ℜ 3 is the angular rate of the UAV relatively to the inertial frame represented in the body frame, J i ∈ ℜ 3×3 is the UAV's inertia matrix w.r.t. the body frame, g is the gravitational constant, γ i ∈ ℜ 3 is the attitude torque control input, δ i , ζ i ∈ ℜ 3 are the aerodynamic perturbations, and S(⋆) : ℜ 3 → so(3) is the skew-symmetric operator defined s.t. for a, b ∈ ℜ 3 , S(a)b = a × b. Some examples of this vectored-thrust UAV include: autonomous helicopters [START_REF] Mahony | Robust trajectory tracking for a scale model autonomous helicopter[END_REF], [START_REF] Tanaka | Wireless visionbased stabilization of indoor microhelicopter[END_REF], VTOL aircraft [START_REF] Oishi | Switched nonlinear control of a vstol aircraft[END_REF], ducted-fan UAVs [START_REF] Pflimlin | Position control of a ducted fan vtol uav in crosswind[END_REF] and, quadrotors [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF], which we use for the experiment in Sec.

III.

Let us also consider a 3-degree-of-freedom (DOF) nonlinear Lagrangian haptic device as modeled by [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF], [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF] M (q)q + C(q, q) q = τ + f

where q ∈ ℜ 3 is the configuration, M (q) ∈ ℜ 3×3 is the positive-definite/symmetric inertia matrix, C(q, q) ∈ ℜ 3×3 is the Coriolis matrix, and τ, f ∈ ℜ 3 are the control and human forces, respectively. It is well-known that this haptic device is (energetically) passive: ∀T ≥ 0, ∃d ∈ ℜ s.t.,

T 0 [τ + f ] T qdt = κ(T ) -κ(0) ≥ -κ(0) =: d 2
which can be easily shown by using that Ṁ -2C is skewsymmetric [START_REF] Spong | Robot modeling and control[END_REF].

Our goal is then to enable a single remote user to teleoperate N UAVs' Cartesian motions x := [x 1 ; x 2 ; ...; x N ] ∈ ℜ 3N via the single 3-DOF master haptic device (3) simultaneously, while providing the user with some useful haptic feedback to convey information of the N UAVs' state and their surrounding environments. There are several interesting aspects/challenges to achieve this: 1) large slave DOF: humans can usually control well only a small number of DOFs at the same time (e.g., 3-DOF master device (3)), yet, the slave N UAVs possess a large number of DOFs; 2) kinematic/dynamic dissimilarity [START_REF] Lee | Feedback r-passivity of lagrangian systems for mobile robot teleoperation[END_REF]: usual master device (3) has a bounded workspace with full-actuation (e.g., joystick), yet, the UAVs' workspace is unbounded (e.g., E(3)) and their dynamics (e.g., (1)-( 2)) typically not so favorable to control as the fullactuated master device (e.g., under-actuated 3-DOF translation dynamics (1) with 1-DOF thrust control ρ i ); and 3) control distribution among UAVs: information flow (either through communication or sensing) among the UAVs is desired to be distributed (i.e., each UAV requires information only from their neighbors and possibly from the master site), particularly when the number of UAVs is large.

To address these challenges, we propose a novel semiautonomous haptic teleoperation control architecture for multiple UAVs, consisting of the three control layers, UAV control layer, VP control layer, and PSPM-based teleoperation layer, each to be detailed in the following three subsections.

B. UAV Control Layer

In this paper, we are interested in teleoperating the Cartesian motions x = [x 1 ; x 2 ; ..., x N ] ∈ ℜ 3N of the N UAVs, whose dynamics (e.g., [START_REF] Valavanis | Advances in Unmanned Aerial Vehicles: State of the Art and the Road to Autonomy[END_REF]), yet, is typically too complicated to be directly handled with by the standard teleoperation techniques (e.g., under-actuation of ( 1)). It is also desirable in many cases to "hide" this complex underlying dynamics of the UAVs from the human user so that s/he can focus more on the high-level teleoperation of the multiple UAVs without being distracted to simultaneously taking care of their low-level dynamics.

To circumvent this issue, we endow each UAV with a 3-DOF Cartesian VP (virtual point), p i ∈ ℜ 3 . The human user will then teleoperate these N VPs instead of the real UAVs, while the real UAV (i.e., x i ) is tracking its own VP (i.e., p i ). See Fig. 1. Abstracting each UAV by their VP and formulating the semi-autonomous teleoperation objectives on these (simpler) N VPs, we can greatly simplify the design of the VPcontrol layer (Sec. II-C) and of the teleoperation layer (Sec. II-D), while encapsulating the issue of the UAVs' complex low-level dynamics within the UAV control layer. This also implies that our semi-autonomous teleoperation architecture is applicable to any types of (possibly heterogeneous) mobile robots (e.g., humanoids, unicycles) as long as they possess some adequately-functioning tracking controller to follow their own VPs.

The vectored-thrust UAVs (1)-( 2) assume many of such well-performing tracking control laws (e.g., [START_REF] Mahony | Robust trajectory tracking for a scale model autonomous helicopter[END_REF], [START_REF] Pflimlin | Position control of a ducted fan vtol uav in crosswind[END_REF], [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF], [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF]). Although one of these schemes can certainly be used, to facilitate readers' implementation of our framework, here, we present a simple tracking control law for (1)-( 2), which turns out to be fairly reliable and robust during our experiments in Sec. III and also many other demonstrations performed at the authors' institutions using quadrotor UAVs.

Our control law is based on the natural decoupling property of ( 1)-( 2), that is, the attitude dynamics (2) is independent from the translation dynamics [START_REF] Valavanis | Advances in Unmanned Aerial Vehicles: State of the Art and the Road to Autonomy[END_REF]. We then design our controller to have the following inner-outer loop structure: 1) a slower outer-loop position tracking controller is designed for (1) to drive x i to track p i , while specifying thrust and attitude commands; whereas 2) a faster inner-loop attitude controller is designed for (2) to attain the attitude commands given from the outer-loop. In the following derivations of the controller, for notational convenience, we omit UAV's index i in ( 1)-( 2).

First, let η := [φ, θ, ψ] T ∈ ℜ 3 be the RPY Euler angle representation of the rotation matrix R, with φ, θ, ψ being respectively the roll, pitch and yaw angles of the UAV along the NED-directions. We can then rewrite the attitude dynamics (2) using η s.t.

Σ 1 : J ẇ = -S(w)Jw + γ + ζ η = T (η)w (4)
where T (η) ∈ ℜ 3×3 is the transformation matrix from w ∈ so(3) to the Euler angle rates η. The translational dynamics (1) can also be written by using η s.t.

Σ 2 : mẍ = mge 3 -ρR(η)e 3 + δ (5)
where

R(η) = R e3 (ψ)R e2 (θ)R e1 (φ)
, with R ei being the elementary rotation matrix about the e i -axis [START_REF] Spong | Robot modeling and control[END_REF].

Notice that the attitude dynamics Σ 1 is independent from Σ 2 (and ρ, which will contain the control action for x), while the reverse does not hold, i.e., Σ 2 depends on Σ 1 due to the term R(η). The goal of the controller is then to make use of the four control inputs, γ ∈ ℜ 3 and ρ ∈ ℜ, to separately control the UAV's position x = (x 1 , x 2 , x 3 ) T to track a (smooth) VP's reference trajectory p = (p 1 , p 2 , p 3 ) T and the yaw angle ψ to (possibly time-varying) target value ψ d .

Let us start with the position controller first. For this, using

R(η) = R e3 (ψ)R e2 (θ)R e1 (φ), we can write (5) as R T e3 (ψ)mẍ = mge 3 -R e2 (θ)R e1 (φ)ρe 3 + R T e3 (ψ)δ (6) whose last row reads s.t. mẍ 3 = mg -cos(φ) cos(θ)ρ + δ 3 with δ = (δ 1 , δ 2 , δ 3 ) T . This then suggests the following thrust control ρ = - m cos φ cos θ [-g+ p3 +k d P ( ṗ3 -ẋ3 )+k p P (p 3 -x 3 )] (7)
which ensures local exponential stability of p 3x 3 , as long as the system is away from the singularity cos φ cos θ = 0 and δ 3 = 0.

On the other hand, the first two rows of ( 5) are given by

m ẍ1 ẍ2 = -ρ cos φ cos ψ sin ψ cos φ sin ψ -cos ψ =:Q(φ,ψ)∈ℜ 2×2 sin θ sin φ + δ 1 δ 2
where Q is always analytically invertible as long as cos φ = 0. This then shows that (p 1x 1 , p 2x 2 ) will be locally exponentially stable (with δ 1 = δ 2 = 0), if the attitude controller (to be defined below) can attain pitch and roll commands θ d , φ d given by

sin θ d sin φ d = mQ -1 -ρ p1 + k d P ( ṗ1 -ẋ1) + kp P (p1 -x1) p2 + k d P ( ṗ2 -ẋ2) + kp P (p2 -x2) . (8)
Here, although [START_REF] Pflimlin | Position control of a ducted fan vtol uav in crosswind[END_REF] defines a nonlinear equation for φ, we found it works pretty well in practice to obtain (θ d , φ d ) while assuming Q be a function of (φ, ψ), particularly with fast enough attitude control servo-rate (e.g., 500Hz in Sec. III). Now, define the desired attitude set-point by

η d := [φ d , θ d , ψ d ],
where φ d , θ d are given above from the (outerloop) position controller [START_REF] Pflimlin | Position control of a ducted fan vtol uav in crosswind[END_REF], while ψ d can be set arbitrarily (e.g., ψ d = 0). Given η d , we then design the (fast/inner-loop) attitude regulation control as follows. First, differentiating the second row of (4), we can get η = T (η) ẇ + Ṫ (η)w which, by using the first row of (4), becomes

η = T (η)J -1 (-w × Jw + γ + ζ) + Ṫ (η)w.
We then choose the attitude regulation control s.t.

γ = JT -1 (η)(-k d A η + k p A (η d -η)) (9)
with which we have the following closed-loop dynamics:

ëη + [k d A + T J -1 S(w)Jw -Ṫ T -1 ] ėη + k p A e η = T J -1 ζ
with e η := ηη d , implying that, if ζ = 0, (e η , ėη ) will be locally exponentially stable, if (e η (0), ėη (0)) is small enough and k d A large enough. The presented controller ( 7)-( 9), although simple in its structure (i.e., easier to implement than [START_REF] Mahony | Robust trajectory tracking for a scale model autonomous helicopter[END_REF], [START_REF] Pflimlin | Position control of a ducted fan vtol uav in crosswind[END_REF], [START_REF] Hua | A control approach for thrust-propelled underactuated vehicles and its application to vtol drones[END_REF], [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF]), turns out to be sufficient for our experiments (and many on-site demonstrations) as evidenced in Sec. III. Thus, from now on, we assume that: 1) we have implemented for each UAV a trajectory tracking control to make x i to follow p i ; and 2) this control performs reasonably well, by keeping ||x ip i || and || ẋiṗi || small enough (||x|| 2 := x T x). How to control the motion of N VPs is then the subject of the next Sec. II-C.

C. Distributed VP Control Layer

We consider N first-order kinematic VPs to abstract each UAV, with their Cartesian position denoted by p i ∈ ℜ 3 (i = 1, 2, ..., N ). Our goal is to render these N VPs as a N -nodes flying deformable object in a distributed manner, so that their shape autonomously deforms reacting to the presence of obstacles, while their collective motion is tele-controlled by a single remote human user, with the information-flow among the VPs distributed. See Fig. 1. As long as the UAV control layer in Sec. II-B ensures ||x ip i || and || ẋiṗi || be small, these VPs' behaviors will be then faithfully duplicated among the real UAVs.

To describe how the VPs are connected (via communication or sensing) to form the N -nodes flying object, we define the undirected connectivity graph G, with the N VPs as its nodes and their connection (i, j) as its edges. Since G is undirected,

(i, j) ∈ E(G) iff (j, i) ∈ E(G)
, where E(G) is the edge set of G. We assume G is connected and also dense enough so that, with some suitably-defined inter-VP attractive/repulsive potentials on E(G), it can define the un-deformed shape of the N -nodes deformable object with no inter-VP separation or collision (e.g., rigid graph [START_REF] Aspnes | A theory of network localization[END_REF], [START_REF] Zelazo | Rigidity maintenance control for multi-robot systems[END_REF]). We also assume G to be time-invariant (e.g., no creation/elimination of edges over time). For some applications, a time-varying G may be useful (e.g., separation of N -nodes flying object to penetrate narrow passages and merge afterwards). See [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF], where such a time-varying G is achieved for second-order dynamic VPs under leader-follower connectivity graph.

We implement the following kinematic evolution of VP on each UAV: for the i th UAV,

ṗi (t) := u t i + u c i + u o i ( 10
)
where 1) u c i ∈ ℜ 3 embeds the inter-VP collision avoidance and connectivity preservation, as defined by

u c i := - j∈Ni ∂ϕ c ij (||p i -p j || 2 ) T ∂p i (11)
where ϕ c ij is a certain artificial potential function to create attractive action if ||p i -p j || is large, and repulsive action if ||p ip j || small (see Fig. 2), and

N i := {j | (j, i) ∈ E(G)}
i.e., the connectivity neighbors of the i th VP on G; 2) u o i ∈ ℜ 3 is the obstacle avoidance action as given by

u o i := - r∈Oi ∂ϕ o ir (||p i -p o r ||) T ∂p i ( 12
)
where O i is the set of obstacles of the i-th VP with p o r being the position of the r th obstacle in O i , and ϕ o ir is a certain artificial potential, which produces repulsive action if ||p ip o r || is small, smoothly converges to zero as ||p ip o r || → d, and stays zero for ||p ip o r || ≥ d, to make the effect of obstacles for each VP gradually emerge/disappear when they move closer/farther from the VP than d > 0 (see Fig. 2); 3) u t i ∈ ℜ 3 will contain the teleoperation command for the UAVs in the "control set" N t ⊂ {1, 2, ..., N }, to enable a remote human user to directly tele-drive the Cartesian velocity of these VPs in N t and, consequently, the collective velocity of the N -nodes flying deformable object (to be designed in Sec. II-D). Each UAV numerically integrates this kinematic evolution equation [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF] over the sampling time (e.g., 150Hz for Sec. III) to obtain (and track) the position p i of their own VP.

Here, both the inter-VP potential ϕ c ij and the VP-obstacle )), to allow for the rotational symmetry [START_REF] Ogren | Cooperative control of mobile sensor networks: Adaptive gradient climbing in a distributed environment[END_REF] of the N -VPs deformable flying object (see Sec. III-D); and 2) they rapidly increase when inter-VP/VP-obstacle collisions or inter-VP separation (e.g., with limited communication range) are impending to prevent that. For a design example of ϕ c ij , ϕ o ir , see [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF], [START_REF] Dimarogonas | Connectedness preserving distributed swam aggregation for multiple kinematic robots[END_REF]. The next Prop. 1 summarizes some key properties of the swarm behavior of the N VPs [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF], with ϕ c ij , ϕ o ir and bounded u t i . For that, define the total potential energy s.t.

V (t) := 1 2 N i=1 j∈Ni ϕ c ij (||p i -p j ||) + N i=1 r∈Oi ϕ o ir (||p i -p o r ||)
and also assume that ϕ c ij and ϕ o ir are constructed s.t.: 1) there exists a large enough M > 0 s.t. V (t) ≤ M implies inter-VP connectivity preservation and no inter-VP/VP-obstacle collisions (e.g., rigid graph [START_REF] Aspnes | A theory of network localization[END_REF], [START_REF] Zelazo | Rigidity maintenance control for multi-robot systems[END_REF]

N i=1 ṗi = N t • u t i , ∀t ≥ 0.
Proof: Here, we only provide a sketch of proof. Refer to [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF] for some more details. Differentiating the above V (t) with

N i=1 j∈Ni ∂ϕ c ij ∂pi ṗi = N i=1 j∈Ni ∂ϕ c ij ∂pj ṗj , we have V = N i=1 W T i (-W i + u t i ) ≤ N i=1 (-||W i || 2 + δ it ū||W i ||)
where [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF] for all VPs with u o i ≡ 0. The assumption (13) of Prop. 1 is mild, since it just rules out the practically improbable (e.g., zero-measure) situation, where, although V (t) is very large, none of the VPs can detect that, with all of their (also very large) forces, ∂ϕ c ij /∂p i and ∂ϕ o ir /∂p i , somehow exactly aligned with each other to make their sum nonetheless to be small. Prop. 1 also says that: 1) all the N VPs is guaranteed to be stable (i.e., bounded ṗi and no collisions/separations) for any bounded u t i regardless of whether u t i are heterogeneous among the VPs or applied to some or all of them; and 2) the more VPs implement u t i , the easier it would be to drive the N -VPs deformable flying object. In the next Sec. II-D, this Prop. 1 will allow us to enforce master-passivity/slave-stability of the total closed-loop teleoperation system, with passive set-position modulation (PSPM [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF]) robustly guaranteeing the boundedness of u t i and master-side passivity. See [START_REF] Cao | Distributed coordinated tracking via a variable structure approach -part ii: swarm tracking[END_REF], where a result similar to Prop. 1 was achieved by using a sliding mode control approach.

W i := -u c i -u o i . Now, suppose that V (t) ≥ M . Then, from (13), ∃ s ∈ {1, 2, ..., N } s.t. ||W s || ≥ ū( √ N t + δ st )/2. We can thus obtain V ≤ -i / ∈Nt∪{s} ||W i || 2 - ||W s || 2 + δ st ū||W s || - i∈Nt\{s} (||W i || 2 -ū||W i ||) ≤ -||W s || -δst ū 2 2 + N t ū2 4 ≤ 0,
Our usage of VPs is inspired by [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF]. Yet, instead of the second-order dynamic VPs in [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF], here, we choose the simpler first-order kinematic VPs [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF], since 1) we can significantly simplify/strengthen the VPs' swarm control design/analysis (e.g., stability of Prop. 1 valid for any bounded u t i applied to any VPs); and 2) we can avoid some performance-limiting aspects encountered with the dynamic VPs (e.g., operator's continuous exercising against system/control damping [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF]). As shown in Sec. II-D, this kinematic VP will also allow us to achieve master-passivity/slave-stability, which is less conservative than master-passivity/slave-passivity typical for the dynamic VPs [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF], [START_REF] Franchi | A passivity-based decentralized approach for the bilateral teleoperation of a group of uavs with switching topology[END_REF], [START_REF] Franchi | Bilateral teleoperation of groups of mobile robots with timevarying topology[END_REF], [START_REF] Secchi | Bilateral teleoperation of a group of UAVs with communication delays and switching topology[END_REF], [START_REF] Giordano | A Passivity-based decentralized strategy for generalized connectivity maintenance[END_REF], thus, would likely provide a sharper performance than achievable with the dynamic VPs. See [START_REF] Sarlette | Cooperative attitude synchronization in satellite swarms: a consensus approach[END_REF], [START_REF] Palafox | Bilateral teleoperation of a formation of nonholonomic mobile robots under constant time delay[END_REF] for other "kinematic" abstractions. Our usage of VPs (10) also shares a similarity with the multi-nodal deformable object modeling or distributed behavioral swarm modeling in computer graphics (e.g., [START_REF] Gibson | A survey of deformable modeling in computer graphics[END_REF], [START_REF] Reynolds | Flocks, herds, and schools: a distributed behavioral model[END_REF]), in which, yet, (useful/important) theoretical guarantees as obtained here (e.g, stability with no collision/separation of Prop. 1; masterpassivity/slave-stability of Th. 1) are typically missing.

D. PSPM-Based Teleoperation Layer

For the teleoperation layer, we utilize passive set-position modulation (PSPM [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF]), whose passifying action theoretically guarantees master-passivity/slave-stability of the closed-loop teleoperation system, while whose flexibility allows us to accommodate kinematic/dynamic dissimilarity between the master device (3) and the VPs [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF] and also various forms of haptic feedback. PSPM also exhibits better performance than other "time-invariant' teleoperation schemes (e.g. wave/PD) due to its less conservative "selective" passifying action. Our treatment on PSPM here is brief: see [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF], [START_REF] Huang | Implementation and experiments of passive setposition modulation framework for interactive robotic systems[END_REF] for more details.

First, to enable a remote human user to teleoperate the VPs in the control set N t , we define u t i (t) in (10) s.t., u t i (t) := λH[q(k)], ∀i ∈ N t [START_REF] Lee | Bilateral teleoperation of multiple cooperative robots over delayed communication networks: theory[END_REF] for t ∈ [t k , t k+1 ), where q(k) is the master configuration q(t) ∈ ℜ 3 received from the communication (e.g., Internet) by the i th UAV at the reception time t k , λ > 0 is to match different scales between q and ṗi , and H[•] is a continuoustime (BIBO) stable low-pass filter with a fast enough timeconstant to ensure u t i (t) be smooth (to obtain p for, e.g., ( 7)-( 8)) while tracking λq(k) quickly enough. By providing coupling between the VPs' velocity ṗi and master's position q, this control ( 14) allows us to address the issue of master-slave kinematic dissimilarity [START_REF] Lee | Feedback r-passivity of lagrangian systems for mobile robot teleoperation[END_REF] (i.e., mobile VPs with unbounded workspace; master device with bounded workspace).

On the other hand, to allow the remote user to tele-sense some (or all) of the UAVs and their surrounding obstacles, we design the haptic feedback y(t) ∈ ℜ 3 s.t.

y(t) := 1 λN s i∈Ns ( ẋi + u o i ) =:yi(t) (15)
where N s ⊂ {1, 2, ..., N } is the "sensing" set among the N UAVs with N s > 0 being its cardinality, ẋi is the i th UAV's velocity (1), and u o i is the i th VP's obstacle avoidance control [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF]. This y(t) is designed to allow the user: 1) to directly perceive the state of the real UAVs, thereby, completes the "information closed-loop" (i.e. master → VPs → UAVs → master) to overcome the unilaterality of the UAV control layer (e.g., prevent the user from tele-driving the VPs without knowing the UAVs are left behind due to, e.g., actuator failures) and also; 2) to tele-sense the presence of obstacles through their collective effects (i.e., (1/(λN s )) i∈Ns u o i ) on the VPs/UAVs in the sensing set N s .

Each UAV in the sensing set N s then sends its y i (t) (15) to the master site over communication network (e.g., Internet). Let us denote by y(k) := (1/(λN s )) i∈Ns y i (k) the haptic signal y(t) constructed in the master side at the reception time t k . We incorporate this y(k) into the teleoperation control τ in (3) s.t.: for t ∈ [t k , t k+1 ),

τ (t) := -B q -K f q -K(q -ȳ(k)) (16)
where B, K, K f ∈ ℜ 3×3 are the positive-(semi)definite diagonal gain matrices, and ȳ(k) is the PSPM-modulation of the received haptic feedback y(k) (to be defined below). This haptic control ( 16) is designed s.t.: 1) if the UAVs fleet reaches the commanded velocity, the user will perceive this steadystate UAVs' velocity haptically (via f in (3)) and/or visually (by seeing q); and 2) if the UAVs fleet approaches obstacles, the user will haptically perceive these obstacles through their collective action N i=1 u o i through the haptic feedback y. See item 2 of Th. 1.

At each t k , PSPM-modulation ȳ(k) in ( 16) is defined by

min ȳ(k) ||y(k) -ȳ(k)|| subj. E(k) = E(k -1) + D min (k -1) -∆ P (k) ≥ 0 that is, ȳ(k)
is chosen as close to y(k) as possible for performance (first line), yet, only to the extent permissible by the passivity constraint (second line). Here, E(k) ≥ 0 is the virtual energy reservoir (simulated in software); D min (k) :=

1 t k+1 -t k 3 i=1 b i (q i (k) -q i (k))
2 is (conservative) estimate of the (otherwise-wasted) energy dissipation via the damping B in [START_REF] Lee | Passive decomposition of mechanical systems with coordination requirement[END_REF], which is re-harvested into E(k), with b i > 0 being the i th diagonal element of B (i = 1, 2, 3), q i the i th element of q, and qi (k) and q i (k) respectively the maximum and minimum of q i (t) during [t k , t k+1 ); and, with ||x|| 2

A := x T Ax, ∆ P (k) := ||q(t k ) -ȳ(k)|| 2 K /2 -||q(t k ) -ȳ(k -1)|| 2 K /2
is the energy jump at t k , which is passified by choosing ȳ(k) to satisfy the passivity constraint (second line).

This PSPM is implemented only on the master side. Since the human operator usually keeps injecting energy into the master system, E(k) may keep increasing as well. To avoid this excessive energy accumulation in E(k), we ceil off E(k), by discarding any energy over a certain upper limit Ē. Note that, if we utilize y(k) directly in ( 16), passivity would be in general violated, since the switchings of the (discrete) signal y(k) can induce (possibly-destabilizing) energy jumps in the system [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF] and also the frequency/phase contents of y(k) itself may not define a passive mapping with q via K, K f . See [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF], [START_REF] Huang | Implementation and experiments of passive setposition modulation framework for interactive robotic systems[END_REF] for more details on PSPM. We now present the main result of this paper, whose proof we omit here and refer readers to a similar proof in [START_REF] Lee | Haptic teleoperation of multiple unmanned aerial vehicles over the internet[END_REF], that can be easily applied to this Th. 1 while recognizing the BIBO stability of H(s) in ( 14) and non-zero K f in [START_REF] Lee | Passive decomposition of mechanical systems with coordination requirement[END_REF].

Theorem 1 Consider N VPs [START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF] and master device (3) with PSPM-modulated teleoperation control [START_REF] Lee | Passive decomposition of mechanical systems with coordination requirement[END_REF]. ∀T ≥ 0, all the VPs are stable (with bounded ṗi ) with no VP-VP/VP-obstacle collision and VP-VP separations, and ( q, q, qȳ(k)) are also all bounded ∀t ≥ 0.

2) Suppose further that (q, q) → 0, E(k) > 0 ∀k ≥ 0, (x i (t), ẋi (t)) → (p i (t), ṗi (t)) and N s = {1, 2, ..., N }.

Then, (a) if u o i = 0, i = 1, 2, ..., N (i.e., no obstacles), the human user will have collective haptic velocity perception with

q(t) → 1 λN t N i=1 ẋi (19)
f (t) → K f λN t + K λ N s -N t N t N s N i=1 ẋi (20)
and (b) if ẋi = 0, i = 1, 2, ..., N (e.g., stopped by obstacles), human will have collective obstacle haptic perception with

q(t) → - 1 λN t N i=1 u o i (21) f (t) → - K f λN t + K λ N s + N t N t N s N i=1 u o i . (22)
The results of Th. 1 still hold even when the masterslave (discrete) communication is imperfect (e.g., Internet with varying-delay, packet-loss, data-swapping, etc), since: 1) the PSPM can passify any discrete data sequence y(k) regardless of how imperfect its transmission is; and 2) Prop. 1 requires only the boundedness of u t i (t), which is guaranteed by the PSPM's ensuring boundedness of q(t) robustly at the master site (with passive human assumption [START_REF] Lee | Feedback r-passivity of lagrangian systems for mobile robot teleoperation[END_REF] and stability of H ( 14)).

Due to the same reason, with PSPM, we can utilize any arbitrary forms of haptic feedback y(t) other than (15) while preserving master-passivity/slave-stability (e.g., y := ẋl + 1/(λN ) N i=1 u o i , where l represents a leader agent; or even nonlinear function y := y( ẋi , u o i )). Such a flexibility is unusual with typical passivity-enforcing schemes (e.g. wave/PD) and has been exploited in [START_REF] Franchi | Modeling and control of UAV bearing-formations with bilateral high-level steering[END_REF], [START_REF] Son | Human-centered design and evaluation of haptic cueing for teleoperation of multiple mobile robots[END_REF], [START_REF] Lee | Backstepping control of quadrotor-type uavs: trajectory tracking and teleoperation over the internet[END_REF]. We also performed, in [START_REF] Son | Human-centered design and evaluation of haptic cueing for teleoperation of multiple mobile robots[END_REF], a psychophysical study on three different haptic feedback forms for (15) (i.e., y = (1/(λN s )) i∈Ns ẋi , y = (1/(λN s )) i∈Ns u o i , and ( 15)) using novel humanperspective performance measures (i.e., maneuverability and perceptual sensitivity), and revealed their respective effectiveness for specific task objectives. See [START_REF] Son | Human-centered design and evaluation of haptic cueing for teleoperation of multiple mobile robots[END_REF] for more details on this psychophysical evaluation on different haptic feedback for [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF]. See also [START_REF] Son | Haptic feedback can improve human perceptual awareness in multi-robots teleoperation[END_REF] for a psychophysical experiment on the contribution of haptic feedback itself to the user's performance improvement.

The item 2 of Th. 1 illustrates that, on top of haptic feedback, the human user will also have visual feedback, which turns out to be quite useful in some situations (i.e., seeing the device-tip position q, proportional to

N i=1 ẋi (19) or N i=1 u o i (21)
). For instance, if we choose N s = N t and K f = 0 as done in Sec. III-D, although f (t) → 0 with u o i = 0 from (20) (i.e., no force feedback in steady-state), with the visual feedback of q(t) providing the velocity command information and the transient (non-zero) force feedback of f (t) conveying the information of the mismatch between the velocity command and N i=1 ẋi , the human user can still adequately perceive the collective behavior of multiple UAVs and utilize it for teleoperating them. See Sec. III-D. Of course, if desired, by setting K f = 0, we can easily recover non-zero steady-state force feedback even with N t = N s (see [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF]).

Due to the limited master device DOF, there is an unavoidable ambiguity in our haptic feedback, that is, the same pair of (q(t), f (t)) may correspond to the velocity information option is adopted for our experiment in Sec. III-D. We also believe that the master-passivity/slave-stability (with kinematic VPs (10)) of item 1 of Th. 1 is more suitable for our purpose than the usual master-passivity/slavepassivity, since 1) it does not require the human user to, e.g., continuously overcome damping dissipation in the (dynamic) VPs simulation and UAVs dynamics (e.g., wind drag [START_REF] Stramigioli | A novel approach to haptic tele-operation of aerial robot vehicles[END_REF]); and 2) enforcing slave-passivity is not so important here (since VPs are not physically interacting with unknown environment) and rather likely detrimental (i.e., unnecessarily enforce conservative passivity both for the master and slave sides).

III. EXPERIMENTAL TEST-BED AND RESULTS

A. Hardware Setup

The experiments reported in this section were run on a customized version of the MK-Quadro 3 , an open-source platform (see Fig. 3). Four propellers of diameter 0.254[m] are attached to four Roxxy 2827-35 motors, each of which driven via a PWM signal by a BL-Ctrl V2.0 brushless controller. The average power consumption sustainable by the controller is 160[W] and the peak current is 40 [A]. The motors are mounted at the end of 4 aluminum rods joined together in a cross-shape by 2 plastic center-plates. The total span and weight of the frame are 0.5[m] and 0.12[kg], respectively.

We identified the static and dynamical characteristics of the motor/propeller system by means of Nano17 force/torque sensor 4 . As expected, the relation between rotational speed and generated force/torque can be well approximated by a quadratic function. The maximum attainable force and torque are 9.0[N] and 0.141[Nm] respectively. The response from speed command to actual propeller speed was well representable by a first order linear system with a time constant of 0.047 [sec].

The avionics of the UAV is composed of two main parts: an onboard computer and a microcontroller. The 8-bit At-mega1284p microcontroller, clocked at 20MHz, is able to send the desired motor speeds to the brushless controller by means of an I 2 C bus. It also receives data from a 3D LIS344alh accelerometer (0.0039g 0 [m/s 2 ] resolution and ±2g 0 [m/s 2 ] range) and 3 ADXRS610 gyros (0.586[ • /s] resolution and ±300[ • /s] range). The onboard computer is a small Q7 board 5with a Z530 Intel Atom processor, 1 GB DDR2 533 MHz RAM, an 8 GB Flash Disk and a WiFi card, with a total We used a commercial Omega.36 as haptic master device, with 3 fully-actuated translational degrees of freedom. See Fig. 3. The maximum device force is about 10[N] and its workspace is cube-shaped with an edge of 0.12[m]. The device is connected to a computer through USB with 2.5kHz servorate. This computer then communicates to the UAV's onboard compute over an Internet communication.

B. Quadrotor Control and Estimation

The inner/outer-loop controller explained in Sec. II-B is used for each quadrotor to track its own VP. The faster inner-loop attitude control ( 9) is run by the microcontroller at a frequency of 500 Hz while the slower outer-loop position control ( 7)-( 8) is run on the onboard computer at a frequency of about 120 Hz.

The position/orientation data provided by the motion tracking system is directly used by the position controller. However, the update rate of the roll/pitch measurements from the motion tracking system is too slow to be fed back to the (faster) attitude controller. To address this issue, we utilized the standard complementary filters (see, e.g., [START_REF] Mahony | Complementary filter design on the special orthogonal group SO(3)[END_REF]) to produce a high-rate estimate of the roll and pitch angles by data-fusing the gyroscope and accelerometer readings. The dynamics of the employed complementary filters (valid for small angles and accelerations) are given as follows:

φ = wφ + k(ā φ -φ), θ = wθ + k(ā θ -θ)
where w and ā are the gyroscope and accelerometer readings influencing the roll and pitch dynamics, and k is a positive gain. Typical performance of this filter is shown in Fig. 4.

C. Software Setup

The software for our semi-autonomous teleoperation system consists of several processes interconnected though custom interfaces, as depicted in Fig. 5. A C++ algorithmic library provides the signal processing and control methods needed by The microcontroller runs a single process implementing the attitude controller and complementary filter. The onboard computer runs a process implementing the position controller, VP simulation, and another processes for the collective control among the UAVs and communication with the human operator. This onboard computer uses WiFi to communicate via Socket IPC with: 1) the other UAVs' onboard computers; 2) the master haptic device computer; and 3) the external motion tracking system. The haptic device computer runs a local control loop implementing the PSPM algorithm and computing the force cues for the human operator at a frequency of 2.5 kHz.

The Safe Module process, a compact and well-tested custom-built program, is also implemented to mediate the communication between the position and attitude controllers with the aim of taking full control of the UAV in the case of detection of some malfunctioning (e.g., erroneous frequencies, excessive jitter, etc.).

D. Illustrative Experiments

Using the testbed described so far, we conducted experiments to illustrate the theoretical framework presented in this paper. For this, we set N t = N s and K f = 0 as explained after Th. 1 in Sec. II-D (i.e., in steady-state, haptic feedback solely due to obstacles or visual feedback of collective velocity). In the following, we present the results of two representative7 experiments. We also invite readers to watch the attached video where these experiments are shown together with additional materials.

For the first experiment, we design the inter-VP potentials ϕ c ij in [START_REF] Lam | Haptic feedback in uninhabited aerial vehicle teleoperation with time delay[END_REF] to let the four UAVs make a square formation with 2[m] edge in free space (i.e., no obstacles). A narrow passage with 2.5[m] clearance is also installed in the middle of the arena. The human user then tele-pushes the team of four UAVs through this passage with haptic feedback several times to show the overall system behavior. See Fig. 6 for some screenshots from a similar experiment. 3)) provided to the user, with the three lines (red, green and blue) of each plot representing their components in the three orthogonal axes. From there, we can then observe that high force-feedback corresponds to rapid changes in the velocity command (i.e., haptic perception of the velocity mismatch -see Sec. II-D) or to the high values of the obstacle gradient (i.e., haptic obstacle perception: note the opposite signs of u t i and (1/4) 4 i=1 u o i as predicted in ( 22)). Also, note that, in steady-state around 10[sec] with almost-zero obstacle actions, as shown in the item 2-(a) of Th. 1, the average velocity follows the human command with zero control torque (i.e., visual feedback [START_REF] Lee | Passive decomposition and control of nonholonomic mechanical systems[END_REF] with zero haptic feedback [START_REF] Lee | Semi-autonomous teleoperation of multiple wheeled mobile robots over the internet[END_REF]).

The very bottom plot of Fig. 7 contains the evolution of the inter-distances among the UAVs, ||x ix j ||, i, j ∈ {1, 2, 3, 4}. Given the chosen square formation, in free space, the four distances should be 2[m] while the remaining two (diagonals) distances 2 √ 2[m]. These nominal values are plotted with dashed horizontal lines. We can then see there that, due to the presence of obstacles and the teleoperation commands, the actual inter-UAV distances deviate from the nominal ones during the operation, yet, with no collisions/separations among the UAVs. Notice also the correspondence between the phases of large inter-UAV distance errors and high force-feedback in Fig. 7 (i.e., haptic obstacle perception).

Fig. 8 shows the VP-UAV position tracking error ||p ix i ||, i = 1, .., 4, which are fairly small (i.e., less than 5% of the undeformed inter-distance 2[m] among the UAVs in Fig. 7). Similar UAV-VP coordination errors have been observed in all the other trials of the experimental campaign. This small VP-UAV error then implies that our (practical) trajectory tracking controller in Sec. II-B works properly and our assumption of small || ẋi -p i || in Sec. II-B is indeed valide for the experiment (i.e., VP behaviors and UAV behaviors are equivalent). This, along with Fig. 7, also manifests the stable behavior of our multi-UAV teleoperation system.

In Fig. 9, we also present the trajectories of four UAVs during the first 10[sec] of the experimental trial. For better presentation while avoiding unnecessary overlaps, here, we report only the results where the fleet of the UAVs is forced by the human user to pass through the narrow opening only once. From Fig. 9, we can then see that the UAVs' formation shape deforms during the transition and comes back to the undeformed one after traversing the narrow passage. We can also consider this phase together with the first 10[sec] of Fig. 7, where it is clear how the inter-UAV distances are first deviated from the nominal values and then restored (at approximately 9[sec]). Finally notice from Fig. 9 (or Fig. 6) that, due to the rotational symmetry of our controller, the square formation shape of the UAVs rotates in E(3) while interacting with the environment.

To conclude the section, Fig. 10 shows screenshots of the other representative experiment where we designed the VP-VP potential ϕ c ij to generate a tetrahedron formation at rest with a ground obstacle. From the four snapshots in Fig. 10, we can then see that, as the human user tele-drives the UAVs over the obstacle, the whole UAVs' formation rolls over the obstacle, again due to the rotational symmetry of the VP-VP potentials.

IV. SUMMARY AND FUTURE RESEARCH

We proposed a novel haptic teleoperation control framework for multiple UAVs, consisting of three layers: 1) UAV control layer to drive each UAV to follow its own VP; 2) VP control layer to render N VPs as a deformable flying object with inter-VP/VP-obstacle collision avoidance and inter-VP connectivity preservation; and 3) PSPM-based teleoperation layer to allow a human user to tele-control the bulk motion of N VPs with some useful haptic feedback over the Internet. Master-passivity/slave-stability and some asymptotic performance measures are proved. Experiment results are also presented.

Some possible future research directions include: 1) reduction of the number of UAVs directly communicating with the master while retaining the same level of performance (e.g., the same level of controllability [START_REF] Egerstedt | Interacting with networks[END_REF]); 2) elimination of VPs Fig. 10. Screenshot from the second experiment: potentials are designed to make a tetrahedron formation; human user is tasked to guide the UAVs over a ground obstacle. altogether (see [START_REF] Lee | Distributed backstepping control of multiple thrust-propelled vehicles on balanced graph[END_REF] for preliminary results in this direction); 3) application to a real task with the haptic feedback [START_REF] Rodriguez-Seda | Bilateral teleoperation of multiple mobile agents: coordinated motion and collision avoidance[END_REF] perceptually-optimized for that task (by using the method of [START_REF] Son | Human-centered design and evaluation of haptic cueing for teleoperation of multiple mobile robots[END_REF], [START_REF] Son | Haptic feedback can improve human perceptual awareness in multi-robots teleoperation[END_REF]); and 4) experimental comparison with other semiautonomous teleoperation control techniques.

Fig. 1 .

 1 Fig. 1. Semi-autonomous haptic teleoperation with four UAVs and their VPs: gray arrows represent information flow of local autonomous UAV/VP control; blue dashed arrows velocity command for tele-control; and red doted arrows haptic feedback for tele-sensing. Here, the control set is Nt = {1, 3, 4} while the sensing set Ns = {2, 3}.

1 ) 1 ( 17 )

 1117 Closed-loop master system is passive, i.e., ∃d 1 ∈ ℜ s.t., T 0 f T qdt ≥ -d 2∀T ≥ 0. Moreover, if the assumptions of Prop. 1 hold and human user is passive, i.e., ∃ c 2 ∈ ℜ s.t.,

NFig. 3 .

 3 Fig. 3. On the left: our quadrotor UAV setup with its avionics parts. On the right: haptic interface used in the experimental testbed.

Fig. 4 .

 4 Fig. 4. Roll and pitch angles estimate by the complementary filter (red) and their ground truth values (blue) in a typical experiment.

Fig. 5 .

 5 Fig. 5. Software implementation architecture

Fig. 6 .

 6 Fig. 6. Screenshot from the first experiment: potentials are designed to render a square formation; human user is tasked to guide the UAVs into a narrow passage.

Fig. 7 .

 7 Fig. 7. Human velocity command u t i , collective obstacle avoidance gradients 4 i=1 u o i , collective UAVs' velocity 4 i=1 ẋi , control torque of the master device τ , and inter-distances among the UAVs ||x ix j ||.

Fig. 8 .

 8 Fig. 8. VP-UAV position tracking errors for each UAV p i -x i , i = 1 . . . 4.

Fig. 9 .

 9 Fig.9. Trajectories of UAVs projected on the XY and XZ planes during the time interval [0 s,10 s]: the dashed lines and big dots illustrate the formation and locations of the UAVs at 0 s, 5 s, and 10 s, while the black thick lines represent the narrow passage gap.

  p j ||)), not vector-based (i.e. ϕ c ij (p ip j

	9	x 10 4		inter-VPs and VP-obstacle potentials
										inter-VP
										VP-obstacle
	6
	3
	0.1 0	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	Fig. 2. Examples of ϕ c ij and ϕ o ir , to ensure VP-VP/VP-obstacle collision distance > 0.1 and VP-VP separation distance < 1 [15], [34].
	potential ϕ o ij are designed s.t.: 1) they are distance-based (i.e. ϕ c ij (||p i

  ); and 2) ∂ϕ c ij /∂p i and ∂ϕ o ir /∂p i are bounded, if ϕ c ij and ϕ o ir are bounded. Suppose u t i is bounded with ||u t i || ≤ ū ∀t ≥ 0, ∀i ∈ N t ⊂ {1, 2, .., N } and V (0) < M . Suppose further that, if V (t) ≥ M , there exists at least one VP, say the s th VP, s ∈ {1, 2, ..., N }, s.t. Moreover, if u o i ≡ 0 ∀i = 1, 2, ..., N (i.e., no obstacles) and u t i is the same for all the UAVs in N t , we have

	Proposition 1 j∈Ns	∂ϕ c sj ∂p s	+	r∈Os	∂ϕ o sr ∂p s	≥	√	N t + δ st 2	ū	(13)

where N t is the cardinality of N t ; and δ st = 1 if s ∈ N t , and δ st = 0 otherwise. Then, all the N VPs are stable with bounded ṗi ; VP-VP/VP-obstacle collisions are avoided; and VP-VP connectivity is preserved.

  st = δ st , -||W i || 2 + ū||W i || ≤ ū2 /4,and δ 2

	st ū2 . This then implies V (t) ≤ M ∀t ≥ 0, ū2 + proving no inter-VP/VP-obstacle collision and inter-VP sep-i∈Nt\{s} ū2 = N t
	aration. Boundedness of ṗi follows from (10), with bounded
	u c i , u o i and u t i . The last assertion of Prop. 1 can be shown by
	summing up

where we use the facts Preprint -final, definitive version available at http://ieeexplore.ieee.org/ 5 accepted for IEEE T-Mech , Apr. 2013 that δ 2

Even for highly autonomous UAVs, some form of (multisensory) feedback of the remote UAVs' state and their environment would still be beneficial for the human user (e.g., better situational awareness, telepresence).

For brevity of the paper and due to the easiness of inferring how communication unreliability would affect our semi-autonomous teleoperation architecture via PSPM from[START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF],[START_REF] Huang | Implementation and experiments of passive setposition modulation framework for interactive robotic systems[END_REF], in this paper, we omit experimental results with imperfect communication and instead refer readers to[START_REF] Lee | Passive-set-position-modulation framework for interactive robotic systems[END_REF],[START_REF] Huang | Implementation and experiments of passive setposition modulation framework for interactive robotic systems[END_REF]. The obtained theoretical results (e.g., Prop. 1, Th. 1) yet equally hold for the imperfect communication.Preprint -final, definitive version available at http://ieeexplore.ieee.org/ 2 accepted for IEEE T-Mech , Apr. 2013

Preprint -final, definitive version available at http://ieeexplore.ieee.org/

accepted for IEEE T-Mech , Apr. 2013

Preprint -final, definitive version available at http://ieeexplore.ieee.org/ 4 accepted for IEEE T-Mech , Apr. 2013

Preprint -final, definitive version available at http://ieeexplore.ieee.org/ 6 accepted for IEEE T-Mech , Apr. 2013

http://www.mikrokopter.de

http://www.ati-ia.com/

http://www.seco.it/en/, http://www.qseven-standard.org/

http://www.forcedimension.com Preprint -final, definitive version available at http://ieeexplore.ieee.org/

Numerous demonstrations of our teleoperation framework have been performed at the authors' institutions and also at the 2012 International Conference on Intelligent & Autonomous Systems between Korea and Germany, with the system behaving as postulated by the theory (e.g., Th. 1).

Research supported in part by the Korea NRF-MEST 2009-0087640 and R31-2008-000-10008-0, and the Max Planck Society, Germany.