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

Cristian Secchi
email: cristian.secchi@unimore.it

Heinrich H Bülthoff

Antonio Franchi
email: antonio.franchi@laas.fr

A Nonlinear Force Observer for Quadrotors and Application to Physical Interactive Tasks

In order to properly control the physical interactive behavior of a flying vehicle, the information about the forces acting on the robot is very useful. Force/torque sensors can be exploited for measuring such information but their use increases the cost of the equipment, the weight to be carried by the robot and, consequently, it reduces the flying autonomy. Furthermore, a sensor can measure only the force/torque applied to the point it is mounted in. In order to overcome these limitations, in this paper we introduce a Lyapunov based nonlinear observer for estimating the external forces applied to a quadrotor. Furthermore, we show how to exploit the estimated force for shaping the interactive behavior of the quadrotor using Interconnection and Damping Assignment Passivity Based Controller (IDA-PBC). The results of the paper are validated by means of simulations.

I. INTRODUCTION

The physical interaction of flying systems is a challenging control and design problem and it became recently the interest of many researchers. Especially quadrotor UAVs are becoming popular tools for physical interaction tasks [START_REF] Arcas | EU Collaborative Project ICT-287617[END_REF] and flying robots are exploited for accomplishing different goals while interacting with the environment: surface inspection [START_REF] Airobots | EU Collaborative Project ICT-248669[END_REF], tool operations [START_REF] Lee | Mechanics and control of quadrotors for tool operation[END_REF], [START_REF] Nguyen | Hybrid force/motion control and internal dynamics of quadrotors for tool operation[END_REF], desired forces application [START_REF] Gioioso | Turning a near-hovering controlled quadrotor into a 3D force effector[END_REF] or operation of an on board robotic manipulator [START_REF] Lippiello | Exploiting redundancy in cartesian impedance control of UAVs equipped with a robotic arm[END_REF].

One of the main challenge when flying machines are interacting with their environment is the stabilization and the control of this interaction in a meaningful way, i.e., they can exploit interactive forces for achieving a desired task. Passivity based controllers are known to be a very powerful tool for controlling the interaction of a mechanical system with (possibly unknown) environments [START_REF] Secchi | Control of Interactive Robotic Interfaces: a port-Hamiltonian Approach[END_REF] and, in particular, Interconnection and Damping Assignment Passivity Based Controller (IDA-PBC) is an effective technique to control the physical interaction of an underactuated system with the environment, by reshaping its physical properties (e.g., energetic and dissipative structure) [START_REF] Ortega | Interconnection and damping assignment passivity-based control of portcontrolled Hamiltonian systems[END_REF], [START_REF] Wang | Regulation control of underactuated mechanical systems based on a new matching equation of portcontrolled Hamiltonian systems[END_REF]. In [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] it is shown how to use and improve IDA-PBC for controlling the physical interaction of a quadrotor.

Another crucial issue for aerial physical interaction is measuring the external wrenches (i.e., forces and torques) acting on the body of the flying system. This can be done using force/torque sensors that, on one side, give a good and reliable measure but, on the other side, increase the cost of the equipment, the need of power supply and the weight of the robot and, consequently, can decrease the flying time of the aerial vehicle. Another viable solution is the use of a wrench estimator, which represents a cheaper solution and that can give a sufficiently accurate estimate of the wrench if the estimator is properly designed and the necessary measurements (velocity, pose and, if available, acceleration) are accurate enough. Furthermore, a wrench estimator is cheaper, lighter and more power efficient than a force/torque sensor and, therefore, it is a more viable solution for small quadrotors, considering their low load capabilities. Another benefit of using a wrench estimator is that it can estimate an external wrench due to forces applied in any point of the quadrotor, while a sensor measurement is localized to the place of its transducer.

In [START_REF] Augugliaro | Admittance control for physical human-quadrocopter interaction[END_REF] an unscented Kalman Filter for linearized model of a quadrotor is proposed. Observers are a common tool for estimating forces applied on the end-effector (see e.g. [START_REF] Phong | External force estimation using joint torque sensors for a robot manipulator[END_REF]). In [START_REF] Hacksel | Estimation of environment forces and rigid-body velocities using observers[END_REF] nonlinear observers for velocity and force estimation on a rigid body are presented.

