Quadrotor Manipulation System: Development of a Robust Contact Force Estimation and Impedance Control Scheme Based on DOb and FTRLS

TL;DR

This paper presents a robust contact force estimation and impedance control scheme for quadrotor manipulation systems, combining disturbance observer and recursive least squares techniques to improve accuracy and stability in dynamic, nonlinear, and payload-carrying scenarios.

Contribution

It introduces a novel force estimation method based on DOb and FTRLS, addressing limitations of existing techniques and enhancing manipulation accuracy.

Findings

Effective force estimation demonstrated through numerical simulations.

Improved stability and robustness in dynamic and payload variations.

Successful trajectory tracking with impedance control.

Abstract

The research on aerial manipulation systems has been increased rapidly in recent years. These systems are very attractive for a wide range of applications due to their unique features. However, dynamics, control and manipulation tasks of such systems are quite challenging because they are naturally unstable, have very fast dynamics, have strong nonlinearities, are very susceptible to parameters variations due to carrying a payload besides the external disturbances, and have complex inverse kinematics. In addition, the manipulation tasks require estimating (applying) a certain force of (at) the end-effector as well as the accurate positioning of it. Thus, in this article, a robust force estimation and impedance control scheme is proposed to address these issues. The robustness is achieved based on the Disturbance Observer (DOb) technique. Then, a tracking and performance low computational linear controller is used. For teleoperation purpose, the contact force needs to be identified. However, the current developed techniques for force estimation have limitations because they are based on ignoring some dynamics and/or requiring of an indicator of the environment contact. Unlike these techniques, we propose a technique based on linearization capabilities of DOb and a Fast Tracking Recursive Least Squares (FTRLS) algorithm. The complex inverse kinematics problem of such a system is solved by a Jacobin based algorithm. The stability analysis of the proposed scheme is presented. The algorithm is tested to achieve tracking of task space reference trajectories besides the impedance control. The efficiency of the proposed technique is enlightened via numerical simulation.

Tables

Table 1. Table 1: System Parameters

Par.	Value	Unit	Par.	Value	Unit
$m$	$1$	$kg$	$L_2$	$85\times {10}^{-3}$	$m$
$d$	$223\times {10}^{-3}$	$m$	$m_0$	$30\times {10}^{-3}$	$kg$
$I_x$	$13.2\times {10}^{-3}$	$N.m.s^2$	$m_1$	$55\times {10}^{-3}$	$kg$
$I_y$	$12.5\times {10}^{-3}$	$N.m.s^2$	$m_2$	$112\times {10}^{-3}$	$kg$
$I_z$	$23.5\times {10}^{-3}$	$N.m.s^2$	$I_r$	$33.2\times {10}^{-6}$	$N.m.s^2$
$L_0$	$30\times {10}^{-3}$	$m$	$L_1$	$70\times {10}^{-3}$	$m$
$K_{F_1}$	$1.6\times {10}^{-5}$	$kg.m.ra{d}^{-2}$	$K_{F_2}$	$1.2\times {10}^{-5}$	$kg.m.ra{d}^{-2}$
$K_{F_3}$	$1.7\times {10}^{-5}$	$kg.m.ra{d}^{-2}$	$K_{F_4}$	$1.5\times {10}^{-5}$	$kg.m.ra{d}^{-2}$
$K_{M_1}$	$3.9\times {10}^{-7}$	$kg.m^2.ra{d}^{-2}$	$K_{M_2}$	$2.8\times {10}^{-7}$	$kg.m^2.ra{d}^{-2}$
$K_{M_3}$	$4.4\times {10}^{-7}$	$kg.m^2.ra{d}^{-2}$	$K_{M_4}$	$3.1\times {10}^{-7}$	$kg.m^2.ra{d}^{-2}$

Table 2. Table 2: Controller parameters

$Parameter$	$Value$	$Par.$	$Val.$
$M_{n_{\zeta }}$	$diag\{0.02,0.02,2,0.05,0.01,0.01\}$	$A_{e_1}$	$0.16$
$S_{c,d}$	$diag\{20,20,30,50,100,500\}$	$D_c$	$0.01I_6$
$D_{c,d}$	$diag\{15,15,25,100,100,100\}$	$g_{v_i}$	$100$
$M_{n_{\sigma }}$	$diag\{0.05,0.05\}$	$A_{e_2}$	$0.032$
${\eta }_h^{min}$	$0.8$	${\gamma }_g$	$5$
$S_c$	$0.1I_6$	$z_0$	$1$
$K_{p_{\sigma }}$ / $K_{d_{\sigma }}$	$20I_2$	$V_{w_{z_0}}$	$3$

Equations

$$F_j=K_{f_j}{\Omega }_j^2,$$

$$M_j=K_{m_j}{\Omega }_j^2,$$

$$R_b=\left[\begin{array}{ccc}{C}_{\psi }{C}_{\theta }& S_{\varphi }{S}_{\theta }{C}_{\psi }-S_{\psi }{C}_{\varphi }& S_{\psi }{S}_{\varphi }+C_{\psi }{S}_{\theta }{C}_{\varphi }\\ S_{\psi }{C}_{\theta }& C_{\psi }{C}_{\varphi }+S_{\psi }{S}_{\theta }{S}_{\varphi }& S_{\psi }{S}_{\theta }{C}_{\varphi }-C_{\psi }{S}_{\varphi }\\ -S_{\theta }& C_{\theta }{S}_{\varphi }& C_{\theta }{C}_{\varphi }\end{array}\right],$$

$$p_e=p_b+R_b{p}_{eb}^b,$$

$$R_e=R_b{R}_e^b,$$

