I Introduction
Over the past decades, aerial manipulation has received great attention in the robotics research community, with many different systems in use [Ruggiero2018, Bonyankhamseh2018]. The tasks solved by aerial manipulators range from grasping, fetching, and transporting arbitrary objects to pushing against fixed surfaces. Potential use cases are numerous: inspection of infrastructure like bridges or manufacturing plants [Ikeda2017, Trujillo2019, Bodie2019], physical interaction through tools like grinding, welding, drilling and other maintenance work in hardtoreach places [Papachristos2014, Bodie2019], and the autonomous pickup and transport of objects [Kessens2016, Augugliaro2014]. All of these tasks require the MAV to be equipped with an additional mechanism which we refer to as the end effector. The coupling of the MAV dynamics with the moving end effector poses an interesting challenge from a control perspective given the inherent instability of MAVs.
The requirements for high precision in real world aerial manipulation applications further increases the difficulty of the task. Efforts made into this direction include [Darivianakis2014, Bodie2019], where the authors use a fixed end effector on a underactuated and an omnidirectional MAV respectively. Compared to the former method, our approach achieves higher accuracy and flexibility in terms of potential use cases due to our moving end effector. In comparison to the second approach, we achieve on par precision while relying on a simpler, underactuated platform. In summary, we claim to show the following contributions:

We present a hybrid model which captures the nonlinear dynamics of the MAV and considers the quasistatic forces introduced by the attached manipulator.

We use this generic model in an NMPC jointly controlling the MAV and arm motion.