The contribution of this paper is a novel observer for estimating the external wrench applied on a quadrotor UAV that is based on the full Lagrangian model of the robot and is, therefore, very general. Furthermore, we exploit the estimated wrench for improving the IDA-PBC controller presented in [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] in a way that the flying system can apply a desired force along any direction.

A summary of the paper follows. In Sec. II we present the nonlinear wrench observer based on Lagrangian dynamics of the quadrotor system. It requires quadrotor position/orientation and their derivatives with applied control inputs to estimate the external wrench acting on the quadrotor Center of Gravity (CoG). Similar to [START_REF] Nikoobin | Lyapunov-based nonlinear disturbance observer for serial n-link robot manipulators[END_REF], we propose a stability analysis using Lyapunov's direct method. In Sec. III we show how to derive a port-Hamiltonian model for quadrotors from its Lagrangian dynamics and how to change its interactive properties by applying IDA-PBC framework. In Section IV we extend the IDA-PBC controller for exerting a desired force at steady state. This is done by finding equilibria conditions and shaping the rotational potential accordingly. Section V presents extensive simulation results, where the capabilities of force observer and force controller are presented. The results show that regardless the presence of noise in state estimation, the force observer gives good approximation, if the external forces are not rapidly varying. We also tested the observer in the non-nominal case of high variations in the external force profile. In this case, as expected, the estimator is less performant than in the nominal case. However, we show that the proposed controller can still stabilize the system and control the physical interaction despite non-perfect output from the observer. This is due to the ability of passivity based control techniques to cope with noise in the estimates and discontinuities of the disturbing forces. Section VI concludes this work with remarks and future works.

II. NONLINEAR WRENCH OBSERVER FOR QUADROTORS

The knowledge of the external force and torque acting on a quadrotor is essential to effectively control the physical interaction with the environment. As we discussed in the introduction it is very convenient in many cases to use a wrench estimator rather than a force/torque sensor mounted, e.g., at the end effector. For this reason, we designed a force/torque observer that is inspired by the nonlinear disturbance observer presented for robotic manipulators in [START_REF] Chen | A nonlinear disturbance observer for robotic manipulators[END_REF] and more general in [START_REF] Nikoobin | Lyapunov-based nonlinear disturbance observer for serial n-link robot manipulators[END_REF].

Consider an inertial frame and a body frame attached to the quadrotor and centered at its center of mass. Both frames follow the NED (North-East-Down) convention. The translational and rotational dynamics of quadrotor can be written using the Lagrangian formalism as follows [START_REF] Carrillo | Quad Rotorcraft Control Vision-Based Hovering and Navigation[END_REF], [START_REF] Raffo | An integral predictive/nonlinear H-infinity control structure for a quadrotor helicopter[END_REF]

mẍ b = -ρR(η)e 3 + mg 0 e 3 + f e J(η)η = -C(η, η) η + τ + τ e (1)
where m > 0 is the mass of the quadrotor; g 0 is the gravitational constant; e 3 = (0 0 1) T ; x b = (x y z) T ∈ R 3 represents the position of the center of mass in the inertial frame; η = (φ θ ψ) T are the roll, pitch, and yaw angles representing the orientation of the body frame w.r.t the inertial frame, and R(η) is the corresponding rotational matrix; ρ ≥ 0 and τ = (τ 1 τ2 τ3 ) T ∈ R 3 are the control thrust amplitude and the control torque (expressed in the body frame), respectively; f e = (f ex f ey f ez ) T ∈ R 3 is the external total force acting on the quadrotor center of mass expressed in the inertial frame; and τ e ∈ R 3 is the external torque expressed in body frame; finally, J(η) is the rotational inertia matrix, and C(η, η) is the matrix representing the Coriolis terms.

Stacking the generalized coordinates, the generalized input and external forces we obtain the following compact expression w e = B(q)q + C(q, q) q + g -G(q)ū

where