$${\dot{p}}_e={\dot{p}}_b-Skew(R_b{p}_{eb}^b){\omega }_b+R_b{\dot{p}}_{eb}^b,$$

$${\omega }_e={\omega }_b+R_b{\omega }_{eb}^b,$$

$$v_{eb}^b=J_{eb}^b\dot{\Theta }.$$

$$v_e=J_b{v}_b+J_{eb}\dot{\Theta },$$

$$v_e=J_b{Q}_b{\chi }_b+J_{eb}\dot{\Theta },$$

$$T_b({\Phi }_b)=\left[\begin{array}{ccc}0& -S(\psi )& C(\psi )C(\theta )\\ 0& C(\psi )& S(\psi )C(\theta )\\ 1& 0& -S(\theta )\end{array}\right].$$

$$v_e$$

$$=J_{\zeta }\dot{\zeta }+J_{\sigma }{\dot{\sigma }}_b,$$

$${\dot{\chi }}_e$$

$$=Q_e^{-1}({\Phi }_e)[J_{\zeta }\dot{\zeta }+J_{\sigma }{\dot{\sigma }}_b],$$

$$M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)+{\tau }_w+{\tau }_l=\tau ,\tau =Bu,$$

$$N=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 1& 1& 1& 1& 0& 0\\ {\gamma }_1& -{\gamma }_2& {\gamma }_3& -{\gamma }_4& 0& 0\\ -d& 0& d& 0& 0& 0\\ 0& -d& 0& d& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],$$

$$H=\left[\begin{array}{ccc}{R}_b& O_3& O_2\\ O_3& T_b^T{R}_b& O_2\\ O_{2x3}& O_{2x3}& I_2\end{array}\right].$$

$$\begin{array}{r}{\tau }_l=J^T{F}_e,\\ F_e=S_c{\chi }_e+D_c{\dot{\chi }}_e,\end{array}$$

$$V_{wz}=V_{w_{z_0}}\frac{z}{z_0},$$

$$F_w=0.61\ast A_e{V}_{wz}^2,$$

$$\begin{array}{r}{F}_{wx}=f_{w{x}_1}{z}^{2}sin(\theta )+f_{w{x}_2}{z}^{2}cos(\theta ),\\ F_{wy}=f_{w{y}_1}{z}^{2}sin(\varphi )+f_{w{y}_2}{z}^{2}cos(\varphi ),\end{array}$$

$$\begin{array}{r}{\tau }^{dis}=(M(q)-M_n)\ddot{q}+{\tau }^d,\\ {\tau }^d=C(q,\dot{q})\dot{q}+G(q)+d_{ex}.\end{array}$$

$$M_n\ddot{q}+{\tau }^{dis}=\tau .$$

$$\begin{array}{r}\tau =\frac{1}{(1-Q(s))}[M_n{\ddot{q}}^{des}-Q(s)M_n\ddot{q}],\\ =M_n{\ddot{q}}^{des}+M_{n}P{e}_v,e_v={\dot{q}}^{des}-\dot{q}.\end{array}$$

$$\begin{array}{r}M(q){\dot{e}}_v+C(q,\dot{q})e_v+K_v{e}_v=\delta ,\\ K_v=P{M}_n,\end{array}$$

$$\begin{array}{r}\delta =\Delta M(q){\ddot{q}}^{des}+C(q,\dot{q}){\dot{q}}^{des}+G(q)+d_{ex},\\ \Delta M(q)=M(q)-M_n.\end{array}$$

$$V=\frac{1}{2}{e}_v^{T}M(q)e_v.$$

$$\dot{V}=e_v^{T}M(q){\dot{e}}_v+\frac{1}{2}{e}_v^T\dot{M}(q)e_v.$$

$$\dot{V}=e_v^T\delta -e_v^T{K}_v{e}_v+\frac{1}{2}{e}_v^T(\dot{M}(q)-2C(q,\dot{q}))e_v.$$

$${\lambda }_{min}{\Vert \nu \Vert }^2\le {\nu }^{T}M(q)\nu \le {\lambda }_{max}{\Vert \nu \Vert }^2,$$

Full text

1 Introduction

Recently, Unmanned Aerial Vehicles (UAVs) especially multi-rotors type, receive great attention due to their higher degree of mobility, speed and capability to access to regions that are inaccessible to ground vehicles. However, UAV as a standalone vehicle has a limited functionality to the search and surveillance applications.

Due to their superior mobility, much interest is given to utilize them for aerial manipulation and thus the application of UAV manipulation systems have been expanded dramatically. Applications of such systems include inspection, maintenance, structure assembly,firefighting, rescue operation, surveillance, or transportation in locations that are inaccessible, very dangerous or costly to be accessed from the ground.

Research on quadrotor-based aerial manipulation can be divided into different approaches based on the tool attached to the UAV including gripper based [1], cables based [2, 3], multi-DoF robotic manipulator based[4, 5], multi-DoF dual-arms manipulator based [6], compliant manipulator -based [7], Hybrid rigid/elastic-joint manipulator [8].

In the gripper/ tool-based approach, the attitude of the payload/tool is restricted to that of the quadrotor, and hence, the resulting aerial system has independent 4 DOFs; three translational DOFs and one rotational DOF (Yaw), i.e., the gripper/tool cannot posses pitch or roll rotation without moving horizontally. The second approach is to suspend a payload with cables but this approach has a drawback that the movement of the payload cannot be always regulated directly. To cope up with these limitations, another approach is developed in which a quadrotor is equipped with a robotic manipulator that can actively interact with the environment. Very few reports exist in the literature that investigate the combination of aerial vehicle with robotic manipulator. Kinematic and dynamic models of the quadrotor combined with arbitrary multi-DOF robot arm are derived using the Euler-Lagrangian formalism in [9]. In [10], a quadrotor with light-weight manipulators, three 2-DOF arms, are tested. In [4], an aerial manipulation using a quadrotor with a 2-DOF robotic arm is presented but with certain topology that disable the system from making arbitrary position and orientation of the end-effector. In this system, the axes of the manipulator joints are parallel to each other and parallel to one in-plane axis of the quadrotor. Thus, the system cannot achieve orientation around the second in-plane axis of the quadrotor without moving horizontally.