We experimentally evaluate our method in ‘aerialwriting’ tasks using our custom built system. Our results demonstrate high repeatability and accuracy in the order of millimetres across multiple trajectories of varying difficulty.
This paper is organised as follows: In Section II we give an overview of the related work on aerial manipulation. In Sections III and IV we describe our notation, geometric arrangement, and the software architecture of our system. We explain the method in detail in Section V followed by the experimental results in Section VI. Finally, in Section VII we discuss our findings and conclude in Section VIII.
Ii Related Work
Aerial manipulation systems can be broadly distinguished based on the MAV type (as being omnidirectional or underactuated) and the end effector (as being fixed or moving). In general, using an omnidirectional MAV to fulfill complex aerial tasks does not require a moving end effector as the necessary 6 Degrees of Freedom are provided by the MAV itself. Examples include the works presented by [Brescianini2016, Brescianini2018, Ryll2019]. Brescianini2016 show an omnidirectional MAV that achieves 6DoF motion by using eight fixed rotors in a noncoplanar configuration. In a subsequent study [Brescianini2018], the same platform is used with a fixed end effector to fetch moving objects. Using a similar approach with a fixed configuration of six tilted rotors, Ryll2019 propose a novel paradigm to control all 6 DoFs of the MAV while using a rigid end effector to exert forces and torques independently. The system is demonstrated in numerous experimental tasks that include contact, e.g. surface sliding and tilted peginhole task.
In a different approach [Kamel2018], a setup consisting of six rotors which actively control the thrust direction is proposed. Bodie2019
leverage this system to solve a variety of aerial manipulation tasks with a rigidly mounted, low complexity end effector. The authors further show precise force control when in contact with unstructured environments while relying on visualinertial state estimation. While this platform allows for accurate 6
DoF flight and longitudinal force exertion with a relatively simple control method, it is mechanically more complex and thus more costly compared to classical multicopter platforms. Recently, a similar approach was followed by Trujillo2019, who introduce AeroX, an omnidirectional octacopter for contactbased inspection. Their end effector design minimises the torque caused by contact forces and features wheels on its base to allow moving along a surface while remaining in contact. In [Papachristos2014, Papachristos2014ICRA] the authors show a less complex but highly capable tritiltrotor MAV for surface grinding and obstacle manipulation. The control model consists of two disjoint modes: one for freeflight and another for physical interaction.Employing an underactuated MAV to perform aerial manipulation typically increases the complexity of the end effector since the latter has to provide the additional DoFs. To investigate this, many different end effector designs have been proposed over the last years. We categorise these works by the increasing complexity of the end effector. In [Darivianakis2014] an underactuated MAV with a fixed end effector was used for performing contactbased tasks. The authors use a switching mode linear Model Predictive Controller (MPC) with different control models for free flight and incontact operation. Similarly to our work, they benchmark their approach by performing ‘aerialwriting’. In [Mellinger2011], the authors use different lightweight, low complexity grippers to perch, pick up, and transport payload. Using a least squares approach, they estimate the inertial parameters of the grasped objects and use this information to adapt the controller and improve tracking. Moving up in terms of complexity, Kim2013 suggest mounting a 2DoF robotic arm on a MAV to allow grasping and transporting of objects. The authors propose an adaptive sliding mode controller for the combined system. In [Orsag2014] the authors present an aerial manipulator with two robotic 2DoF arms to open a valve. The MAV and arms are controlled as a coupled system which is modelled as a switched nonlinear system during valve turning. In an approach very similar to ours, Lunni2017 propose an NMPC which jointly controls the motion of an underactuated MAV and the attached manipulator. However, the control model is only formulated for free flight operation and the presented experiments do not include contact. In a more recent work, Suarez2018 propose a lightweight, humansized dual arm system designed to minimise the inertia transferred to the MAV. Each of the two arms add 5 DoFs to the system and the arm control law applied takes into account that lowcost servo motors do not allow torque control but require position commands. Further, a torque estimator is used to predict the torques produced by the servos and inform the MAV control algorithm accordingly. In order to minimise such disturbances coming from the end effector, Nayak2018 propose a lightweight design mounted on top of an MAV. While attaching a serial robotic arm on an MAV increases the number of tasks it can perform, they only provide limited precision when using low cost and lightweight actuators. Some previous efforts try to mitigate this by mounting a parallel delta arm on an MAV instead. [Kamel2016_1, Kamel2016_2] demonstrate a multiobjective dynamic controller which also considers dynamic effects between the MAV and its 3DoF delta arm. The same system is used in [Steich2016] to inspect tree cavities with a camera mounted at the end effector.
In our work we use a low mechanical complexity setup consisting of an underactuated hexacopter and a 3DoF parallel arm. We propose a general method which could be transferred to other aerial manipulation platforms. We extend the NMPC framework presented in [Tzoumanikas2020] for manipulation tasks and use a hybrid model that captures the effect of contact and coupled MAVarm dynamics. Our main focus is on the developed software while we use the presented hardware platform to showcase our method. To our best knowledge, we achieve unprecedented accuracy in experiments requiring contact by using an underactuated MAVdelta arm system.
Iii Notation and Coordinate Frames
We denote vectors as bold lower case symbols, e.g.
. We use lefthand subscripts, e.g. , to indicate the coordinate representation in the frame of reference. The rotation matrix changes the representation of the vector from to as . Analogously to the rotation matrix , we use the quaternion with denoting the quaternion multiplication. We useto denote the skew symmetric matrix of the vector
. The motion of the MAV body frame (x: forward, y: left, z: upward) is expressed with respect to the World frame (z: upward).Iv System Overview
An overview of the different software components of the proposed system is outlined in Figure 2. The NMPC is given full state trajectory commands for the MAV and the end effector corresponding to a given aerial manipulation task. Based on those references and the estimated system state, it produces the desired MAV
body moments, collective thrust, and end effector position. The control allocation block is responsible for converting the moments and thrust into individual motor commands while the inverse kinematics block computes the desired link angles for the given end effector position. All algorithms run onboard at a rate of 100
while the estimate of the MAV position and orientation is provided externally. The different coordinate frames are displayed in Figure 3.V Method
Va Hybrid Modelling
The standard NewtonEuler equations are used to model the combined MAVarm dynamics. We model the MAV as a single rigid body object and only consider quasistatic forces introduced by the arm dynamics and its interaction with the environment. Overall, the combined dynamics take the following form:
(1a)  
(1b)  
(1c)  
(1d)  
(1e) 
where , are the combined MAV
arm mass and inertia tensors, respectively and
the acceleration due to gravity. Regarding the forces and moments , we use the subscript to distinguish the ones generated by the MAV motors from the ones caused by the end effector movement and its potential contact with the environment. In our system, the MAV motorgenerated forces and moments are given by:(2a)  
(2b) 
with the thrust produced by the motor, its position with respect to , the known thrust to moment coefficient and . Equation (2) can be summarised as with the allocation matrix related to the MAV geometry as described in [Tzoumanikas2020]. and are given by:
(3a)  
(3b) 
where is the contact force acting on the end effector expressed in its frame and the nominal end effector position which results in no Centre of Mass (CoM) displacement. The two terms in (3b) represent the moments due to contact and due to the displacement of the CoM respectively. The combined mass is the sum of the MAV and end effector mass, respectively, while the combined rotational inertia can be computed as with the inertia tensor of the MAV (including the arm in nominal position) and the corresponding diagonal matrix.
Regarding the contact force, we assume that this can be approximated with a linear spring model as:
(4) 
where is a known spring coefficient and is the normal component of the contact surface penetration. This way, the controller can anticipate contact before it even happens and there is no need for a switching mode controller (one for free flight and another one for contact dynamics).
VB Model Based Control
For the control formulation we define the following control state and input:
(5a)  
(5b) 
Note that we use for the formulation of the control model, while is used in the control input. We use the constant and known homogeneous transformation to change the coordinate representation of these position vectors.
We use the following error functions for the position of the MAV, the position of the end effector, the MAV linear and angular velocity, the orientation, the contact force and the control input, respectively:
(6a)  
(6b)  
(6c)  
(6d)  
(6e)  
(6f)  
(6g) 
with and the superscript used to denote the timevarying reference quantities. The optimal input sequence is obtained by the online solution of the following constrained optimisation problem:
(7a)  
s.t. :  (7b)  
(7c)  
(7d) 
with the discrete horizon length, the discrete version of the dynamics given in Equations (1) – (3), the latest state estimate and appropriate lower and upper bounds for the control input defined in (5). For the intermediate and final terms we use quadratic costs of the form as defined in (6) where the gain matrices were experimentally tuned.
The optimal control problem is implemented using a modified version of the CT toolbox [ct] with a 10 discretisation step and a 2 constant prediction horizon. We use a RungeKutta 4 integration scheme followed by a renormalisation for the quaternion. As common in receding horizon control, the first input is applied to the system and the whole process is repeated once a new state estimate becomes available. The motor commands for the MAV are obtained by solving the following Quadratic Program (QP):
(8a)  
s.t. :  (8b) 
where are the minimum and the maximum motor thrust, is a gain matrix and a regularisation parameter. The end effector position commands are mapped into servo angle commands by solving the inverse kinematics problem for the delta arm explained in Section VC. In the case of an infeasible (e.g. outside the arm’s workspace) or unsafe end effector position command (e.g. one that results collision between the MAV propellers and the arm’s links), the position command is reprojected onto the boundary of the feasible and safe to operate workspace. In practice, this was rarely the case as the MAV and end effector reference trajectories are designed so that the end effector operates close to its nominal position. In this way the usable workspace is maximised while the effect of the CoM displacement (which is captured by our control model) is minimised.
Our C++ implementation of the above, requires 6.7
with a standard deviation of 0.57
per iteration. On average, 98% of the computation time is spent on the optimisation problems defined in (7) and (8).We would like to highlight that our method is generic enough to be applied to different types of vehicles such as omnidirectional ones and/or other types of manipulators. In these cases, the control input, the control allocation and the arm kinematics have to be adapted based on the vehicle and manipulator type. Similarly, the model can easily be extended to capture aerodynamic friction, gyroscopic moments, handle multiple contact points, or use more sophisticated contact models (e.g. ones that include a combination of linear springs and dampers). Similarly, the writing task which we use for the experimental evaluation, is just an example application that requires precision. We believe that our algorithms are adaptable to other tasks such as inspection through contact.
VC Delta Arm Kinematics
Our MAV is equipped with a custom built 3DoF delta arm [delta_arm_paper]. Its main advantages are speed, as its few moving parts are made of lightweight materials, precision and the easy to solve forward and inverse kinematics. The forward kinematics problem (i.e. determining the position of the end effector given the joint angles ) can be solved by computing the intersection points of three spheres (shown in Figure 4) of radius with the following centres:
where correspond to the arm physical parameters shown in Figure 4, , and rotation matrices of and degrees around . The maximum number of intersection points is two which corresponds to an end effector position above () and below () the arm base. In our setup, solutions above the arm base are mechanically impossible and thus rejected. For the inverse kinematics the intersection between a sphere with radius and a circular disk with radius has to be computed for every joint angle. For the first joint, as shown in Figure 4, the centre of the sphere is with while the centre of the circular disk is . Given the intersection point , the joint angle can be recovered as . The joint angles and can be computed by performing the same procedure for the spheres with centres , same radius and the unit disks centred at and with radius . The points can be easily computed as follows:
(10a)  
(10b)  
(10c)  
(10d) 
VD Trajectory Generation
We use a trajectory generator to map arbitrary sets of characters to end effector trajectories. We use a constant acceleration motion model to generate trajectories with a smooth velocity profile. This is of special importance when the reference path contains sharp edges. Through our software we can adjust the velocity and acceleration profile by changing the maximum and . Once the trajectory for the end effector has been computed, we proceed with the computation of the reference position and velocity for the MAV as follows:
(11a)  
(11b) 
The reference MAV orientation is chosen such that the end effector is always perpendicular to the contact surface, assuming perfect position tracking, while the reference force is obtained using the spring model and the nominal displacement of the end effector into the contact surface’s normal direction. In our framework each trajectory is accompanied by an appropriate flag which disables or enables the position tracking for the end effector. This is achieved by setting the appropriate gains to zero. In that case the NMPC may decide to move the arm to assist the reference tracking of the MAV due to the CoM displacement. This potentially unwanted behaviour can be avoided by further penalising the arm displacement from its nominal position (i.e. by increasing the input gains). However, it is an interesting capability enabled by our hybrid modelling.
Vi Experiments
Via Experimental Setup
The experiments presented in this section were performed using a custom built hexacopter equipped with a sideways mounted delta arm manipulator. The MAV features a frame with a 550 diameter, a Pixracer flight controller running a modified version of the PX4 firmware, and an Intel NUC7567U onboard computer running Ubuntu 16.04. It uses 960 KV motors and DJI 9450 propellers. The delta arm uses magnetic universal joints for the connection of the servos with the end effector which maximizes the workspace, minimises backlash and allows the arm to disassemble during possible crashes preventing it from breaking. The arm uses three Dynamixel AX18A servo motors which are comparably fast and accurate but have limited maximum torque. The system is powered by a 4S 4500 mAh battery and has total weight of 2.6 . The end effector holds the pen which is mounted on a spring to provide additional compliance. We set the coefficient of the contact model in (4) to match the used spring. The applied force is measured by a SingleTact force sensor mounted at the end of the spring. We estimate the spring coefficient by measuring the applied force for known tip displacements. The dimensions of the delta arm were obtained from a highly detailed CAD file and were verified manually. We measured the inertia of the MAV by measuring its angular response to constant input torque while it is hanging to freely rotate. The thrust to moment coefficient is measured using a thrust stand. A table with the numeric values of the system parameters is given in Table I and a photo showing the platform and its different components is shown in Figure 5.