B = mI 0 0 J ∈ R 6×6 C = 0 0 0 C ∈ R 6×6 (3)
g = -mg 0 e 3 0 ∈ R 6 G(q) = R(η)e 3 0 0 I ∈ R 6×4
where w e = (f T e τ T e ) T represents the external wrench acting on the quadrotor, q = (x T η T ) T = {q 1 , • • • , q 6 } ∈ R 6 is the quadrotor configuration and ū = (ρ τ T ) T is the control input vector. We propose the following disturbance observer based on [START_REF] Chen | A nonlinear disturbance observer for robotic manipulators[END_REF] and [START_REF] Nikoobin | Lyapunov-based nonlinear disturbance observer for serial n-link robot manipulators[END_REF], ẇe = L(q, q)(w e -ŵe ) = -L(q, q) ŵe + + L(q, q) (B(q)q + C(q, q) q + g -G(q)ū) [START_REF] Nguyen | Hybrid force/motion control and internal dynamics of quadrotors for tool operation[END_REF] where ŵe = ( f T e τ T e ) T ∈ R 6 is the estimated wrench and L(q, q) ∈ R 6×6 will be designed in order to ensure convergence of the observer. Since we do not assume any specific model for the external force, we have no prior information about the derivative of the external forces (or disturbances). Therefore it is assumed that

ẇe = 0, (5)
which is inevitable if one does not know anything about the environment geometry and dynamics. This is an important limitation of this force observer, i.e., it will give poor estimation in case of a rapidly varying external wrench. This can be avoided by designing a good model for the environment, which is however a topic of another study. By defining the observer error

e = w e -ŵe (6)
we can now calculate ė = ẇe -ẇe = L(q, q) ŵe -L(q, q)w e [START_REF] Secchi | Control of Interactive Robotic Interfaces: a port-Hamiltonian Approach[END_REF] which can be expressed by considering [START_REF] Lippiello | Exploiting redundancy in cartesian impedance control of UAVs equipped with a robotic arm[END_REF] as below;

ė + L(q, q)e = 0 (8)
This means that the choice of L(q, q) will directly affect the asymptotic stability of the error dynamics. In order to implement (4) one needs the knowledge of q, q, and q. Measuring or estimating q and q is quite standard in current platforms (see, e.g., [START_REF] Abeywardena | Improved state estimation in quadrotor MAVs: A novel drift-free velocity estimator[END_REF]). However, for many applications, a reliable measurement of the acceleration q (i.e., both the linear and angular acceleration) is not available. For this purpose define the auxiliary vector:

Ψ = ŵe -p( q). (9)
Taking the time derivative we have

ẇe = Ψ + ∂p( q) ∂ q q. ( 10
)
By equating ( 4) and ( 10) we get Ψ + ∂p( q) ∂ q q = -L(q, q)(Ψ + p( q))+ + L(q, q)(B(q)q + C(q, q) q + g -G(q)ū). [START_REF] Augugliaro | Admittance control for physical human-quadrocopter interaction[END_REF] By choosing

∂p( q) ∂ q = L(q, q)B(q) ( 12
)
the dynamics of the nonlinear observer (see Fig. 1) can be written as following Ψ = -L(q, q)Ψ + L(q, q) (C(q, q) q + g -G(q)ū -p( q)) ŵe = Ψ + p( q) (13)

Observer 1 s Ψ Ψ q, q, ū p( q) Ψ ŵe + + Fig. 1. Nonlinear Force Observer.
which is not depending anymore on q. Therefore this scheme (Fig. 1) can be implemented without measuring the acceleration of the generalized coordinates q on commonly available quadrotor platforms.

As it is seen from ( 8), we must choose L(q, q) such a way that the error dynamics become asymptotically stable. Moreover, the decision made in ( 12) brings a strict dependency of L(q, q) on the choice of p( q). Consider the following choice

p( q) = c    q1 . . . q6    (14)
where c > 0 is an observer gain. The choice of p( q) is different from the one made for robot manipulators as shown in [START_REF] Chen | A nonlinear disturbance observer for robotic manipulators[END_REF] and [START_REF] Nikoobin | Lyapunov-based nonlinear disturbance observer for serial n-link robot manipulators[END_REF], since we are dealing with quadrotor dynamics. We obtain that