From the above discussion, the current introduced systems in the literature that use a gripper suffers from the limited allowable DOFs of the end-effector. The other systems have a manipulator with either two DOFs but in certain topology that disables the end-effector to track arbitrary 6-DOF trajectory, or more than two DOFs which decreases greatly the possible payload carried by the system.

In [11, 12, 5], the authors propose a new aerial manipulation system that consists of 2-link manipulator, with two revolute joints whose axes are perpendicular to each other and the axis of the first joint is parallel to one in-plane axis of the quadrotor. Thus, the end-effector is able to reach arbitrary position and orientation without moving horizontally with minimum possible actuators.

In order to achieve position holding during manipulation, uncertainties and disturbances in the system such as wind, contact forces, measurement noise have to be compensated by using a robust control scheme. Disturbance Observer (DOb)-based controller is used to achieve a robust motion control [13, 14]. The DOb estimates the nonlinear terms and uncertainties then compensates them such that the robotic system acts like a multi-SISO linear systems. Therefore, it is possible to rely on a standard linear controller to design the controller of the outer loop such that the system performance can be adjusted to achieve desired tracking accuracy and speed. In [15, 16, 17], DOb-based motion control technique is applied to robotic-based systems and gives efficient results.

In the motion control of the aerial manipulator, achievement of the compliance control is very important because the compliance motion makes possible to perform flexible motion of the manipulator according to desired impedance [18]. This is very critical demand in applications such as demining and maintenance. In the compliance control, end-effector position and generated force of the manipulator are controlled according to the reaction force detected by the force sensor. In this method, the desired impedance is selected arbitrary in the controller. However, the force sensor is essential to detect the reaction force as presented in [19, 20, 21]. On line identified environment impedance has also been used for transparency in teleoperation systems [22]. These problems are more severe when environment displays sudden changes in its dynamic parameters which cannot be tracked by the identification process. In [23], it is found that in order to faithfully convey to the operator the sense of high frequency chattering of contact between the slave and hard objects, faster identification and structurally modified methods were required. However, these methods need the measurement of force.

Several techniques are proposed to estimate the contact force and the environment dynamics. In [24], the DOb and Recursive Least Squares (RLS) are used to estimate the environment dynamics. However, in this method, two DObs are used besides the RLS, and the estimation of contact force is activated only during the instance of contacting, thus there is a need to detect the instant at which the contact occurs. However, this is not practical approach especially if we target autonomous system. In [25, 26, 27, 28, 29, 30], several techniques are proposed to achieve force control without measuring the force. However, these techniques are based on ignoring some dynamics and external disturbances which will produce inaccurate force estimation. In [31, 32], an impedance control is designed for aerial manipulator without the need to measure/estimate the contact force. However, in such work, the authors neglect some dynamics as well as external disturbances, in addition to, the proposed algorithm is model-based and it does not have a robustness capability. In [33], a scheme is proposed which allows a quadrotor to perform tracking tasks without a precise knowledge of its dynamics and under the effect of external disturbances and unmodeled aerodynamics. In addition, this scheme can estimate the external generalized forces. However, as the authors claim, this estimator can work perfectly with constant external disturbances. In addition, the estimated forces contain many different types of forces such as wind, payload, environment impacts, and unmodeled dynamics. Thus, it can not isolate the end-effector force only from the others. The authors in [34] present a model-based method to estimate the external wrench of a flying robot. However, this method assumes that there are no modeling errors and no external disturbance. Moreover, it estimates the external force as one unite and it can not distinguish between external disturbance and the end-effector force which we need to calculate for teleoperation purposes. In addition, it uses a model based control which needs a full knowledge of the model.

In this article, a new scheme is proposed to cope up with these limitations of the currently developed techniques to solve the issues of this complicated multibody robotic system. Firstly, a DOb inner loop is used to estimate both the system nonlinearities and all external forces to compensate for them, as a result, the system acts like a linear decoupled MIMO system. Secondly, a fast tracking RLS algorithm is utilized with the linearization capabilities of DOb to estimate the contact force, in addition to, it enables the user to sense the contact force at the end-effector that it is not available in the current developed schemes. Thirdly, a model-free robust impedance control of the quadrotor manipulation system is implemented. The DOb is designed in the quadrotor/joint space while the impedance control is designed in the task space such that the end-effector can track the desired task space trajectories besides applying a specified environment impedance. Thus, Fourthly, a Jacobian based algorithm is proposed to transform the control signal from the task space to the quadrotor/joint space coordinates. The rigorous stability analysis of the proposed scheme is presented. Finally, the system model is simulated in MATLAB taking in to considerations all the non-idealities and based on real parameters to emulate a real system.

2 System Modeling

Fig. 1 presents a 3D CAD model of the proposed quadrotor-based aerial manipulator. The system is composed of a manipulator mounted on the bottom center of a quadrotor.