We use a Vicon motion capture system to provide external pose estimates. The contact surface is a whiteboard for which we estimate its pose based on Vicon measurements. Each experiment consisted of the following different trajectory stages: (i) approach, (ii) write, and (iii) return home. The end effector was enabled, using the appropriate flags as mentioned in Section VD, for the trajectory writing in (ii) and disabled for the rest. Our analysis mainly focuses on the trajectory writing which includes contact whereas for the other two parts (approach/return) the MAV performs simple position tracking. We evaluate the accuracy of our system by comparing the reference trajectories to those estimated by the Vicon motion capture system. In addition, we use a visionbased error as a performance metric. This is because we observed inaccuracies in the Vicon measurements stemming from either bad calibration, poor object visibility, or marker reflections on the whiteboard surface. The visual error is computed by running a 2D Iterative Closest Point (ICP) method [Besl1992] on a filtered and rectified photo of the final writing and a rendering of the planned path. After registration of the two point sets, we use the nearest neighbour distance to evaluate the accuracy.
ViB Trajectory Tracking
Figure 7 shows the tracking of the RSS trajectory visualised in the contact frame for the end effector and the MAV. The maximum reference velocity was set to and the maximum acceleration to . The trajectory consists of four contact segments with a combined duration of 65 . Based on the Vicon estimates the tracking error, shown in Figure 7, of the end effector almost always remains in the range during the contact segments while the same quantity for the MAV is in range. This highlights the efficacy of using a manipulator with faster dynamics than the MAV’s for precision tasks such as ‘aerialwriting’.
Similarly to the above, Figure 9 shows the trajectory tracking for the more challenging experiment which contains ten contact segments with a combined duration of 63 . Tracking accuracy is similar as before with the end effector and MAV tracking error in the and range. The accuracy can be visually verified since the overlapping segments of the ‘R’ and ‘m’ coincide almost perfectly. Additionally, the consistent approaching and retracting from the contact surface leads to identical starting points of individual letter segments, e.g. the three horizontal lines of the letter ‘E’. In both cases, the maximum error based on the visual error analysis is 10 mostly originating from temporary loss of contact. Possible reasons for this are bad estimation of the orientation part of the contact frame transformation , the assumption of a perfectly flat contact surface being incorrect but also the finite accuracy of the delta arm. The imperfect tracking along the contact frame normal direction (shown in blue in Figures 7, 9) is also reflected in the reference force tracking.
To prove the repeatability of our approach, we conducted each experiment thrice. We give the relevant tracking statistics for the MAV and arm separately in Figure 10, in which the textured box plots correspond to MAV data and the plain ones to that of the end effector. The median and extreme values for the end effector are significantly lower than the ones for the MAV and consistent with the values reported above. This further shows the benefit of using an aerial manipulator for precise tasks including contact. The MAV tracking accuracy along the axis was the lowest amongst all axes, as this was most affected by the interaction forces and unmodelled torque disturbances due to the movement of the servos.
ViC Velocity Sweep
The aim of the next experiments is to demonstrate the effects of the input velocity and acceleration on the writing accuracy. We performed five iterations of the same Hello trajectory experiment with different velocity and acceleration profiles. The different iterations correspond to maximum velocity and maximum acceleration .
Figure 11 shows the box plots for the MAV and end effector tracking accuracy based on the Vicon measurements. The plots show consistent tracking results in all the different velocity and acceleration settings tested. The numeric values of the tracking error are similar to the ones presented in Section VIB with the end effector achieving sub centimetre accuracy (per axis) while the MAV error is consistently less than . However, by observing the visual error, we see that as the reference velocity increases, the system struggles more with the segments containing curvature e.g. ‘e’ and ‘o’. In contrast, the performance on the straight line segments remains similar.
We believe that the tracking error of the MAV can be further reduced by giving the NMPC dynamically feasible trajectories not only for position and velocity but also acceleration, jerk, and snap. Regarding the end effector tracking error, we generally expect this to increase for reference velocities beyond the ones tested here. This is due to the mismatch between the formulated control model and the real one.
ViD Text Size Sweep
In Figure 12 we show the visual error of our system for the same trajectory in four different text sizes ranging from 10 to 40 . The consistent accuracy shows that the system can handle the fast direction changes imposed by the small scale.
Vii Discussion
Overall, our system achieves accurate and consistent results over a series of different trajectories. The tracking error of the end effector is significantly lower than that of the MAV; highlighting the accuracy boost due to the utilisation of the arm. We would like to mention that our system was built using relatively lowcost offtheshelf components and 3D printed parts. This leads to errors in the manufacturing with respect to the reference model, e.g. errors in the true inverse kinematics of the arm due to nonidentical dimensions of its links.
Another important issue that we faced during our experiments was the reliance on the motion capture system for localisation. Apart from issues related to WiFi delays, which resulted in temporary loss of tracking, we faced problems with poor object visibility resulting in unreliable estimates of both the static objects, such as the contact frame, and moving ones such as the MAV. In fact, during our data analysis we realised that there are segments where Vicon returned mechanically impossible configurations for our system e.g. end effector position below the surface of the contact frame. Despite these problems which further propagate into tracking errors, our system successfully handled multiple transitions to contact during the same experiment.
We experimentally verified that for contact tasks, where the main objective is accuracy instead of speed, using a planner respecting full state dynamic feasibility is not an absolute necessity. Despite our simplified motion planner, our system achieves subcentimetre accuracy. However, we argue that for more aggressive maneuvers, a full state dynamically feasible planner would be required.
Viii Conclusion
We have presented a hybrid modelbased algorithm for aerial manipulation using an underactuated MAV with an attached end effector performing ‘aerialwriting’ tasks on a whiteboard. We demonstrated our system in a series of experiments with trajectories requiring multiple transitions from freeflight to contact and viceversa. The end effector tracking error consistently remains in the range for trajectories with maximum velocities ranging from 7.5 to 27.5 , maximum acceleration from 3.75 to 13.75 and text sizes from 10 to 40 . All algorithms run in realtime onboard the MAV. We believe that our method is generic enough to be applied to different types of MAVs and manipulators. We further believe that our framework can be extended to more practical applications such as physical interaction with surfaces e.g. drilling, welding, grinding or inspection through contact.
Regarding future work, online estimating and closing the manipulation task loop with visual feedback would constitute a major improvement to correct for errors stemming from bad calibration, Vicon delays and inaccuracies of the arm. This will allow us to better handle the system’s sensitivity to imperfect estimates of transformations such as the contact frame and arm to body transformation . The current method further assumes an instantaneous end effector position response, which we plan to substitute with a more realistic model or ultimately formulate the full multibody dynamics model which would result in a performance boost mainly for faster maneuvers. Finally, we aim to implement and test our method using omnidirectional MAVs which due to their ability of generating lateral forces have superior force exertion capabilities compared to the one presented in this work.
Comments
There are no comments yet.