L(q, q) = cB(q) -1 . (15)
Proposition 1: Consider the wrench estimator [START_REF] Hacksel | Estimation of environment forces and rigid-body velocities using observers[END_REF] and assume that the roll and pitch velocities are bounded, i.e. | φ| < φ and | θ| < θ, where φ, θ ∈ R + . If (5) holds and if L(q, q) is defined as in [START_REF] Chen | A nonlinear disturbance observer for robotic manipulators[END_REF], then ŵe → w e .

Proof: We will prove that the estimation error defined in (6) will asymptotically vanish by showing that the error dynamics ( 8) is asymptotically stable. Let

V (e, q) = e T B(q)e (16)
be a positive definite candidate Lyapunov function. Considering ( 8) and ( 15), we can write:

dV (e, q) dt = 2e T B ė + e T Ḃe = -2e T BLe + e T Ḃe = -2ce T e + e T Ḃe. (17)
The first component of the right hand side of ( 17) is negative definite for c ∈ R + . The second component has an indefinite sign. Nevertheless, since B(q) = B T (q), Ḃ(q, q) is symmetric and, therefore, its eigenvalues are real. From (3) we can compute:

Ḃ(q, q) = 0 0 0 J(η, η) (18)
where the 0 are the zero matrices of proper dimensions. Considering that ( [START_REF] Carrillo | Quad Rotorcraft Control Vision-Based Hovering and Navigation[END_REF], [START_REF] Raffo | An integral predictive/nonlinear H-infinity control structure for a quadrotor helicopter[END_REF]) J = T T JT with

T (η) =   1 0 -s θ 0 c θ s φ c θ 0 -s φ c φ c θ  
where s * = sin( * ), c * = cos( * ), and J = diag{J xx , J yy , J zz } is the inertia matrix in body frame (which is constant w.r.t. q), it is possible to write

J(η, η) = Φ(η) φ + Θ(η) θ (19)
where

Φ =   0 0 0 0 Φ 22 Φ 23 0 Φ 32 Φ 33   Θ =   0 0 Θ 13 0 Θ 22 Θ 23 Θ 31 Θ 32 Θ 33   (20) with Φ 22 = 2J zz c φ s φ Φ 23 = Φ 32 = c θ (-J zz c 2 φ + J yy c φ c θ + J zz s 2 φ ) Φ 33 = 2c 2 θ c φ s φ (J yy -J zz ) Θ 13 = Θ 31 = -J xx c θ Θ 22 = -2J yy c θ s θ Θ 23 = Θ 32 = (J zz c φ -2J yy c θ )s φ s θ Θ 33 = (J xx -J zz c 2 φ -J yy s 2 φ )s θ . (21)
Since | φ| < φ and | θ| < θ, from (21) it is easy to find two finite numbers α, β ∈ R such that α < Ḃij < β, i, j ∈ {1, . . . , 6}. Thus, as shown in [START_REF] Zhan | Extremal eigenvalues of real symmetric matrices with entries in an interval[END_REF], it is always possible to find a finite upper bound λ B for all the possible eigenvalues of Ḃ(q, q):

max q, q λ M Ḃ(q, q) ≤ λ B < ∞ (22)
where λ M Ḃ(q, q) is the maximum eigenvalue of Ḃ(q, q). Thus, we have that

e T Ḃe ≤ λ B e T e (23)
It is therefore possible to choose a c > λ B 2 which implies V is negative definite and that, therefore, e(t) → 0 which proves the statement. Stability (and accuracy) of the observer thus relies on the observer gain, which should be carefully chosen depending on the task.

Remark 1: If measurement for linear acceleration is available, such as quadrotors with Inertial Measurement Unit (IMU), then it is possible to divide the observer in two parts. In first part, the external force acting on the quadrotor can simply be estimated using the first equation given in [START_REF] Arcas | EU Collaborative Project ICT-287617[END_REF]. In second part, the external torque would be estimated the same way as presented above, and the same stability conditions would also hold for it.

III. PORT-HAMILTONIAN MODEL AND IDA-PB CONTROL OF A QUADROTOR

In order to design a controller that is able to stably perform interaction with the environment we decided to use an energetic approach, which is more suited for this purpose. Therefore it is convenient to write the quadrotor equations using the port-Hamiltonian modeling, which provides explicit view on the energetic flow of the system. As shown in [START_REF] Ortega | Interconnection and damping assignment passivity-based control of portcontrolled Hamiltonian systems[END_REF], one can write a port-Hamiltonian system as

   ṡ = J (s) -R(s) ∂H ∂s + G(s)u y = G(s) T ∂H ∂s (24)