System geometrical frames, which are assumed to satisfy the Denavit-Hartenberg (DH) convention, are illustrated in Fig. 2. The manipulator has two revolute joints. The axis of the first revolute joint, $z_{0}$, is parallel to the quadrotor $x$-axis. The axis of the second joint, $z_{1}$, is normal to that of the first joint and hence it is parallel to the quadrotor $y$-axis at the extended configuration. Therefor, the pitching and rolling rotation of the end-effector is allowable independently from the horizontal motion of the quadrotor. Hence, with this proposed aerial manipulator, it is possible to manipulate objects with arbitrary location and orientation. Consequently, the end-effector can make motion in 6-DOF with minimum possible number of actuators/links that is critical factor in flight.

The quadrotor components are designed to achieve a payload capacity of $500$ g. Asctec pelican quadrotor [35] is utilized as a quadrotor platform. The maximum thrust force for each rotor is $6$N. The arm is designed so that the total weight of the arm is $200$ g, it has a maximum reach in the range of $22$, and it can carry a payload of $200$ g. It has three DC motors, (HS-422 (Max torque = $0.4$ N.m) for gripper, HS-5485HB (Max torque = $0.7$ N.m) for joint $1$, and HS-422 (Max torque = $0.4$ N.m) for joint $2$).

The angular velocity of each rotor $j$ is $\Omega_{j}$ and it generates thrust force $F_{j}$ and drag moment $M_{j}$ that are given by

[TABLE]

[TABLE]

where $K_{f_{j}}$ and $K_{m_{j}}$ are the thrust and drag coefficients.

2.1 Kinematics

Let $\Sigma_{b}$, $O_{b}$- $x_{b}$ $y_{b}$ $z_{b}$, represents the quadrotor body-fixed reference frame with origin at the quadrotor center of mass, see Fig. 2. Its position with respect to the world-fixed inertial reference frame, $\Sigma$, $O$- $x$ $y$ $z$, is given by the $(3\times 1)$ vector $p_{b}=[x,y,z]^{T}$, while its orientation is given by the rotation matrix $R_{b}$ which is given by

[TABLE]

where $\Phi_{b}$=$[\psi,\theta,\phi]^{T}$ is the triple $ZYX$ yaw-pitch-roll angles. Note that $C(.)$ and $S(.)$ are short notations for $cos(.)$ and $sin(.)$. Let us consider the frame $\Sigma_{e}$, $O_{2}$- $x_{2}$ $y_{2}$ $z_{2}$, attached to the end-effector of the manipulator, see Fig. 2. Thus, the position of $\Sigma_{e}$ with respect to $\Sigma$ is given by

[TABLE]

where the vector $p^{b}_{eb}$ describes the position of $\Sigma_{e}$ with respect to $\Sigma_{b}$ expressed in $\Sigma_{b}$. The orientation of $\Sigma_{e}$ can be defined by the rotation matrix

[TABLE]

where $R^{b}_{e}$ describes the orientation of $\Sigma_{e}$ w.r.t $\Sigma_{b}$. The linear velocity $\dot{p}_{e}$ of $\Sigma_{e}$ in the world-fixed frame is obtained by the differentiation of (4) as

[TABLE]

where $Skew(.)$ is the $(3\times 3)$ skew-symmetric matrix operator [36], while $\omega_{b}$ is the angular velocity of the quadrotor expressed in $\Sigma$. The angular velocity $\omega_{e}$ of $\Sigma_{e}$ is expressed as

[TABLE]

where $\omega^{b}_{eb}$ is the angular velocity of the end-effector relative to $\Sigma_{b}$ and is expressed in $\Sigma_{b}$.

Let $\Theta=[\theta_{1},\theta_{2}]^{T}$ be the $(2\times 1)$ vector of joint angles of the manipulator. The $(6\times 1)$ vector of the generalized velocity of the end-effector with respect to $\Sigma_{b}$, $v^{b}_{eb}=[\dot{p}^{bT}_{eb},\omega^{bT}_{eb}]^{T}$, can be expressed in terms of the joint velocities $\dot{\Theta}$ via the manipulator Jacobian $J^{b}_{eb}$ [37], such that

[TABLE]

From (6) and (7), the generalized end-effector velocity, $v_{e}=[\dot{p}^{T}_{e},\omega^{T}_{e}]^{T}$, can be expressed as

[TABLE]

where $v_{b}=[\dot{p}^{T}_{b},\omega^{T}_{b}]^{T}$, $J_{b}=\begin{bmatrix}I_{3}&-Skew(R_{b}p^{b}_{eb})\\ O_{3}&I_{3}\end{bmatrix},$ $J_{eb}=\begin{bmatrix}R_{b}&O_{3}\\ O_{3}&R_{b}\end{bmatrix}J^{b}_{eb}$,

where $I_{m}$ and $O_{m}$ denote $(m\times m)$ identity and $(m\times m)$ null matrices, respectively. If the attitude of the vehicle is expressed in terms of yaw-pitch-roll angles, then (9) will be

[TABLE]

with $\chi_{b}=\begin{bmatrix}p_{b}\\ \Phi_{b}\end{bmatrix},$ $Q_{b}=\begin{bmatrix}I_{3}&O_{3}\\ O_{3}&T_{b}\end{bmatrix},$ where $T_{b}$ describes the transformation matrix between the angular velocity $\omega_{b}$ and the time derivative of Euler angles $\dot{\Phi}_{b}$, and it is given as

[TABLE]

Since the vehicle is an under-actuated system, i.e., only $4$ independent control inputs are available for the 6-DOF system, the position and the yaw angle are usually the controlled variables. The pitch and roll angles are used as intermediate control inputs to control the horizontal position. Hence, it is worth rewriting the vector $\chi_{b}$ as follows $\chi_{b}=\begin{bmatrix}\eta_{b}\\ \sigma_{b}\end{bmatrix},$ $\eta_{b}=\begin{bmatrix}p_{b}\\ \psi\end{bmatrix},$ $\sigma_{b}=\begin{bmatrix}\theta\\ \phi\end{bmatrix}.$

