Coupled Dynamic Modeling and Control of Aerial Continuum Manipulation Systems

: Aerial continuum manipulation systems (ACMSs) were newly introduced by integrating a continuum robot (CR) into an aerial vehicle to address a few issues of conventional aerial manipulation systems such as safety, dexterity, flexibility and compatibility with objects. Despite the earlier work on decoupled dynamic modeling of ACMSs, their coupled dynamic modeling still remains intact. Nonlinearity and complexity of CR modeling make it difficult to design a coupled ACMS model suitable for practical applications. This paper presents a coupled dynamic modeling for ACMSs based on the Euler–Lagrange formulation to deal with CR and the aerial vehicle as a unified system. For this purpose, a general vertical take-off and landing vehicle equipped with a tendon-driven continuum arm is considered to increase the dexterity and compliance of interactions with the environment. The presented model is independent of the motor’s configuration and tilt angles and can be applied to model any under/fully actuated ACMS. The modeling approach is complemented with a Lyapunov-wise stable adaptive sliding mode control technique to demonstrate the validity of the proposed method for such a complex system. Simulation results in free flight motion scenarios are reported to verify the effectiveness of the proposed modeling and control techniques.

In the presence of established methods, the dynamics of AMSs with rigid manipulators (AMS-RMs) have already been analyzed [8][9][10] and a number of control methods ( [7,8,[15][16][17][18][19]) have been proposed. However, AMS-RMs have several issues, highlighted in Table 1. For instance, AMS-RMs are not compliant with their environment or shape of the objects, hence requiring exact information about the object shape and position, environmental constraints, and so on [20]. Therefore, AMS-RMs are often operated at a reduced speed to ensure successful grasps [21]. Accordingly, preference is given to manipulators with few DOFs, limiting their capability to compensate for disturbances caused by the vehicle positioning errors. Moreover, external contact forces are avoided as much as possible in order to maintain stability [21]. Such conventional AMSs can hardly be utilized for interactive missions in confined and unstructured environments [22]. As a remedy for the aforementioned issues, AMS with continuum robotic (CR) arms, namely aerial continuum manipulation system (ACMS), [20] (Figure 1) and soft robots [21] have been introduced in the literature. The advantages of CRs and ACMSs have been highlighted in Table 1 and the recent report [20]. For instance, mechanical compliance of ACMS can boost safety when manipulators physically interact with the environment [23][24][25]. Despite the preliminary analysis of ACMSs in [20], there are still many open problems related to their modeling. For example, a previously reported modeling approach [20] for the continuum arm and UAV platform was based on the decoupled dynamics assumption, which decreases the accuracy of the model and introduces certain challenges to the whole ACMS control design. In addition, some simplifications become unavoidable within a decoupled approach when considering the effects of the arm on UAV and vice versa. The focus of this work will be placed on coupled modeling of ACMS dynamics. While coupled modeling of AMS-RMs can be found in the literature ( [16,26]), coupled dynamic modeling of ACMS and its verification remain to be done and constitute the main contribution of this paper.
Generally, there are four methods for dynamic modeling of continuum arms in the literature including those relying on models of Cosserat rod [27], constant deformation [28], center of gravity [29] and Lagrangian [30][31][32]. Our previous published paper on decoupled modeling and control of ACMS [20] was based on Cosserat rod theory. Although this method can be considered as an accurate modeling approach for continuum robots, it is computationally expensive because it needs an iterative solver to satisfy the robot's boundary conditions at each iteration of the program. Therefore, its hardware and software implementation for real time purposes is really problematic. Furthermore, previously proposed Cosserat rod-based controllers in the literature [33][34][35] can only control the robot in x and/or y directions using two tendons per each direction. Thus, the previous method could not handle 3D motion control of the arm. In addition, the previous Cosserat rod-based work can only model the arm. Therefore, Cosserat rod theory cannot be chosen for ACMS coupled dynamic modeling.
As highlighted in [29], control design for Cosserat rod theory is really challenging and only a few papers [33][34][35] have been published to address the control design of the continuum arm using this method. Writing continuum robots formulations to facilitate the control design process is one of the main motivations of introducing a constant deformation method [29]. Recently, some control approaches utilized this method to control continuum arms [36][37][38]. However, we used this method earlier to model the continuum arm and noticed that it is even computationally more expensive than Cosserat rod theory. Furthermore, similar to the Cosserat rod theory, it is presented for arm modeling. Therefore, it can be selected for decoupled ACMS modeling design not for coupled approach. Finally, the center of gravity approach [29] has the same assumption as our method. However, presented formulations for this method are too complicated, resulting in difficulty of adding the effects of the arm on UAV and vice versa. Therefore, it can be concluded that in the first three continuum arm modeling methods, execution time, hardware implementation, real time 3D control and coupling with UAV are problematic. The only remaining method that is also suitable for coupled modeling is the Lagrangian method. However, existing Lagrangian techniques [30][31][32] did not present a general formulation and mostly used Taylor expansion and many curve fittings to simplify equations. Therefore, an appropriate Lagrangian method is needed to address these CR modelling problems in the first step. Additionally, the coupling effect within an aerial manipulation platform should be considered to derive the whole dynamic model of ACMSs. Our previous work [20] focused on the decoupled modelling and control of ACMSs with quadcopters as the base platform. Within this research, we concentrated on the decoupled modelling approach, which did not consider all interactions between the arm and its base platform. To solve this limitation and to provide a more accurate and general dynamic model, here we consider coupled modelling of ACMS with a multi-rotorbased platform.  Another limitation of the previous work relates to the type of UAV used in modeling. Commonly used UAVs, in general, and quadrotors, in particular, are non-holonomic and underactuated. These multi rotors with parallel thrusts cannot independently control UAV's position and orientation (pose) simultaneously [39]. Moreover, they cannot be used for dexterous application with arbitrary generalized forces [40]. Therefore, using an under-actuation platform may degrade the performance of the ACMS missions. To tackle this issue, one strategy is to exploit tilt rotors to achieve non-parallel thrust vectors besides having at least six actuators. A fully actuated UAV has capacity to instantaneously resist arbitrary forces and torques [39] besides the ability to hover at arbitrary orientations [41]. For aerial manipulation tasks, resisting any arbitrary wrench with fast response of exerting forces is a significant factor. Since variable tilting propellers may not be fast enough to deal with external wrenches [39], in this work, fixed-angle propellers are selected for aerial continuum manipulation. While most of previous works on dexterous aerial platforms adopted six propellers, we formulated a general model for a dexterous multirotor platform that is independent of the number and angles of propellers.
In this work, the modeling approach is complemented with a control technique to demonstrate the effectiveness of the proposed modeling technique for not only dynamic simulations but also control purposes. It is a commonly held view that AMS models, even for rigid robotic arms, are highly nonlinear, coupled and complex. Additional complication comes from the integration of continuum arm models that are costly for realtime control. As highlighted in Table 1, control of CRs alone is more difficult than rigid arms and CR's modeling and control are currently defined as a pending challenge in the robotics community. Therefore, the coupled control of ACMSs with these high levels of nonlinearities and uncertainties are quite challenging and can broaden a new horizon and subject of research to the robotic community. In this paper, we also developed an adaptive sliding mode method for the coupled dynamic model to enable coupled control of ACMSs. The control emphasis is on the free motion control with no grasped object to enable coordinated position control of CRs and UAV for reaching the defined target. The designed controller enables incorporating several advantages such as low level of complexity, ease of implementation, whole dynamic consideration, proper execution time, stability and accuracy. Finally, the proposed modeling and control techniques are verified through simulations.
The rest of the paper is organized as follows. The coupled kinematic model of an ACMS is presented in Section 2 while Sections 3 and 4 include coupled nonlinear dynamic model of ACMS using the Euler-Lagrange method. In Section 5, the suggested controller is discussed. Simulation results are presented in Section 6. Finally, Section 7 concludes the work.

Kinematic Modeling
A tendon driven continuum arm attached to the center of a multi-rotor UAV is considered as the ACMS ( Figure 2). The CR arm contains an elastic continuum backbone, a number of tendon spacer disks, and a number of tendons, which can be in the form of cables or wires. The spacer disks are fixed to the elastic backbone and the pinholes on them are required to actuate the tendons. By pulling the tendons, the backbone bends toward the contracted tendon.
To describe the ACMS motion, different coordinate frames are introduced ( Figure 1), with the inertial reference frame ℱ : , , having the origin O I while ℱ : , , possessing the origin indicates the body frame attached to the center of mass of the multi-rotor.
The position of the multi-rotor with respect to ℱ is indicated by , while its attitude is described by Euler angles, = [ ] containing roll, pitch, and yaw angles, respectively. The relationship between the linear velocity of multi-rotor with respect to ℱ , , and its linear velocity with respect to ℱ , , is defined as: where is a rotation matrix of ℱ with respect to ℱ . The angular velocity of the multirotor with respect to ℱ , , can be determined based on the derivative of Euler angles, i.e., = .
(2) Figure 2. A hexa-rotor ACMS with a tendon-driven CR Using a rotation matrix and the relationship between angular velocity in the ℱ coordinate system and the derivative of Euler angles given in Equation (2), the angular velocity of the multi-rotor with respect to ℱ , , can be found as where = . In addition, = [ ] is the angle vector for the arm's end-effector with respect to the UAV. As shown in Figure 1, the end-effector's position can be found using two angles, and , where is the angle of the end-effector with respect to the in the UAV's horizontal plane ( plane), and denotes the angle of CR's endeffector with respect to the perpendicular to the UAV's horizontal plane. Thus, the generalized states of the ACMS, , can be expressed using UAV's position and attitude and CR's end-effector angles as: To calculate the absolute position of each point of the continuum arm at section s in ℱ (where ∈ [ 0 ] denotes the arc length of each section and is the length of CR), the following equation holds, indicating that the absolute position of that special point of the arm equals the summation of the multi-rotor's position in ℱ and the relative position of that point of arm with respect to the center of mass of the multi-rotor in ℱ , where is the relative position of each section (point) of arm with respect to ℱ , which can be found from forward kinematics of the continuum arm [42]. By calculating the position Jacobian matrix at each section ( ( )), the linear velocity of the arm at section s in ℱ can be determined as: In addition, the rotation matrix between the arm point at section s and the multirotor, , can be found in [42] and expressed as [30] = [ ], where , and are columns of the rotation matrix of . Having the rotation matrix, the angular velocity of each section of the arm in ℱ , , is calculated as follows [30] = ( ) = ( ) .
Here, ( ) is the contribution of the Jacobian matrix relative to the angular velocity of the arm at each section in the ℱ and Skew is the skew symmetric matrix of the vector. Finally, the absolute linear velocity of the arm for each section, , can be calculated by taking the derivate of Equation (5) as: Similarly, the angular velocity of the arm for each section in the ℱ frame is demonstrated in Equation (10), based on the angular velocity of the multi-rotor and relative angular velocity of the arm with respect to the multi-rotor.

Dynamic Modeling
To derive the coupled dynamic model of ACMS, the Euler-Lagrange formulation must be considered. The Lagrangian, , is defined based on the total kinetic energy, , and potential energies, , of the system as [43] = − .
Since the damping effect in continuum arm modelling is not negligible, the modified Euler-Lagrange equation [43] is applied as: Total kinetic energy of the ACMS, , is composed of the contribution of three kinetic energy components related to the multi-rotor, , main backbone, . , and disks, . The total kinetic energy and its components are calculated in Equations (13)- (16). In Equation (14), the kinetic energy of the multi-rotor is computed, while in Equation (15), that of the main backbone for each section is presented. Therefore, total kinetic energy for the backbone is determined by integration of Equation (15) along the length of the continuum arm. From Equation (16), the kinetic energy for the j-th disk of the continuum arm is obtained. As a result, the total kinetic energy for whole disks is the superposition of energy for each disk. The equations in Equation (17) indicate that position, rotational velocity and rotation matrix of each disk are the same as those of each section of the continuum arm where disks are fixed to the main backbone. Finally, total kinetic energy of the whole system is determined in Equation (18).
Here, is the mass of the multi-rotor, while denotes the mass of each disk. In addition, , and are angular velocity, position and rotation matrix of each disk, respectively. Here, and present the moment of inertia for the main backbone at each section and that for each disk, respectively, while denotes the moment of inertia for the multi-rotor. Moreover, denotes the density of the main backbone at section , while is the cross-section area for each section of the continuum arm. Furthermore, is the number of disks and ℎ is the distance between disks. Here, the assumption is that the distance between all disks is equal.
By considering Equation (1), Equations (9) and (10), and Equations (13)- (18), the total kinetic energy of the ACMS determined in Equation (18) can be expressed as where is the positive definite symmetric inertia matrix with elements as follows. where × is the identical matrix of size three and = . Similar to the kinetic energy, total potential energy, , reserved in the system is composed of the contributions of potential energy related to the multi-rotor, , main backbone, . , and disks, , as follows.
The potential energy for the multi-rotor is calculated as Equation (21), while potential energy related to the main backbone and disks are given in Equations (22) and (23), respectively. Main back bone potential energy consists of the gravitational potential energy and its elastic potential energy [30].
where is the gravity acceleration and denotes Young's modulus of the main backbone, while is a 3 by 1 unit vector indicating the direction of gravity acceleration with respect to the ℱ .
Finally, having the total kinetic energy (Equation (19)) and total potential energy (Equation (20)) of the system, and using the Euler-Lagrange equation (Equation (12)), the whole dynamic of the ACMS can be written as where is the generalized force, while is the external force exerted on the system. In addition, ( ) and the elements of matrix are given as [22] ( ) = ( ) , The only remaining term in Equation (24) is matrix , which can be determined based on the actuation model. Here, = [ ] , where , and are control signals for UAV position and its attitude, and tension of tendons for the continuum arm, respectively. These terms can be determined as: = , where and are forces and torques for UAV motors in ℱ , respectively, while is the tension vector for tendons of the continuum arm. These matrices and matrix will be defined in the next section.

Actuation Modeling
In this section, an actuation mechanism for a general n-rotor UAV with an attached tendon-driven continuum arm is presented.
As discussed in the introduction, 6 DOFs' control of multi-rotors can be achieved using a nonparallel thrust configuration. In such configurations, each propeller is rotated a cant angle, , around its radius. As shown in Figure 1, is the angle between the normal vector of each propeller plane and UAV's vertical plane. These rotors are laid out along the edge of the disk canted tangentially to the edge of the disk [41]. This means that motors with odd and even numbers are rotated toward the left and right side of their radii, respectively. The position of n motors and their rotations and forces in the horizontal plane of the UAV are defined in Figure 3. In the body reference frame, it is assumed that the n th motor is located at the front of the UAV on the axis of the body frame. The other motors are located at angles , 2 , 3 , …, ( − 1) with respect to the axis. Angle depends on the number of motors ( ) and is given by = . To control the torque applied to the system, motors with odd numbers are rotating clockwise, while motors with even numbers are rotating counterclockwise [44]. Torque produced by each motor can be expressed as = , where is the thrust of the i-th motor, and k is the ratio of the torque to thrust produced by each motor.
To compute the net force/torque acting on the UAV from all thrusters, first each motor's thrust and torque is decomposed into components on the body frame as given in Equation (30). Then the total force/torque for motors can be obtained [44].
where is the cant angle from vertical, and represents the rotor's position in the horizontal UAV frame and equals . In addition, d is the distance from the rotor center to the central axis of UAV (the radius of the n rotor). Moreover, is a sign parameter that is +1 for even motors and −1 for odd motors. If all cant angles ( ) are equal to zero, it means that all motors are perpendicular to the UAV horizontal plane and the generalized model will reduce to a conventional multi-rotors model.
After decomposition of force/torque produced by each motor, the net force and torque of the n-rotor UAV can be computed as = ∑ and = ∑ , respectively. The only remaining term is , which is discussed in the following lines. Here, the assumption is that we have three tendons for actuation of the continuum robot. The first tendon's actuator is located on the of the horizontal body reference frame of the multirotor and two others are located with ± radians with respect to the first one. To control the end effector, two tendons out of three should be actuated at the same time depending on the angle [31]. Table 2 shows the tendon's actuation policy. The displacement of each tendon with respect to the main backbone can be found by [32].
where is the distance between each tendon actuator and the main backbone. If each tendon is actuated with tension , = 1,2,3, then the generalized force can be obtained depending on which tendon is activated based on Table 2, i.e., if tendons 1 and 2 are actuated, one can obtain Finally, generalized force for tendons, , can be rewritten as

Control Design
The formulated controller is based on the adaptive sliding mode approach, which considers the whole dynamics of ACMS and is robust against modeling uncertainties and external disturbances.
In the control design process, the error vector , sliding surface , and the reference signal rate are considered as: where and are the desired states for ACMS and their desired derivatives, and is a positive diagonal matrix. With these definitions, the control signal can be proposed as: where ^ shows the nominal values of matrices while and are positive definite gain matrices. In addition, represents the estimated uncertainty, which is calculated using the adaptation law by updating as: where is a positive definite gain. Substituting the adaptive sliding mode control law of Equation (39) into the ACMS dynamics Equation (24), the closed loop dynamics of the system can be found as: where ∼ is the error on the modeling and uncertainty estimation. Stability analysis of the formulated control for ACMS with coupled dynamics can be conducted using a Lyapunov candidate function as follows.
The derivative of the Equation (42) can be found as By using the closed loop dynamics of Equation (41) and the skew symmetricity feature of − , the derivative of the Lyapunov function can be written as Considering the definition of sliding surface, the above equation can be rearranged as Substitution of the adaptive law of Equation (40) into the previous Lyapunov derivative with the assumption that changes slower than and consequently, ≈ ,one can obtain Being a positive definite of Equation (42) and a negative semi-definite of Equation (45) means that , and are bounded. Therefore, is bounded because = − . Furthermore, using Equation (41), it can be shown that and = − are bounded as well. Taking derivation of Equation (45), can be found as = − − 2 . Since all terms of the right-hand side of the equation are bounded, is bounded and is uniformly continous. Therefore, according to the Barbalat's lemma, converges to zero. As a result, , and converge to zero. According to Equation (41), if remains regulated by the controller and remains close to zero, it can be concluded that approches zero.

Experimental Simulation Results
In this section, the developed coupled dynamic modeling approach for the considered ACMS configurations (a fully actuated octa-rotor [45] with a tendon driven continuum arm actuated by three tendons [20]) is validated using experimental simulation results to verify dynamics and control methods. In this regard, first, the dynamic model of the coupled ACMS is derived using the MATLAB 2020a symbolic toolbox. Then, the ACMS model and control are implemented in MATLAB 2020a, while sampling time is chosen as 0.01 s and a 5th order Runge Kutta solver was used.
The objective of the experimental simulation scenario is to verify the performance of the control method designed for the coupled system dynamics in free flight motions when there is not any interaction with the environment. For this purpose, similar to [22], we analyzed two scenarios: (1) tracking a helix trajectory while the continuum arm is fixed, and (2) tracking an agile infinity trajectory while the arm moves.
In these two control cases, an octa-rotor with the capability of carrying 3 motors required for arm actuation, i.e., Dynamixel MX-28-AT geared servo motor (ROBOTIS, Lake Forest, CA, USA), each with a mass of 91 gr was considered. The selected octa-rotor weight was 3.2 kg [45] with a cant angle of 25 deg and the selected continuum arm specifications were the same as [20].
For the first case, we assumed that the arm was fixed in the desired position with respect to the UAV ( = 2 ⁄ rad and = 4 ⁄ rad) while the UAV tracks a helical trajectory of 4 − 4 4 0.2 (all in m) with zero desired Euler angles. Since a fully actuated UAV is selected, it should be able to control its position and attitude independently. That is why desired Euler angles are considered zero.
In the second case, UAV was driven to follow an agile infinite path of 2.5 2.5 2 (all in m) with zero desired Euler angles while the end effector of the arm had a circular motion of ( = rad and = rad) with respect to the UAV.
The initial condition of the arm was = rad and = rad in all cases while the initial condition of UAV in the first case was [0 0 0 0 0 0] (in and ) and that of the second case was [0 0 2 0 0 0] (in and ) . A minimal pose representation with position and RPY Euler angles was used here. Figure 4 illustrates the trajectories of the drone for these two simulation cases.
(a) (b) For the first control case, the helical trajectory with fixed arm angles, the states of ACMS and their desired values are shown in Figure 5. In this figure, ACMS states are shown by blue lines while red lines illustrate desired values of states. As shown in Figure  5, ACMS effectively tracked desired values. Figure 6 presents tracking error for this case, which is around 0.1 m for position and zero for attitude tracking. Octa-rotor thrusts to track the given helical trajectory and required tension of tendons to reach their fixed desired values are illustrated in Figures 7 and 8, respectively. Figures 7 and 8 show that UAV motor thrusts had an oscillatory behavior around 5 N to produce the required force and moments in all directions to track the desired helical trajectory with zero Euler angles. In addition, only tendons 1 and 2 were actuated because is between 0 and 120 deg. Figures 9-12 present those results for the second control case, i.e., an infinity trajectory with a moving arm. Figure 9 illustrates ACMS states and their desired values. As it can be concluded from Figures 9 and 10, similar to the previous case, ACMS tracks its desired values with a reasonable error. The error for tracking desired X is around 0.1 m while those of Y and Z directions are around 0.2 m and zero, respectively. Like the previous simulation scenario, UAV's motors have to oscillate to generate required thrust to track the given infinity trajectory since desired Euler angles are considered zero. However, motor numbers 6 and 8 oscillate less than other motors according to Figure 11. Finally, required tensions of tendons to track a circular motion are presented in Figure 12. As expected from the tendon actuation policy (Table 2), only two tendons are activated at the same time depending on values. As shown in two control scenarios, Figures 5-12, the proposed control method can successfully control the ACMS with the coupled dynamic model. Although the coupled model is more complicated than the decoupled model, the simulation took 0.02 s per iteration using a PC with an intel core i7-8700 processor, 3.2

Z (m)
GHz, and 16 GB RAM, within MATLAB 2020a, which is 10 times faster than decoupled model execution time. This is in large part due to the fact that the coupled model is not based on the iterative solver used for the decoupled approach. For real-time implementation, a real-time computing language such as C must be used, which could reduce execution time up to 500 times [46].

Conclusion
This paper presents the first coupled dynamic modeling for aerial continuum manipulation systems (ACMSs), which considers both UAV and the continuum arm as a unified system. The proposed model was developed using a modified Euler-Lagrange framework to take into account damping and elasticity effects of the arm. The approach is enriched by considering a general tiltrotor UAV with fix tilt angles as a base platform, which can model both under-actuated and fully actuated ACMSs. The integrated continuum arm attached to the underneath of the UAV is a general tendon-driven type with three tendon actuators. The model can be extended to different continuum arm configurations such as series, parallel and collaborative types. Furthermore, in the future, the actuating mechanism can be replaced with other methods such as pneumatics actuators in concentric tube continuum arms. A nonlinear adaptive sliding mode control method was also formulated to assess the behavior of a fully actuated ACMS model with an octa-rotor base platform in free flight missions when there is no interaction with the environment. The stability of the proposed controller was proved using Lyapunov stability theorem. Experimental simulation results were reported to verify the proposed dynamic modeling and control approaches. Results taken from a developed ACMS simulator in MATLAB environment demonstrate the impressive performance of the controller in tracking desired helix and infinity trajectories of the UAV while the arm is fixed or moving. It was also shown that the proposed coupled model is ten times faster than the previously reported decoupled model for ACMSs. Our future work will include extending the present work to interactive aerial missions like impedance and hybrid position-force control of the system.