where s ∈ R n is the state and H(s) : R n → R represents the total amount of energy stored in the system. Square matrices J (s) = -J (s) T and R(s) ≥ 0 represent the internal energetic interconnections and the dissipation of the port-Hamiltonian system, respectively, G(s) is the input matrix depending on the technology of the system, and the input-output pair (u, y) represents a power port. Any system written in the form of ( 24) is passive w.r.t to this power port, provided that H(s) is lower bounded. In [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] it is shown how to transform Newton-Euler model of a quadrotor into a port-Hamiltonian system with dissipation. A similar strategy can be pursued for transforming the Lagrange dynamics given in (1) into port-Hamiltonian dynamics with dissipation by choosing

τ = C(η, η) η -k d J η + Jτ + (J -I)τ e ( 25
)
where

k d ∈ R + ; τ = (τ 1 τ 2 τ 3 ) T ∈ R 3
and I is the identity matrix of proper dimension. Hence, the rotational dynamics becomes

η = -k d η + τ + τ e (26)
which can be written as

q ṗ = 0 I -I 0 - 0 0 0 R ∂H ∂q ∂H ∂p + 0 0 G I u we (27)
where p = M q ∈ R 6 be momentum variable. R = diag(0 3 k d I 3 ) ∈ R 6×6 is the dissipation in rotational dynamics introduced by (25). Furthermore, u = (ρ τ T ) T ∈ R 4 is the control input vector. The total energy function and the input submatrix G are given by:

H(q, p) = 1 2 p T M -1 p + V (q) = 1 2 p T M -1 p -mgq 3 (28) G = G 1 0 0 I ∈ R 6×4 with G 1 = -R(η)e 3 ∈ R 3 (29)
where

M = mI 0 0 I ∈ R 6×6 (30)
The system in ( 27) is cyclo-passive [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF], which means it is passive as long as the energy function is lower bounded. This perfectly matches with the physical properties of a quadrotor system.

In [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] we propose a controller to reshape the physical properties of the the quadrotor dynamics that has the following form

u = (G T G) -1 G T ∂H ∂q -M M -1 d ∂H d ∂q - -K v G T M -T M T d ∂H d ∂ p + + G + (M M -1 d (I -M d M -1 )w e ) + u h . ( 31
)
Quadrotor System IDA-PBC ū q, q w e Desired Quadrotor System u h Under the action of (31) the original quadrotor behaves as a different one whose dynamics that can be described by the following port-Hamiltonian system

q ṗ = 0 I -I 0 - 0 0 0 R d ∂H d ∂q ∂H d ∂ p + + 0 I we + 0 M d M -1 G u h ( 32
)
where p = M d q is the new momentum, associated to the new inertia matrix:

M d = m d I 0 0 N ( 33
)
where m d ∈ R + and N ∈ R 3×3 is a symmetric positive definite matrix representing the new mass and rotational inertia, respectively; H d is a new energy function

H d = 1 2 pT M -1 d p + V d (q), (34)
where V d is the new potential energy of the form V d (q) = -m d gq 3 + Vd (q 4 , q 5 , q 6 );

R d is the new dissipation matrix; we is the partially compensated external wrench; finally, u h is new control input for controlling the physically modified quadrotor. A schematic representation of IDA-PBC controlled quadrotor system is shown in Fig. 2.

The exhaustive analytical derivation of the controller can be found in [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF]. In order to be implemented, controller (31) requires the measurement of the external wrench w e . The method explained in Sec. II for estimating the external wrench represents thus a valuable tool to allow the actual implementation of the IDA-PBC controller in the real world, as it will be mentioned in Sec. V.

IV. USE OF IDA-PB CONTROL TO EXERT ARBITRARY 3D FORCES

The linear equation in (1) can be rewritten for the reshaped quadrotor as