Thus, the differential kinematics (10) will be

[TABLE]

where $\zeta=[\eta_{b}^{T},\Theta^{T}]^{T}$ is the vector of the controlled variables, $J_{\eta}$ is composed by the first 4 columns of $J_{b}Q_{b}$, $J_{\sigma}$ is composed by the last 2 columns of $J_{b}Q_{b}$, and $J_{\zeta}=[J_{\eta},J_{eb}]$.

If the end-effector orientation is expressed via a triple of Euler angles, $ZYX$, $\Phi_{e}$, the differential kinematics (12) can be rewritten in terms of the vector $\dot{\chi}_{e}=[\dot{p}^{T}_{e},\dot{\Phi}_{e}^{T}]^{T}$ as follows

[TABLE]

where $Q_{e}$ is the same as $Q_{b}$ but it is a function of $\Phi_{e}$ instead of $\Phi_{b}$.

2.2 Dynamics

The equations of motion of the proposed robot have been derived in details in [11]. The dynamical model of the quadrotor-manipulator system can be reformulated in a matrix form as

[TABLE]

where $q=[x,y,z,\psi,\theta,\phi,\theta_{1},\theta_{2}]^{T}$ $\in R^{8}$ represents the vector of the generalized coordinates, $M$ $\in R^{8\times 8}$ denotes the symmetric and positive definite inertia matrix of the system, $C$ $\in R^{8\times 8}$ represents the Coriolis and centrifugal terms, $G$ $\in R^{8}$ represents the gravity term, $\tau_{w}$ $\in R^{8}$ is vector of the external disturbances, $\tau_{l}$ $\in R^{8}$ is vector of the contact force effect, $\tau$ $\in R^{8}$ is the generalized input torques/forces, $u=[F_{1},F_{2},F_{3},F_{4},\tau_{m_{1}},\tau_{m_{2}}]^{T}$ is vector of the actuator inputs, and $B=HN$ is the input matrix which is used to produced the body forces and moments from the actuator inputs. The control matrix, $N$, is given as

[TABLE]

where $\gamma_{j}=K_{m_{j}}/K_{f_{j}}$, and $H$ $\in R^{8\times 8}$ is matrix that transforms body input forces to be expressed in $\Sigma$ and is given by

[TABLE]

The environment dynamics, contact force, $\tau_{l}$, can be modeled as following:

[TABLE]

where $S_{c}=diag\{S_{c_{1}},S_{c_{2}},S_{c_{3}},S_{c_{4}},S_{c_{5}},S_{c_{6}}\}$ and $D_{c}=diag\{D_{c_{1}},D_{c_{2}},D_{c_{3}},D_{c_{4}},D_{c_{5}},D_{c_{6}}\}$ represent the environment stiffness and the environment damping, receptively.

The wind dynamics, $\tau_{w}$, can be modeled as following [38, 39, 40]:

The average wind velocity is determined by

[TABLE]

where $V_{wz}$ is the wind velocity at altitude $z$, $V_{w_{z_{0}}}$ is the specified (measured) wind velocity at altitude $z_{0}$. To simulate wind disturbances, it is worth calculating the wind force, $F_{w}$, which influences the platform than the wind velocity. This force can be determined by

[TABLE]

where $0.61$ is used to convert wind velocity to pressure, and $A_{e}$ is the influence effective area which depends on the quadrotor structure and its orientation.

This force can be projected on the axes of frame $\Sigma$ as

[TABLE]

where $f_{wx_{1}}=0.61*A_{e_{1}}(\frac{V_{w_{z_{0}}}}{z_{0}})^{2}cos(\psi_{w})$, $f_{wx_{2}}=0.61*A_{e_{2}}(\frac{V_{w_{z_{0}}}}{z_{0}})^{2}cos(\psi_{w})$, $f_{wy_{1}}=0.61*A_{e_{1}}(\frac{V_{w_{z_{0}}}}{z_{0}})^{2}sin(\psi_{w})$, $f_{wy_{2}}=0.61*A_{e_{2}}(\frac{V_{w_{z_{0}}}}{z_{0}})^{2}sin(\psi_{w})$, $\psi_{w}$ represents the angle of wind direction, and both $A_{e_{1}}$ and $A_{e_{2}}$ depend on the quadrotor design parameters.

3 Controller Design

3.1 Control Objectives

Our goal is to design of estimation and control system to achieve the following objectives:

Robust Stability: The robotic system in Fig. 3 is robust and stable against the external disturbances, parameters uncertainties, and noises. 2. 2.

Force Estimation: The end-effector contact force has to be estimated with fast response and the estimation error tends to zero as the time tends to $\infty$. 3. 3.

$6$-DOF Impedance Control: In the presence of the applied force/desired impedance at the end-effector, the end-effector tracking error tends to zero as time tends to $\infty$.

To this end, we propose a control scheme as shown in Fig. 3. In this control strategy, the system nonlinearities, external disturbances (wind), $\tau_{w}$, and contact force, $\tau_{l}$, are treated as disturbances, $\tau^{dis}$, that will be estimated, $\hat{\tau}^{dis}$, and compensated by the DOb in the inner loop. The system can be now tackled as linear SISO plants. The output of DOb with system measurements of both joint and task spaces variables are used as the inputs to the FTRLS to obtain the end-effector contact force $\hat{F}_{e}$. The task space impedance control is used in the external loop of DOb and its output is transformed to the joint space through a transformation algorithm.

3.2 Disturbance Observer Loop