m d ẍb = -ρR(η)e 3 + mg 0 e 3 =fne +f e ( 36
)
where we denoted with f ne the sum of the total thrust force and the gravity force, i.e., all the forces contributing to the linear motion of the quadrotor except for the external force f e . If ẍb = 0, i.e., the quadrotor moves at constant (possibly zero) velocity, then we have that f ne = -f e , i.e., f ne represents, in this case, the force that the quadrotor is exerting on the external world. In this section we present a method that allows to use the framework presented in Sec. III in order to regulate f ne to a certain desired value. At the end of the section we show how this action can be used, e.g., to counterbalance an external disturbance like a constant wind or to press against a wall or an object with a certain given force.

Denote with f * = (f * x f * y f * z ) T ∈ R 3 the value of the desired f ne . Imposing f ne = f * we obtain the following nonlinear system of equations

ρ   s φ s ψ + c φ c ψ s θ c φ s θ s ψ -c ψ s φ c θ c φ   =   -f * x -f * y -f * z + m d g 0   (37)
after some straightforward algebra we obtain

ρ   c φ s θ -s φ c θ c φ   =   -f * x c ψ -f * y s ψ f * x s ψ -f * y c ψ -f * z + m d g 0   (38)
which, assuming that also ψ is known and denoting is value with ψ ss , it can be solved in the unknown ρ, φ, and θ resulting in

ρ * = f * x 2 + f * y 2 + (f * z -m d g 0 ) 2 (39)
φ * = arcsin -f * x c ψss -f * y s ψss ρ * (40) θ * = -arcsin f * x s ψss -f * y c ψss ρ * c φ * . (41)
Consider now the rotational dynamics of the reshaped quadrotor. We choose Vd (q 4 , q 5 , q 6 ) in (35) as

Vd (q) = 1 2 η T e K p η e ( 42
)
where

K p = diag{k φ p , k θ p , k ψ p }, η e = η -η =   φ - φ θ - θ ψ   . (43)
and φ, θ are parameters to be designed in order to obtain φ → φ * and θ → θ * . By virtue of (42) the rotational dynamics of the reshaped quadrotor can be expressed as

η = -k d η + K p (η -η) + τe (44)
where τe denotes the external torque acting on the reshaped quadrotor (that, we recall, can be estimated using the technique presented in Sec. II). Assuming that τe is constant we obtain the following equilibrium at steady state,

K p (η ss -η) = τe (45)
where η ss = (φ ss , θ ss , ψ ss ) represent the steady state attitude. It is straightforward to see that

ψ ss = τe k ψ p
that can be used in ( 40) and (41) in order to find the exact values of φ * and θ * that are needed to achieve the desired f * . Given those values of φ * and θ * we choose φ and θ such that φ ss = φ * and θ ss = θ * , i.e.,

ψ = ψ * - τex k φ p (46) θ = θ * - τey k θ p . (47)
In summary, by choosing the thrust as in (39) and the potential as in (42), with φ and θ given by ( 46) and ( 47), respectively, we can let f ne converge to f * even in the presence of a disturbing (but constant) external torque τe . This technique can be applied for example to balance an external constant force f e produced by a wind or any other external agent. To this aim one has to estimate f e using the observer presented in Sec. II and then select f * = -f e , which results in a compensation of the external force acting on the quadrotor. In Sec. V we provide a numerical study of this particular case. Another possible application of the technique presented in this paper is the exertion of a constant force to a wall or to load for the purpose, e.g., of pushing, lifting, and so on.

V. SIMULATIONS AND RESULTS

In this section we present different simulation results to show the capabilities of nonlinear force observer presented in Sec. II for different case studies, in which IDA-PBC is used for controlling the physical interaction of the quadrotor system. In order to reproduce realistic scenarios, we added independent Gaussian noises to the measurements of the quadrotor state, based on our laboratory experiences (see Fig. 3). The first simulation aims at showing the accuracy of the force observer in an environment, where disturbing forces and torques are acting on quadrotor CoG. In the second simulation we consider the relevant case where a rigid tool is attached to quadrotor CoG and the tooltip is sliding on the surface of a ceiling (for a detailed explanation see the Sec. V of [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF]). This case study provides a highly varying external force profile, which is a non-nominal situation for using the force observer presented in this work. g Fig. 4. IDA-PBC controlled quadrotor for aerial physical interaction and force estimation using nonlinear force observer. The external wrench is directly acting on the center of gravity of the quadrotor. First 1 N of external force along the x axis is applied between 5 s and 10 s. Then 1 N of external force along the -z axis is applied between 15 s and 20 s. Later a constant torque of 1 Nm is applied around the x axis between 25 s and 30 s. Finally all of these external disturbances are applied between 35 s and 40 s. Third column shows the exact external force and torques (black plots) and their estimations (red plots). The observer gain for force estimation is set to c = 5, and for torque estimation to c = 0.5.

simulation presents a case where quadrotor is exposed to a constant force (e.g., modeling a constant wind force) which has to be balanced while being subject to other disturbances from the environment at the same time. This is the case in which the quadrotor is applying a desired force in order to balance the external force, as explained in Sec. IV. We encourage the reader to watch the video attached to this submission that provides a visual complement to the plots of the three simulation studies.

In simulations we used a quadrotor system with the following parameters. The mass is m = 1 kg, gravity is g 0 = 9.81 m/s 2 and the inertia matrix is J = diag([0.13, 0.13, 0.22]). The legends QCi in the plots represent the quadrotor with the i-th target dynamics assigned using IDA-PBC [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF]. The environment is modeled with no dissipation, which means the aerial drag acting on the body of the quadrotor is not considered.

A. Testing the Observer with a Slowly Varying Force Profile

The aim of the first simulation is to show the accuracy of force estimation, done by nonlinear force observer proposed in Sec. II. As a case study, we choose an external force/torque profile, where first 1 N of force applied along x axis; then along -z axis; later 0.5 Nm of torque around x axis of the quadrotor; and finally all together at the same time. The force profile can be seen in the last row of Fig. 4. The high level control input u h is chosen such as to scale the gravity effect.

As explained in [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF], IDA-PBC can be used to change the interactive behavior of the quadrotors by assigning different dynamic properties. In this case, our goal is to prove that nonlinear force observer gives good external wrench estimation, which can be used for a quadrotor controlled with IDA-PBC. The first two rows of Fig. 4 present the evaluation of quadrotor position, orientation and necessary control inputs when interacting with the external forces. The last row of the same figure shows the estimated forces and torques. It can be seen that force/torque estimation (red plot) is very accurate considering the exact forces (black plot).

B. Testing the Observer with a Rapidly Varying Force Profile

We tested the observer also in an extreme case where there are rapidly varying external forces. The chosen case study is a challenging one, where a rigid tool is sliding on a rough surface. This tool is rigidly attached on a quadrotor, placed above the center of gravity at a position d = (0.2 0 -0.2) T [m], so the surface can be interpreted as a ceiling. This case is shown in [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] without the use of a force observer. The idea is to apply a constant high-level control input, u h , to the system with a desired attitude, and let it fly while using IDA-PBC for stable contact between this uneven ceiling surface and tooltip, and sliding along the surface for performing such as ceiling painting, cleaning or surface inspection tasks. The results are shown in Fig. 5, where t z represents the tooltip position along z axis. QC1 is assigned with smaller desired inertia; and damping along the thrust direction, compared to QC2. QC2 shows a stiffer behavior with respect to the changes on the surface compared to QC1, which has more . Rigid tooltip on a quadrotor interacting with a rough ceiling surface. The blue plot is presenting the behavior of QC1, and the red plot the behavior of QC2. Second row shows the force estimation of the nonlinear force observer. The observer gain is set to c = 20 for both forces and torques. In case of rapidly varying external force profile, the observer gives a less precise estimation, as expected. However, IDA-PBC controller still stabilizes the interaction, thanks to its ability of preserving the passivity of the controlled system.

oscillations and takes more time to reestablish the contact with the surface. These behaviors are consistent with the results in [START_REF] Yüksel | Reshaping the physical properties of a quadrotor through IDA-PBC and its application to aerial physical interaction[END_REF] where the exact knowledge of the external wrench have been used instead of an estimation of it. The main components of the exact and estimated wrench are shown in third column of Fig. 5. As expected the estimator cannot precisely track the rapidly varying external wrench. However, the controller shows the capability to stabilize the physical interaction even when the real and estimated values of the contact force present some discrepancies. The main reason being that our passivity based control method can cope well with these kind of non-idealities. This result makes us optimistic on the success of future experiments.

C. Exerting 3D Desired Force