A block diagram of the DOb inner loop is shown in Fig. 4. In this figure, $M_{n}$ $\in R^{8\times 8}$ is the system nominal inertia matrix, $\tau$ and $\tau^{des}$ are the robot and desired inputs, respectively, $P=diag([g_{1},...,g_{i},...,g_{8}])$ with $g_{i}$ is the bandwidth of the $i^{th}$ variable of $q$, $Q(s)=diag([\frac{g_{1}}{s+g_{1}},...,\frac{g_{i}}{s+g_{i}},...,\frac{g_{8}}{s+g_{8}}])$ $\in R^{8\times 8}$ is the matrix of the low pass filter of DOb. The DOb requires velocity measurement. Practically, the velocity have to be fed through a low pass filter, $Q_{v}(s)=diag([\frac{g_{v_{1}}}{s+g_{v_{1}}},...,\frac{g_{v_{i}}}{s+g_{v_{i}}},...,\frac{g_{v_{8}}}{s+g_{v_{8}}}])$ $\in R^{8\times 8}$, and with cut-off frequency of $P_{v}=diag([g_{v_{1}},...,g_{v_{i}},...,g_{v_{8}}])$. $\tau^{dis}$ represents the system disturbances, and $\hat{\tau}^{dis}$ is the estimated disturbances.

If we apply the concept of disturbance observer to the proposed system, the independent coordinate control is possible without considering coupling effect of other coordinates. The coupling terms such as centripetal and Coriolis and gravity terms are considered as disturbance and compensated by feed forward the estimated disturbance torque.

The disturbance $\tau^{dis}$ can be assumed as

[TABLE]

Substituting from (21), then (14) can be rewritten as

[TABLE]

The control input, $\tau$, see Fig. 4, is given as

[TABLE]

Applying this control law results in

[TABLE]

where

[TABLE]

Stability of this inner loop can be proved as following:

To simplify the analysis, let us ignore the effect of the velocity filter which will be considered later.

Let us use a Lyapunov function as

[TABLE]

The time derivative of this function is

[TABLE]

Substituting from (24), then (27) becomes

[TABLE]

To complete this proof, the properties of the dynamic equation of motion (14) will be utilized. Theses properties are [41, 36]:

Property 3.1

[TABLE]

Property 3.2

[TABLE]

where $\nu\in R^{8}$ represents a $8$-dimensional vector, and $\lambda_{min}$ and $\lambda_{max}$ are positive real constants.

Substituting from (30), then (28) will be

[TABLE]

From property (29), one can get

[TABLE]

From the analysis presented in [42], (32) can be reformulated as

[TABLE]

Thus, the error dynamics is $L_{p}$ input/output stable with respect to the pair ($\delta$,$e_{v}$) for all $p\in[1,\infty]$ with the assumption that the system states, $q$ and $\dot{q}$, are bounded.

If one considers the effect of using the velocity filter, then the characteristic equation of the inner loop is

[TABLE]

where $\alpha_{i}=\frac{M_{n_{ii}}}{M_{ii}}$.

To improve the robustness, the damping coefficient of this equation, which is $0.5\sqrt{\frac{g_{i}g_{v_{i}}}{\alpha g_{i}}}$, should larger than or equal $0.707$ . Therefore, the following inequality

[TABLE]

should be hold. Recasting (35) with respect to $K_{v}$ gives to

[TABLE]

Summarizing, (32) shows that the stability and robustness of the control system is enhanced by increasing $K_{v}$, i.e., by increasing $M_{n}$ and $P$, but without violating the robustness constraint given in (36).

If the DOb performs well, that is $\hat{\tau}^{dis}$ = $\tau^{dis}$, the dynamics from the DOb loop input $\tau^{des}$ to the output of the system is given as

[TABLE]

Since $M_{n}$ is assumed to be a diagonal matrix, the system can be considered as a decoupled linear multi SISO systems as

[TABLE]

or in the acceleration space as:

[TABLE]

The next step is to design an Impedance tracking based controller in the outer loop for the system of (39).

3.3 Fast Tracking Recursive Least Squares

In this part, we develop a technique which utilizes a Fast Tracking Recursive Least Squares (FTRLS) to estimate the contact force with the aid of the DOb linearization capabilities. The FTRLS algorithm is one of the fast online least squares-based identification methods used for the identification of environments with varying dynamic parameters [43, 44]. To apply FTRLS, the dynamic equations (14- 20) have to be parametrized (i.e., to be product of measurement data regressor and dynamic parameters) as follows:

The system dynamic part, $\tau_{int}=M(q)\ddot{q}+C(q,\dot{q})\dot{q}+G(q)$, can be rewritten as the product of data regressor, $Y_{i}(q,\dot{q},\ddot{q})$, and platform parameters, $h_{i}$. The environment dynamics, $\tau_{l}$, can be reformulated as $Y_{l}(q,\dot{q},\ddot{q},\chi_{e},\dot{\chi}_{e})*h_{l}$, where, $Y_{l}=J^{T}Y_{e}$, $Y_{e}$ is a function of the end effector states, ($\chi_{e}$,$\dot{\chi}_{e}$), and $h_{l}$ is the environment parameters $S_{c}$ and $D_{c}$. Finally, the wind effect is formulated as $Y_{w}(z,\theta,\phi)*h_{w}$, where $h_{w}$ is the wind parameters. Thus, the total dynamics can be reformulated as

[TABLE]

where $Y$ $\in R^{8\times 40}$ and $h$ $\in R^{40}$ are the data regressor and parameters vector of (14), respectively.

The parameter estimation error is

[TABLE]

while the estimation error is

[TABLE]

By minimizing a cost function with respect to the parameter estimation error, one can find the time derivative of the estimated parameters vector, $\hat{h}$, as following

[TABLE]

where $R(t)$ is the parameters’ covariance matrix, and it can be calculated from

[TABLE]

where $\eta_{h}$ is the forgetting factor, and it is given as

[TABLE]

where $\eta_{h}^{min}$ is a constant representing the minimum forgetting factor, $NINT(.)$ is the round-off operator, and $\gamma_{g}$ is a design constant. This adaptive formulation of the forgetting factor enables the RLS to track the non-stationary parameters to be estimated.

The convergence/stability ($\tilde{h}(t)\longrightarrow 0$) proof of this algorithm can be implemented as following:

Let us assume the Lyapunov function as

[TABLE]

If $R^{-1}(t)$ is chosen to be positive definite, then $V(t)$ will be positive definite. To prove the positive definiteness of $R^{-1}(t)$, let us use the solution of the differential equation (44) which is

[TABLE]

where $\Phi_{h}^{T}(t,t_{0})$ is the state transition matrix of a system described by $\dot{\upsilon}(t)=-\frac{1}{2}\eta_{h}\upsilon(t)$. Thus, by choosing $R^{-1}(t_{0})>0$, then the first term in (47) will be positive definite. The second term is also positive definite. As a result, the proposed covariance matrix update formula is positive definite, and thus, the chosen Lyapunov function (46) is positive definite.

The time derivative of Lyapunov function is

[TABLE]

However, by differentiating both sides of (41) with respect to time, one can find that $\dot{\tilde{h}}=-\dot{\hat{h}}$, by substituting from the proposed formula of $\dot{\hat{h}}$ (43) and (42), then

[TABLE]

Substituting from (49) in (48), then $\dot{V}(t)$ will be

[TABLE]

Substituting from the proposed formula (44) for $\dot{R^{-1}}$ into (50), then

[TABLE]

Thus, the time derivative of $V(t)$ is negative definite which ensures the asymptotic stability of the estimation error ($\tilde{h}(t)\longrightarrow 0$ as $t\longrightarrow\infty$)

Finally, for both teleoperation impedance control purposes, the user can calculate the estimated environment impedance, contact force, from

[TABLE]

Therefore, unlike the current the developed schemes, with this technique, one can isolate and estimate the end-effector contact force apart from the whole estimated forces in the systems.

3.4 Impedance Control

The objective of the impedance control is to regulate the end-effector interaction force, which may vary due to the uncertainty in the location of the interaction point and/or the structural properties of the environment, besides achieving task space trajectory tracking. The linear impedance control is designed in the task space. This is based on the linearization effect of the designed DOb in the joint space. The desired acceleration in the task space, $\ddot{\chi}^{des}_{e}$, can be calculated from

[TABLE]

where $S_{c,d}$ and $D_{c,d}$ are the desired values of $S_{c}$ and $D_{c}$ respectively, which determine the desired impedance that the end-effector will apply to the environment.

Let us define the quadrotor/joint space tracking error as

[TABLE]

while the task space tracking error can be defined as

[TABLE]

where $\chi_{e,r}$, $\dot{\chi}_{e,r}$, and $\ddot{\chi}_{e,r}$ are the reference trajectories for the position, velocity, acceleration in the task space, respectively which are chosen to be bounded and continuous. $q_{r}$, $\dot{q}_{r}$, and $\ddot{q}_{r}$ are the reference trajectories for the position, velocity, acceleration in the quadrotor/joint space, respectively. Transformation from the task space to quadrotor/joint space will be implemented via the inverse of system Jacobian.

The relation between the inner loop and the outer loop errors can be obtained as follows.

The DOb loop error can be expressed in the task space, $e_{v,e}$, via the Jacobian by

[TABLE]

where $e_{v,e}=\dot{\chi}_{e}^{des}-\dot{\chi}_{e}$. From the previous analysis, it is proved that $e_{v}$ is bounded as in (33). If we define $\dot{e}_{v,e}=\ddot{\chi}_{e}^{des}-\ddot{\chi}_{e}$, then by substituting from (53), one can get

[TABLE]

Equation (57) can be reformulated in a state space form as

[TABLE]

where $X_{e}=\begin{bmatrix}e_{e}\\ \dot{e}_{e}\\ \end{bmatrix}$, $A_{e}=\begin{bmatrix}O_{6}&I_{6}\\ -S_{c,d}&-D_{c,d}\end{bmatrix}$, $B_{e}=\begin{bmatrix}O_{6}\\ I_{6}\\ \end{bmatrix}$, and $U_{e}=\dot{e}_{v,e}+\hat{F}_{e}$.

By inspecting the matrix $A_{e}$ and based on the boundedness of both $\dot{e}_{v,e}$ and $\hat{F}_{e}$, one can find that the state, $X_{e}=[e_{e}^{T},\dot{e}_{e}^{T}]^{T}$, is bounded and exponentially tends to zero as time tends to infinity as soon as the matrices, $S_{c,d}$ and $D_{c,d}$, are positive definite. As a result, since the Jacobian inverse exists (no singularities), the system errors, $e$ and $\dot{e}$, are also bounded and exponentially tends to zero as time tends to infinity.