In Sec. IV we showed how to use the proposed controller and estimator in a way that the quadrotor can exert an arbitrary constant 3D force on the environment. In order to validate the theory in simulation, we consider the case of in which a constant 1 N force is applied to quadrotor along the x axis, continuously. In addition to the force, we apply an impulse of a disturbance force of 1 N along the -z an impulse of a disturbance torque of 1 Nm about the x axis. Finally all these disturbances are applied in the same time. The controller finds a rotational equilibrium, as explained in Sec. IV, so that the quadrotor stops accelerating along the 3D axes by balancing both the external force and torques. The results are presented in Fig. 6. Two different desired dynamics are shown: QC3 is designed with smaller desired inertia and stiffness in rotational dynamics, compared to QC4, which reacts to external effects faster than QC3. This generates a less travelled distance along the direction of motion for QC4. As it is seen in Fig. 6, the quadrotor exerts a counterbalancing force so that the components of the accelerations along all translational axes are converging to zero in steady state (sub figures in first row of Fig. 6). Another sub-figure in the third row shows that f ne,x rapidly converges to f * x = -f ex in order to counterbalance it. Since no dissipation is present in the environment (no aerial drag is acting on the body of the quadrotor), the system floats with a constant velocity in each direction. Similar to the case shown in Fig. 4, the observer provides good estimation of the external wrench.

VI. CONCLUSIONS

We presented a nonlinear force observer for quadrotors to estimate the external wrench acting on the CoG. We tested this observer with IDA-PBC controller, which uses a passivity based control technique used to reshape the physical properties of a quadrotor in a way that it behaves with a given target dynamics. We have presented a force observer that doesn't require additional sensor or hardware, and gives good approximation for slowly varying torques and forces acting to the center of gravity of the system from any direction. We have also shown an improvement of the IDA-PBC controller for exerting desired forces at steady state that can be used e.g., to balance the external disturbances. The results show that the combination of IDA-PBC controller and the nonlinear force observer succesfully stabilizes the system and controls the interaction with the external environment. Even in case of poor force estimation, the interactive behavior remains stable thanks to the robustness of the IDA-PBC.

In the near future we are working toward the realization of an experimental testbed with a real quadrotor that will Force g Fig. 6. IDA-PBC controlled quadrotor is exerting desired forces and counterbalancing external disturbances. The external wrench is directly acting on the center of gravity of the quadrotor. A constant 1 N of external force along the x axis is applied during the whole simulation for replicating a constant wind effect. Then 1 N of external force along the -z axis is applied between 5 s and 10 s. Later a constant torque of 1 Nm is applied around the x axis between 15 s and 20 s. Finally all of these external disturbances are applied between 25 s and 30 s of the simulation. Third row shows the exact external force and torques (black plot); and their estimations (red plot). It can be seen from the first row that the accelerations in all three axes converge to zero at steady state. The sub-figure presented in the last row is showing the desired force applied by the quadrotor to counterbalance the constant wind effect. The observer gain for force estimation is set to c = 5, and for torque estimation to c = 0.5.

hopefully confirm the promising simulative results of the proposed observer/controller scheme.

Fig. 2 .

 2 Fig. 2. Scheme for Controlled Quadrotor System Using IDA-PBC.

Fig. 3 .

 3 Fig.3. Gaussian Noise added to the state measurements. Blue plots on the left column represent the noise in position and orientation (Euler angles); red plots on the right column show the noise in linear velocity and Euler rates. The bias in all measurements are ignored, since they can be eliminated using suitable techniques.

  Finally, lastPreprint version, final version at http://ieeexplore.ieee.org/

Fig. 5

 5 Fig.5. Rigid tooltip on a quadrotor interacting with a rough ceiling surface. The blue plot is presenting the behavior of QC1, and the red plot the behavior of QC2. Second row shows the force estimation of the nonlinear force observer. The observer gain is set to c = 20 for both forces and torques. In case of rapidly varying external force profile, the observer gives a less precise estimation, as expected. However, IDA-PBC controller still stabilizes the interaction, thanks to its ability of preserving the passivity of the controlled system.

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

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

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

Preprint version, final version at http://ieeexplore.ieee.org/ 4 2014 AIM

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