A complete and detailed block diagram of the proposed control scheme is illustrated in Fig. 5. Quadrotor position and yaw rotation are the controlled variables, while pitch and roll angles are used as intermediate control inputs to achieve the desired $x$ and $y$. Therefore, the proposed scheme has two DOb-based controllers include one for $\zeta=[x,y,z,\psi,\theta_{1},\theta_{2}]^{T}$ (with $M_{n_{\zeta}}$, $P_{\zeta}$, $Q_{\zeta}$) and the other for $\sigma_{b}=[\theta,\phi]^{T}$ (with $M_{n_{\sigma}}$, $P_{\sigma}$, $Q_{\sigma}$). The desired 6-DOF trajectories for the end-effector’s ($\chi_{e,r}$), their actual values calculated by the forward kinematics, and the estimated end-effector force, are applied to the impedance control algorithm, $K_{e}$ that is given in (53). Then, a transformation from task space to joint space is done by using (59) to get $\ddot{\zeta}^{des}$. The desired acceleration in the joint space, $\ddot{\zeta}^{des}$, can be calculated by differentiating (13) with respect to time as

[TABLE]

The desired acceleration in quadrotor/joint space, $\ddot{\zeta}^{des}$, is then applied to the DOb of the independent coordinates, $\zeta$, to produce $\tau_{\zeta}$. The desired values for the intermediate DOb controller, $\sigma_{b,r}$, are obtained from the output of position controller, $\tau_{\zeta}$, through the following simplified nonholonomic constraints relation

[TABLE]

The external controller of second DOb controller, /$\tau_{\sigma}$, is used as a PD controller with velocity feedback, $K_{\sigma}$, as following:

[TABLE]

After that, $\ddot{\sigma}^{des}$ is applied to the second DOb to generate $\tau_{\sigma}$.

It is konwn that the response of the $\sigma$ controller must be much faster than that of the position controller. This can be achieved by the tuning parameters of both DOb and PD of $\sigma$-controller.

The output of two controllers are converted to the required forces/torques applied to quadrotor/manipulator by

[TABLE]

where $B_{6}$ $\in R^{6\times 6}$ is part of $B$ matrix and it is given by $B_{6}=B(3\mathrel{\mathop{\mathchar 58\relax}}8,1\mathrel{\mathop{\mathchar 58\relax}}6)$.

4 Simulation Results

In this section, the presented aerial manipulation robot model with the proposed control technique is implemented in MATLAB/SIMULINK.

4.1 Simulation Environment

For a more realistic simulation studies, the following setup have been made:

• •

Linear and angular position and orientation of the quadrotor are available at rate of $1$ KHz. In [45], a scheme is proposed to measure and estimate the vehicle (Asctec Pelican Quadrotor) states based on IMU and Onboard camera in both indoors and outdoors.

• •

The joints angles are measured at rate of $1$ KHz and angular velocities are estimated by a low pass filter.

• •

The measured signal are affected by a normally distributed measurement noise with mean of $10^{-3}$ and standard deviation of $5\times 10^{-3}$.

• •

$1$ KHz Control loop.

• •

To test the robustness against model uncertainties, a step disturbance is applied at $15$ s to both the inertia matrix, $M(q)$, and the control matrix, $N$, (Actuators’ losses) with $10\%$ error.

The desired trajectories of the end-effector are generated to follow a circular helix, while its orientation follows quintic polynomial trajectories [36]. Parameters of the proposed algorithm are presented in Table 2. The controller is tested to achieve task space trajectory tracking under the effect of the contact force, wind disturbances, and measurement noise.

4.2 Estimation of the End-effector Contact Force

Fig. 6 shows the response of the proposed algorithm to estimate the environment effect/end-effector contact force. From this figure, it is possible to recognize that the norm of the end-effector generalized force has maximum value of $1$ N/N.m at the beginning of operation due to the time taken by the DOb to estimate the system dynamics and external forces. This initial time is about $3$ s. After that period the error norm decreases gradually. The norm of estimation error in both $x$ and $y$ directions have maximum value of $0.03$ N with sinusoidal shape which is due to the sinusoidal motion in theses axes. While the norm value in $z$ direction have maximum value of $0.005$ N. In both $\phi$ and $\theta$ directions, the norms reach value of $0.007$ N.m. The maximum value of the norm in the $\psi$ direction is about $0.0015$ N.m. Thus, it is possible to appreciate the estimation performance of the end-effector generalized forces. Consequently, one can contend that the second control objective is achieved.

4.3 Impedance Control

Fig. 7 shows the response of the system in the task space (the actual end-effector position and orientation can be found from the forward kinematics). From this figure, it is possible to recognize that the controller has good tracking of the desired trajectories of the end-effector (i.e., the tracking error tends to zero as with time). Moreover, it is clear that the capability of the proposed technique to recover the trajectory tracking in the presence of parameters uncertainties which are applied at instant $15$ s. As we see, in the $x$, $y$, $z$, and $\psi$ directions, there is no effect on the tracking. However, in the $\theta$ and $\phi$ directions, the uncertainty effect appears and the controller can recover the tracking quickly.

Fig. 8 shows the motion of the end-effector in the 3D dimension (The markers represent the orientation ). These results show that the proposed impedance motion control scheme provides a robust performance to track the desired end-effector trajectories as well as achieve the desired compliance/impedance effect on the environment taken into consideration the external disturbances and noises. As a result, one can claim that the three control objectives are achieved.

5 Conclusion

The problem of the contact force estimation and impedance control of an aerial manipulation robot is presented with a new solution. A brief presented of the system modeling is given. DOb-based system linearization is implemented in the quadrotor/joint space. A DOb is used in the inner loop is to achieve robust linear input/output behavior of the system by compensating disturbances, measurement noise, and uncertainties. Contact force/environment impedance is estimated based on FTRLS and DOb which appear efficient estimation results and stability guarantee. Then, a linear impedance control is designed and implemented in the task space. The inverse kinematics problem is solved by utilization of the system Jacobian. The controller is tested to achieve trajectory tracking under the effect of external wind disturbances, parameters uncertainty, and measurement noise. Numerical results enlighten the efficiency of the proposed control scheme.