Three-Dimensional Modeling and Control of a Tethered UAV-Buoy System

This work presents the nonlinear dynamical model and motion controller of a system consisting of an unmanned aerial vehicle (UAV) that is tethered to a floating buoy in the three-dimensional (3D) space. Detailed models of the UAV, buoy, and the coupled tethered system dynamics are presented in a marine environment that includes surface-water currents and oscillating gravity waves, in addition to wind gusts. This work extends the previously modeled planar (vertical) motion of this novel robotic system to allow its free motion in all three dimensions. Furthermore, a Directional Surge Velocity Control System (DSVCS) is proposed to allow both the free movement of the UAV around the buoy when the cable is slack, and the manipulation of the buoy’s surge velocity when the cable is taut. Using a spherical coordinates system centered at the buoy, the control system commands the UAV to apply forces on the buoy at specific azimuth and elevation angles via the tether, which yields a more appropriate realization of the control problem as compared to the Cartesian coordinates, where the traditional x-, y-, and z-coordinates do not intuitively describe the tether’s tension and orientation. The proposed robotic system and controller offer a new method of interaction and collaboration between UAVs and marine systems from a locomotion perspective. The system is validated in a virtual high-fidelity simulation environment, which was specifically developed for this work, while considering various settings, operating conditions, and wave scenarios.


Introduction
Teams of unmanned aerial vehicles (UAVs) and unmanned surface vehicles (USVs) harness the advantages of each vehicle to form a superior robotic system. To face challenges brought by flood disasters [29], marine oil spill events [10,22], search and rescue [4,21], monitoring and patrol [9], and water surface cleanup [11], amongst others, heterogeneous UAV-USV systems have been employed.
The advantages of combining UAVs and USVs have been explored in literature [14,29] based on their complementary abilities. UAVs are advantageous given their wide field of vision and bird's-eye view, higher maneuverability, flexibility, and ease of deployment; on the other hand, USVs are advantageous in their durability, extended running time, and large loadcarrying capacity. For instance, to deal with flooding disasters, the advantage of a UAV's wide-angle view of the environment, including the visual horizon, was combined with the long cruising ability of a USV to overcome their individual shortcomings, thus allowing the creation of a system with higher efficiency and lower risk in Search and Rescue (SAR) operations [29].
UAVs and USVs can form a dexterous robotic team in the marine environment, as evidenced by the various applications considered in the literature. A heterogeneous robotic swarm was considered in [21], where UAVs are used to generate and transmit pathplanning information in addition to full area maps for USV agents to perform rescue missions. Further cooperation between the two vehicle systems can be achieved via additional sensing and control schemes, among which is the aerial visual tracking of the USV. One control solution for robust visual tracking was proposed in [9], and going further, even visual pose estimation of the USV is made possible using the UAV's on-board camera as proposed in [4]. If they get within close proximity, UAVs can perform coordinated trajectory tracking of USVs as in [26], and they can target-track them to prioritize the camera tilt over the UAV's motion while maintaining continuous monitoring of the USV [17]. Physical interaction is another aspect of the heterogeneous UAV-USV system. For instance, the USV can serve as a landing platform for the UAV [23], power transmission from a USV to a UAV can be achieved through an umbilical power cable [24], and even floating object manipulation can be made through a tethered UAV system [14].
In addition to using robotic manipulators, a UAV can interact with its environment via a tether [19,20], as it is flexible, extendable, light-weight, and can transmit tensile forces. Recently, there have been several advances and studies on the tethered UAV problem. A tethered power system for UAVs was proposed in [2] and is now commercially available [7]. Such a system can be optimized to suit the special case of a continuously oscillating marine power station as studied in [24]. A tethered UAV is also useful in information transmission. For instance, it plays a keys role in the emergency marine communication network proposed in [28], where it takes off from the communication support ship to secure a wide network coverage, while being fed with bandwidth and power to sustain its location. To address the challenge presented by wind gusts in the open marine environment, which can affect the stability of tethered systems, few works can be found on the effect of wind disturbances on the tether and the UAV's stability [19], and on the tether vibration [27]. To simplify the analysis of the tether dynamics, the system was modeled as a multi-element body in [3].
Beyond the use of a tether for power or information transmission, it can transmit tensile forces whereby a UAV can control the magnitude and direction of the link/tether force, which was proven to be a set of differentially flat outputs of tethered UAV systems in [25]. A practical example of force transmission is found in [18], where the transportation of a payload is cooperative performed by a team of UAVs, which was made possible through calculation of the required wrench set. The clever use of a tether even allowed a single UAV to lift a heavy object after hooking itself at high ground then activating its on-board winch [12]. Furthermore, a UAV was employed in [14] to manipulate a floating object's surge velocity by means of tether in a 2D space (vertical plane), even in the presence of waves and water currents.
Amongst the various applications on a forcetransmitting tethered UAV, the emerging marine locomotive tethered UAV−buoy system introduced in [14] offers a new type of interaction between a UAV and floating objects including USVs, which provides a force interconnection between the elements of the system in addition to power and information transmission. The potential of this marine robotic system is yet to be fully discovered, and features a wide range of applications such as search and rescue, floating sensors manipulation, water surface cleanup, and building and inspecting marine structure, to only name a few.
The 2D-planar model of the tethered UAV−buoy system that was proposed in [14] presented a thorough analysis of the system's motivation, its potential application, the subsystem's dynamics and the systems' stability, the effect of waves, the proposed working bounds, steady-state values of the system states, and several practical considerations for implementing the system. However, the 2D tethered UAV−buoy model does not capture the complete dynamics of the real system. For instance, the 2D model has limited representations of the Euler angles of both the UAV and the buoy, it does not encompass composite waves and currents with varying directions, and it does not allow the system to follow trajectories that are in the 3D space. In addition, the orientation of the buoy is not captured in the 2D-planar model, which can deviate from the actual direction of the tether. Finally, the designed Surge Velocity Control System (SVCS) in [14] does not include the tether's azimuth angle within its controlled state variables, thus the controller is unable to provide directional manipulation of the buoy in the horizontal plane if deployed on a 3D model. Therefore, for real world application, it is necessary to extend the 2D planar model of the tethered UAV−buoy system to the 3D space, so that it includes all position and orientation states of the UAV, the buoy, and the tether, as well as the surface water model.
The contributions of this paper are presented next. First, we extend the model of the 2D-planar tethered UAV−buoy system dynamics to the 3D space by incorporating the full six degrees-of-freedom (6-DOF) rigid-body model for each of the UAV and buoy, and a 2-DOF model of the tether, to arrive at a 11-DOF dynamical model that allows for a more realistic representation of the system's physics. Second, the ocean/sea environment is modeled to include not only the surface water current, but also the full effect of oscillating gravity waves, even in the vertical direction. Third, a rigorous mathematical derivation of the composite system is presented through the Lagrangian formulation, while defining the required conditions and constraints to regulate the system's motion. Forth, a spherical coordinates-based dual controller is designed to allow the free movement of the UAV around the buoy while the tether is slack, and to manipulate the buoy's surge velocity while the tether is taut. This controller allows the UAV to apply directional tension through the tether on the buoy so that it can trigger motion in all directions of the water surface.
The remainder of this paper is organized as follows. Section 2 presents a multi-physics description of the 3D tethered UAV−buoy system components and the water environment. Section 3 introduces the control system design for controlling the relative UAV−buoy position and the buoy's surge velocity. Comprehensive numerical co-simulation results of the derived system model and designed control system are presented in Section 4. Section 5 closes the paper with with a conclusion and identifies future development tracks.

Tethered UAV−Buoy System Dynamical Model
The heterogeneous marine robotic system of a tethered UAV−buoy system has multi-physics elements, which must be integrated with a proper definition of the marine environment that it operates in. This section defines the system components and their interconnection including the water environment, the USV/buoy, the tether, and the UAV, and it finally Figure 1 Three-dimensional model of a tethered quadrotor UAV−buoy system with tensile force interaction in a marine environment.
culminates with the coupled dynamical model of the tethered UAV−buoy system.

Preliminaries
Before defining the problem, we briefly specify general notations to support the model derivation process. For some angle (•), we let c • , s • , and t • respectively be the cosine, sine, and tangent functions. Also, we define R >0 to be the set of positive-real numbers such that {x ∈ R | x > 0}, and R ≥0 to be the set of non-negative real numbers such that {x ∈ R | x ≥ 0}.

Problem Definition
Consider the 3D space above the water surface in which the problem is defined as shown in Fig. 1. Let O I be the origin of the inertial frame of reference, W = {x, y, z}, located at the horizontal plane of the local mean sea level. Let r u = {x u , y u , z u } ∈ R 3 and r b = {x b , y b , z b } ∈ R 3 be the coordinates of the quadrotor UAV's center of mass, (O u ), and that of the buoy, (O b ), in W, respectively. Let B u and B b be the body-fixed reference frames of the UAV at O u , and of the buoy at O b , respectively. The UAV is physically connected to the buoy by means of an inextensible massless cable of length l ∈ R >0 , and the connections to the cable are modeled as two frictionless revolute joints. The quadrotor UAV has a mass m u and a moment of inertia tensor J u = diag(J u,xx , J u,yy , J u,zz ) ∈ R 3 >0 in B u . Let the floating buoy have a volume b ∈ R >0 , a bounded mass m b ∈ (0, ρ w b ), with ρ w being the water density, and a moment of inertia tensor Also, let the orientation of B u and B b with respect to W be described by the Euler and Ω u = [p u , q u , r u ] ⊺ ∈ R 3 be the UAV's linear and angular velocities in B u , respectively; and be the buoy's linear and angular velocities in B b , respectively. Furthermore, let the translational velocity transformation matrix from any body frame, B • , to W be described as: and the rotational velocity transformation matrix from W to any body frame be described as: We also define the rotation matrix about the vertical z-axis only (yaw) from one body frame, B • , to W as: Both the UAV and the buoy are subject to cable tension, T ∈ R ≥0 , and to gravitational acceleration, g. The tensile force is applied on the tether hinge in each platform located at distances, r GH,u and r GH,b , for the UAV and the buoy, measured in their individual body frames, respectively. Moreover, the UAV's propulsion is described by the input vector u u = [u 1 , u 2 , u 3 , u 4 ] ⊺ ∈ R 4 including the thrust force vector ≥0 , and three torques that induce a rotational motion about each of the body frame axes such that τ u = [u 2 , u 3 , u 4 ] ⊺ ∈ R 3 . The dynamics of the UAV actuators (motor propellers) are neglected in the modeling given that their response time is considerably smaller than that of the UAV (> 10X). Finally, the buoy is subjected to hydrodynamic and hydrostatic forces that are subsequently described. The system is left-handed, with r ≥ 0, α ∈ (− π 2 , π 2 ], and ϕ ∈ (−π, π].

Spherical Coordinates
The spherical coordinates system for the tethered UAV−buoy system, W ′ , is illustrated in Fig. 2. The Cartesian position of the UAV in W with respect to W ′ is defined as: r rel = [x rel , y rel , z rel ] ⊺ = r u − r b ∈ R 3 , where the subscript (rel) denotes 'relative. ' We let its spherical coordinates in W ′ , r = [r, 0, 0] ⊺ , be represented by the triplet {r, α, ϕ} with unit vectors {ê r ,ê α ,ê ϕ }, such that the elevation angle α ∈ [−π/2, π/2] is defined between r and the horizontal plane, and the azimuth angle ϕ ∈ (−π, π] is defined between the projection of r on the {x, y} plane and the positive x-axis, which is mathematically represented as: To transform vectors from spherical to Cartesian coordinates, two consecutive rotations are needed such that: where R ϕ and R α are defined in C4 and C5, by which the Cartesian coordinates can be retrieved from the n wave component Figure 3 Water medium model elements and description, including: wave characteristics, wave velocity, floating object orientation, immersed height, Stokes-drift current, and lumped water surface current. spherical coordinates such that: which expands to: Conversely, the spherical coordinates can be retrieved from the Cartesian coordinates such that: where R C2S = R ⊺ S2C . Hence, the UAV's inertial coordinates can be obtained from the buoy's inertial coordinates and the relative spherical coordinates as follows: Similarly, the UAV's velocity and acceleration vectors are calculated via:ṙ whereṙ andr are explicitly given in C7.

Water Medium Model
Consider the water medium of an ocean/sea environment, as visually depicted in Fig. 3, with two main elements that are detailed next: the water surface current and gravity waves.

Gravity Wave Model
Assuming a large water depth compared to the wavelength of gravity waves, at time t and horizontal position {x, y}, the water elevation variation, ζ, due to gravity waves is described statistically as [5]: ζ(x, y, t) = N n A n sin − k n x cos ψ n,w − k n y sin ψ n,w + ω n t + σ n , where A n , k n , ω n ∈ R ≥0 , ψ n,w , and σ n ∈ (−π, π] are respectively the wave amplitude, wave number, circular frequency, wave propagation angle with respect to the inertial frame's x-axis, and random phase angle of wave component number n ∈ S n with S n = {1 ≤ n ≤ N | N ∈ N}, and N being the total number of active wave components. Furthermore, referring to the dispersion relation in deep water, the wave number is given as k n = ω 2 n /g. The wave-induced velocity vector of the fluid particles in W, U w = [u w , v w , w w ] ⊺ ∈ R 3 , has its components described as [5]: u w (x, y, z, t) = N n ω n A n e knz sin(−k n x cos ψ n,w − k n y sin ψ n,w + ω n t + σ n ) cos ψ n,w , v w (x, y, z, t) = N n ω n A n e knz sin(−k n x cos ψ n,w − k n y sin ψ n,w + ω n t + σ n ) sin ψ n,w , w w (x, y, z, t) = N n ω n A n e knz cos(−k n x cos ψ n,w − k n y sin ψ n,w + ω n t + σ n ), where ω n relates to the wave period, T n , via ω n = 2π/T n . Finally, the water surface plane at point {x, y, ζ} has tilt angles in roll and pitch with respect to the inertial frame, which are calculated by differentiating (11) with respect to x and y, respectively: φ w (x, y, t) = atan N n −A n k n sin ψ n,w cos(ω n t − k n x cos ψ n,w − k n y sin ψ n,w + σ n ) , θ w (x, y, t) = atan N n −A n k n cos ψ n,w cos(ω n t − k n x cos ψ n,w − k n y sin ψ n,w + σ n ) .

Water Current
The water surface current acts in the {x, y} (horizontal) plane of W, such that U c = [u c , v c , 0] ⊺ ∈ R 3 , and is given as: [6], with its components in the inertial frame, W, defined as: and is the resulting sum of other water current components, determined as follows: where U l and ψ l are the average velocity and direction of the current in W, respectively.

Buoy's Dynamic Model
Various types of forces affect the buoy's motion, with a major contribution from restoration, damping, and radiation forces. These forces substantially depend on the immersed volume of the buoy, im ∈ [0, b ], which is a function of the buoy's elevation, defined Let the vector from the buoy's center of gravity to its center of buoyancy in W be defined as Additionally, let the tether hinge location with respect to the buoy's center of gravity be , as seen in Fig. 1, where the cable tension is applied to the buoy in W ′ , such that

Assumption 1
The water-buoy friction dominates the airbuoy friction, thus we neglect the air drag on the buoy.
Consider the buoy dynamics in B b , and let its state Applying Newton's second law of motion gives: are respectively the buoy's inertia, damping, and Coriolis matrices ; the gravitational forces and moments vector are included in G ′ b ∈ R 6 ; and the external forces and moments are captured in τ ′ b ∈ R 6 . The inertia matrix expands as The buoy's total damping in B b is expressed as: where the generalized radiation-induced potential damping matrix is expanded as is the skin friction matrix, calculated as: where is the wetted area of the buoy, and D S,4−6 ∈ R ≥0 can be approximated by considering the moments effect of D S,1−3 over the buoy's surface. The effect of the wave drift damping matrix, D W ∈ R 3×3 , is already included in the Stokes drift velocity in (15), thus it will be dropped from (19). The buoy dynamics in (17) are expressed in W with the state vector, where M b , D b , and C b ∈ R 6×6 are the buoy's inertia, damping, and Coriolis matrices, respectively, are respectively the vectors of the gravitational and external forces and moments on the buoy in W given by: We also define:

UAV's Dynamic Model
We let the UAV's thrust vector in the Cartesian frame, be calculated as: with its elements being explicitly represented as: In the spherical frame, the UAV's thrust vector, f u,S = [u r , u α , u ϕ ] ⊺ , is expressed as: with the element u r , u α , and u ϕ being explicitly represented in C8. The tether's tension on the UAV expressed in W ′ , representing the distance from the UAV's center of gravity to its tether hinge. Finally, the local wind speed that disturbs the UAV's motion is defined in The quadrotor UAV system dynamics in W are obtained from Newton's second law of motion with the state vector η u = [η 1,u ; η 2,u ], where η 1,u = r u and η 2,u = Ξ u , yields: where: ≥0 are the UAV's translational and rotational damping friction matrices, respectively. The relative velocity vectors of the UAV areη 1,u 2,u in rotation. τ 1,u and τ 2,u ∈ R 3 are vectors of other external forces in W and moments in B u of the UAV, respectively, and they include the rotors' thrust and tether tension effects, expressed as: The damping matrix elements, D u,i , i ∈ {1, ..., 3}, are approximated as: is the UAV's cross-sectional area in the respective plane, ρ a is the air density andν u,i is the UAV-wind relative velocity in the respective body frame axis. The adopted model captures the major elements required to represent the quadrotor in a tethered UAV−buoy system with slow-to-moderate dynamics, thus acrobatic maneuvers and their influence on the system dynamics are not considered. For more details on the quadrotor UAV model, readers can refer to [30].

System Constraints
In this section, we generalise the system's 2D constraints presented in [14] to the 3D space. These constraints help establish the bounds and operating conditions for when the coupled model of the tethered UAV−buoy system is applicable.

Taut-Cable Constraint
Since the cable is assumed inextensible, it remains under tension (taut) at time t if r(t) = l, which is equivalent to: We label the tethered UAV−buoy system as 'coupled' when (30) is satisfied, and 'decoupled' otherwise.

No Buoy-Hanging Constraint
To keep the buoy floating on top of the water surface, the vertical tension component transmitted from the UAV through the tether must not exceed the weight of the buoy, otherwise it will be lifted into the air and the system will reduce to a UAV with a slung payload. This is prescribed by: as deduced from (21) and (22) at steady-state.

No 'Fly-Over' Constraint
'Fly-over' occurs when the buoy starts hopping over wave crests [8]. This phenomenon is avoided if: which means that the buoy remains partially immersed at all times. The 'fly-over' phenomenon is related to the total surface velocity of the buoy, described as: where V and ψ V are its magnitude and direction in the horizontal plane, respectively. Next, we define the wave encounter frequency for the n th wave component, ω e,n , as [5]: The excitation of the buoy's heave dynamics at ω e,n induces 'fly-over' if it approaches the heave's natural frequency, as interpreted in [14].

The Tethered UAV−Buoy System Model
In a coupled form, the tethered UAV−buoy system dynamics can be formulated by referring to the Euler-Lagrange formulation, which leverages the results in Sections 2.3, 2.5, 2.6, and 2.7. First, we define the Lagrangian function as L(q,q) = K(q,q) − U (q), where K(q,q) ∈ R ≥0 and U (q) ∈ R are the kinetic and potential energies of the system, with being the generalized coordinates vector. The equations of motion of the UAV−buoy system are obtained via: where τ ∈ R 11 is the external forces and moments vector. Per Assumption 1, let D be the global damping matrix formulated based on (23) without including a wind-induced component, and the dissipative forces are captured by the power function, P ∈ R, such that ∂P ∂q := Dq.q is defined as: If tether is assumed to have a negligible mass, the coupled system's kinetic energy is obtained as the sum of the individual energies of the UAV and the buoy as follows: where the global inertia matrix of the UAV−buoy system, M, is formulated by using the elements of M u and M b , as described in (C11). Next, the potential energy of the system can be formulated by referring to (22) and (2.6) as: The details of the Euler-Lagrange formulation in (35) are detailed in Appendix D, which finally leads to the following equations of motion: where C represents the global Coriolis matrix, and G is the global vector of gravity forces and moments.
Assumption 2 The design of a stable buoy is beyond the scope of this work, and only buoys with inherited stability are considered. Thus, we assume that the buoy's center of buoyancy always lies above its center of gravity, and that the roll and pitch dynamics of the buoy are damped and stable, which indicates that D b,44 > 0 and D b,55 > 0. Thus, we assume that the buoy remains tangent to the water surface.
With Assumption 2 and the dominance of waves with relatively long periods and moderate heights, the time derivatives of the buoy's roll and pitch angles, φ b andθ b , are small and thus their effects can be neglected inṀ.
Since the tether joints are considered revolute and frictionless, there is no coupling between the tether's rotational motion and the UAV and buoy's rotational dynamics, whereas their translational dynamics are coupled through (9). Additionally, as can be deduced by inspecting the elements of (35), the UAV and buoy's rotational dynamics are independent of each other, and can still be described by the dynamic models of their individual systems in (27) and (21), respectively. This is further elaborated in Appendix D. Given this, the tethered UAV−buoy system coupled dynamics can be represented by the first fives states: x b , y b , z b , α, and ϕ, which are independent of other rigid body orientation states that concern either the UAV's or the buoy's dynamics alone. Hereafter, we limit the representation of the coupled system to the first five state variables mentioned above, while the UAV and buoy's rotational dynamics can still be described by η 2,u in (27) The Coriolis matrix is explicitly represented as: and the damping matrix is explicitly represented as: Additionally, the gravitation force vector and external forces and torques vector are explicitly represented as: and respectively. If constraints (30) and (32) are satisfied, the coupled form of the dynamic model equations is given by: After solving the above differential equations, the UAV's position and velocity vectors can then be computed from (C9) and (C11), respectively. We note that if the tether's weight is to be considered, the tethered system's model should be updated as per [15].

Control System Design
The control system design problem is defined as manipulating the surge velocity of the buoy, u b , to track a desired reference, while orienting the cable in the desired elevation angle, α, and azimuth angle, ϕ. This allows for applying a tension force in any required direction, while tracking the desired surge velocity.
As an extension to the Surge Velocity Control System (SVCS) presented in [14], which was designed to only control u b and α, the herein controller will be dubbed as Directional Surge Velocity Control System (DSVCS), given that it adds the pulling force direction (azimuth angle, ϕ) to its controlled states. We note that solely controlling these states via setpoint tracking does not yield inertial velocity tracking.
For example, external forces that result in sway motion (in the direction of the buoy's body-fixed y-axis) are not rejected by the controller. That said, the controller with the reduced states can be equipped with a path planner, which can apply lateral tension by choosing a specific azimuth angle to counter the external forces in the sway direction. Such an architecture would give the system the ability to track reference trajectories in the inertial frame, but this is beyond the scope of this work, thus we stick with the design of the DSVCS.

Operational Modes and State Machine
To complement the surge velocity controller, a UAV−buoy relative position controller is required to control the radial distance of the UAV in W ′ , r, instead of u b . Both controllers are part of the DSVCS. In [16], a state machine was proposed to allow the system to switch between the coupled and decoupled states by alternating between the two controllers. In addition, the state machine can trigger a change in the positioning of the UAV with respect to the buoy between front and rear to change the pulling direction. Contrarily, in the 3D problem of this work, we only consider the forward motion of the buoy, while the velocity vector steering will be handled by changing the cable's azimuth angle. The purpose of the state machine is to allow the system to smoothly couple and decouple when requested by the controller, which will be sufficient to achieve the goals of the controller in this work. For that, we adopt the same state-machine of [16], with the exclusion of the 'repositioning' maneuver trigger part.

Controller Design
With the inner-loop/outer-loop cascaded structure shown in Fig. 4, the DSVCS allows tracking of the states u b , α, and ϕ by orienting and scaling the UAV's thrust vector. The setpoint for the controller is (ū b0 ,z u0 ,φ 0 ) in the buoy's surge velocity control mode, and (r 0 ,z u0 ,φ 0 ) in the UAV's relative position mode. A prepocessing unit then generates the resultinḡ α based onz u0 , and a smoothed signal for the tracking setpointsū b0 ,r 0 , andφ 0 . When the control mode changes, a fast and smooth transition function handles the switching between the two control laws [14]. The resulting force vector with its three components alongê r ,ê α , andê ϕ is transformed to the inertial frame to determine the desired UAV's thrust vector magnitude and orientation. At this point, the problem reverts to a basic UAV thrust decoupling problem, where the resulting thrust vector magnitude and tilt angles can be easily computed, as shown later, while the yaw angle remains free to set. Since the UAV is expected to maintain visual tracking of the buoy, the controller sets the desired yaw angle to equal the tether's azimuth angle.
In summary, and as shown in Fig. 5, the proposed DSVCS can control three main variables: 1) the relative position between the UAV and the buoy (outer-loop), 2) the buoy's surge velocity (outer-loop), and 3) the UAV's attitude (inner-loop). A state machine selects which of the two outer-loop controllers to activate based on the coupling between the UAV and buoy.

Reference Signals and Velocity Setpoint
The UAV's desired motion during a buoy manipulation task is limited to a horizontal plane of constant elevation (z u0 ), which reduces the UAV's power consumption [1] and results in safer and more predictable paths. This can be achieved by computing the resulting reference elevation angle, α, based on the desire reference elevation,z u , as: The azimuth angle is chosen to set the pulling direction, and the UAV's yaw angle can be independently manipulated without affecting the other states; here, we set it equal to azimuth angle in order to keep the UAV directed forward while an onboard camera can still point towards the buoy at all times: Note that the steady-state velocity vector of the buoy in the horizontal plane of W,ψ V , does not necessarily point in the same direction asφ, since the buoy is free to slide sideways, such that: where ǫ ψ is an error angle that is related to the sideslip angle of the buoy, β u . Finally, the radial position, r 0 , and the velocity setpoint,ū b0 , are smoothed by low-pass filters of fourth-and second-order, respectively, which results in a smoother performance due to respecting the system dynamics [25].

UAV−Buoy Relative Position Control Law (Outer-Loop)
To design the relative position control law, we must refer to the UAV dynamics in W ′ where the states of interest are explicitly expressed. If we reorder the realization ofr u in (10b) and multiply both sides by m u , we get: The first right-hand side term of (49) can be seen as the mapping of the translational terms of (27) into W ′ (kinetics), and the second one can be seen as the mapping of the buoy's linear accelerations into W ′ (kinematics). Expanding (49) yields: Consider the case of nonzero tension for the relative position dynamics of the UAV−buoy system's in (50), with states vectors X 1 = [r, α, ϕ] ⊺ and X 2 = [ṙ,α,φ] ⊺ , control input vector U = [u r , u α , u ϕ ] ⊺ , subject to unknown external disturbances including water currents, gravity waves, and wind gusts. The equations of motion in the kinetic form are expressed as: where the subscript (dc) refers to the decoupled dynamics,Ẍ dc =Ẍ 1 , the inertia matrix M dc = m u diag(1, r 2 , r 2 c α ), the parameter vector Θ dc = T , the regressor vector Φ dc = [−1; 0; 0], and b dc = diag(1, r, r). δ dc represents the vector of lumped modeling errors and disturbances, and the vector H dc = [H dc,1 ; H dc,2 ; H dc,3 ] represents all of the nonlinear Euler, Coriolis, centrifugal, and gravitational forces and moments, and is given by: Water medium effect Wind effect which results inV 2 = −e ⊺ 1 k 1 e 1 − e ⊺ 2 k 2 e 2 . Via Barbalat's lemma with Assumption 4, the asymptotic convergence of V 2 to zero is guaranteed. Note that if Assumption 4 is violated by the presence of strong wave disturbances or wind gusts, stability and finite tracking error are still achieved by increasing the controller gains, such that they overcome the disturbances mismatch effect onV 2 . Finally, the control law in the PID-like form in (53) is obtained by substituting e 2 andΥ in (54), then setting e I 1 :=δ(γk 1 ) −1 .

Buoy Surge Velocity Control Law (Outer-Loop)
If the cable is taut, its length becomes constant andṙ andr can be set to zero. Substituting l for r in (50) yields the UAV's motion equations in W ′ in spherical coordinates notation: In order to control the buoy's surge velocity, u b , it must be explicitly expressed in (55). The buoy's acceleration in the body frame is described as:ν 1,b = R ⊺ 1,bη 1,b . When waves with relatively long periods and moderate heights are dominant, the resulting buoy's roll and pitch angles are small and cannot be used by the controller unless a proper measurement technique is available for the UAV. On the other hand, the USV/buoy heading can be visually estimated with good accuracy as was demonstrated in [4]. For this reason, we limit the transformation of R 1,b to yaw only, such thatν 1,b ≈ R ⊺ z bη 1,b . Hence,ẍ b andÿ b in (55) can be transformed intou b andv b accordingly.
Consider the UAV dynamics in W ′ in the coupled case, while following the spherical coordinates notation as presented in (55). By choosing the state vector X ′ 1 = [u b , α, ϕ] ⊺ , we can rewrite (55) in the kinetic form as: where the subscript (cp) refers to the coupled dynamics,Ẍ cp =Ẍ ′ 1 , the control vector U ′ = U , the parameter vector Θ cp = T , the regressor vector Φ cp = [−1; 0; 0], and b cp = diag(1, r, r). δ cp represents the vector of lumped modeling errors and disturbances, the inertia matrix is expressed as: The state-space model in (56) is formulated as a timevarying second-order nonlinear system as: The outer-loop velocity controller manipulates the state vector X ′ 1 , and a control law for the surge velocity can be designed in a similar fashion as described in Section 3.2.2, with a difference that only one step is required in the backstepping process for the state u b . LetX ′ 1 = [ū b ,ᾱ,φ] ⊺ be the reference state vector and e ′ 1 = X ′ 1 −X ′ 1 be the respective state vector error. The control law is defined as: >0 , are controller gains that are defined next; andΘ ′ is the estimate of Θ ′ . Theorem 2 Consider the state-space representation in (57) of the UAV−buoy's relative position dynamics in (55) when the system is coupled. If Assumption 3 holds true, the control law in (58) generates the thrust vector elements in the spherical reference frame ur, uα, and uϕ, that can stabilize the outer-loop dynamics of the system, and reduce the tracking error to a small region neighboring the origin in finite time for a set of gains k ′ 1 , k ′ 2 , and γ ′ ∈ R 3×3 ≥0 , defined in a similar way to their counterparts in Theorem 1. Additionally, if Assumption 4 holds, the tracking error is reduced to zero in finite time.
Proof We employ the backstepping control design, which involves one step for the u b state, and two steps for the α and ϕ states. To keep a compact representation of the proof in vector form, the buoy velocity states are tackled in the second step. First, the Lyapunov function V ′ 1 = 1 2 e 2 α + 1 2 e 2 ϕ is proposed, and let its derivative be expressed asV ′ 1 = eαėα + eϕėϕ. We proceed to a second step in the control design process for the states' rates given thatėα andėϕ do not include an explicit control input. To stabilize e 1 , we define a virtual control input as: Next, we define the rates error as: e ′ 2 =Ẋ ′ 1 − Υ ′ . Note that the rates error of the buoy velocity's x-component simply reduces to e ′ 2,1 = u b −ū b . Here, we define a second Lyapunov function as: withδ ′ =δ ′ − δ ′ , then we differentiate to have: Finally, the control inputs vector and the update rates of the lumped disturbances and modeling errors are chosen to satisfy the negative semi-definiteness ofV ′ 2 , as: where the tuning gain k ′ 2 ∈ R 3×3 is a diagonal matrix, and we getV ′

UAV's Attitude Controller (Inner-Loop)
Let φ u,c and θ u,c be the desired roll and pitch angles of the UAV, respectively, which can be obtained from the outputs of the outer-loop controller along with the UAV's total thrust command, u 1,c . The transformation from thrust components in spherical coordinates to the UAV's thrust and tilt angles is performed by referring to (25) and (26) in two steps. We let the command thrust vector in the world frame, W ′ , be U C,c = [u x,c , u y,c , u z,c ] ⊺ = R u,1 [0, 0, u 1,c ] ⊺ , and in the spherical frame, W ′ , be U S,c = [u r,c , u α,c , u ϕ,c ] ⊺ = R C2S U C,c . In the first step, given the force commands in the spherical frame generated from the outer-loop controller, U S,c , we calculate U C,c = R S2C U S,c , which expands to: In the second step, we calculate φ u,c and θ u,c given that U C,c = R u,1 [0, 0, u 1,c ] ⊺ , which expands to: With mathematical manipulation, we get the following relationships: x,c + u 2 y,c + u 2 z,c , φ u,c = arctan(u x,c s ψu − u y,c c ψu )/u z,c , θ u,c = arcsin(u x,c c ψu + u y,c s ψu )/u 1,c .

(62)
Let φ ′ u,c = φ u,m tanh φ u,c /φ u,m and θ ′ u,c = θ u,m tanh θ u,c /θ u,m be smooth and bounded versions of φ u,c and θ u,c , where φ u,m and θ u,m ∈ (0, π 2 ) are the absolute upper limits of the UAV's roll and pitch angles, respectively. The UAV's yaw dynamics can be be controlled independently of the buoy's manipulation.
The state-space model in (63) is formulated as the following time-varying second-order nonlinear system: LetX ′′ 1 = [φ u ,θ u ,ψ u ] ⊺ be the reference state vector, and e ′′ 1 = X 1 −X ′′ 1 be the respective state error vector. The proposed control law that stabilizes the UAV's rotational dynamics is defined as: where k ′′ P , k ′′ I , and k ′′ D ∈ R 3×3 >0 are controller gains, andΘ ′′ is the estimate of Θ ′′ . Equivalent outcomes of Theorem 1 are applicable for the UAV's inner-loop controller, and can be proved by replicating the procedure of Theorem 1's proof, which is omitted for brevity.

Parameter Estimation
If the taut-cable condition in (30) holds, the following realization of the cable tension is given: where ǫ α ∈ R ≥0 is a constant that prevents singularity in a small region near α = π 2 . The cable tension expression in (66) can be determined from the sum of the first and second rows of the buoy dynamics in (21), which shows a direct link with u b for the first case. The singularity near the vertical cable configuration (α = π/2) is handled by the UAV dynamics in (27) to compute the actual cable tension, T , as in the second case of (66). To get an estimateT of T , as required by the control laws of the DSVCS, (66) should be used. Here, we note that even ifT is inaccurate, the controller is able to compensate for the tension effect based on its integral term action.

Simulations
This section presents numerical simulations for validating the 3D tethered UAV−buoy robotic system model, and for evaluating the performance of the proposed control system in allowing the UAV to manipulate the buoy in the horizontal plane as well as the tether's elevation and azimuth angles. We validate the system in realistic environmental and operating conditions including wind, waves, and water current. To increase the model's fidelity, the UAV model includes the motor and propeller dynamics, and the controller uses non-exact feedback signals (e.g. noise and bias), which simulate measured outputs based on state-of-the-art navigation sensors [14].

Simulation Settings
The validation of the proposed system is carried in the MATLAB Simulink ® simulation environment. The system includes a medium-sized quadrotor UAV tethered to a small buoy having a boat-like shape to enable self-alignment along the pulling direction. The parameters of the system are included in Table 1. The UAV's thrust-to-weight ratio is considered about 2.5, which gives a maximum thrust of about 120 N. The motor's model is considered a low-pass filter of the firstorder with a time constant τ m = 0.05 s. The buoy's immersed volume, wetted area, and skin friction coefficients are calculated as per [14]. Additionally, the added mass and damping are calculated based on the strip theory for surface vessels at low oscillating frequencies with their values presented in Table 1 [6].

Simulation Scenarios
The designed controller's performance and the fidelity of the derived system model are validated with different simulation scenarios. We first simulate two cases that include constant water current, wind gust, and moderate waves. To further validate the system in additional conditions, we offer extended simulation scenarios (in Section 4.4) by varying the water current direction, the wave components in the water environment, and other system parameters such as buoy size and cable length.
The first two cases, C1 and C2, both include a wind gust of U wd = [−5, 0, 0] m s −1 and a water current component with U l = 1 m s −1 and ψ l = π. The two scenarios include: • C1: water current and wind gust only.
• C2: water current, wind gust, and moderate waves with one wave components (N = 1), such that: A 1 = 0.75 m, ψ 1,w = 0, T 1 = 5.7 s, and σ 1 = π. In both cases, the UAV is initially hovering around the buoy, it then positions itself in a location that is suitable to start pulling the buoy at mean elevation,z u , of 5.0 m, which corresponds to a mean elevation angle ofᾱ 0 = 45 • . In the first stage, the UAV pulls the buoy until it reaches a surge velocity ofū b = 5 m s −1 . In the second stage, the UAV gradually steers the buoy by commanding an azimuth angleφ that increases from 0 to 360°, thus completing a full rotation. The resulting profile is shown in Fig. 7.
The buoy has an initial velocity that matches its surrounding water as calculated via (12) and (14); by referring to Assumption 2, we assume that buoy stays tangent to the water surface, such that: while its yaw dynamics are governed by (17).

Simulation Results and Discussion
The simulation results for C1 and C2 are shown in Fig. 6a (left) and Fig. 6b (right), respectively. Stable and accurate performance is noticed in both cases, where the UAV is able to successfully drive the buoy at the desired velocity (u b ), while maintaining good tracking of the reference elevation (α) and azimuth (ϕ) angles. As seen in the (r) subplot, the UAV first positions itself near the pulling location, it then enters the pulling state by making the cable taut during the entire duration of the buoy's manipulation between [12,100] s. Since there are no ambient waves in C1, there is no apparent fluctuation in the buoy's vertical position (z b ), immersed volume ( b ), and elevation angle (α); contrarily, the waves in C2 induce fluctuations in these three states. It is noted that the UAV keeps a level flight in both cases, as seen in the (z u ) subplot in Fig. 6a-b. The UAV is commanded to remain oriented forward during the entire maneuver, which induces large pitch angles (θ u ) to counteract disturbances, the roll angle (φ u ) remains near-constant, and the yaw angle (ψ u ) smoothly follows the azimuth angle of the cable (ϕ) in order to complete a full rotation during the course of the motion trajectory. The full rotation in azimuth (ϕ) causes the buoy to encounter the waves in different directions, where the frequency of encounter is low during the initial phase ([20, 40] s) due to following seas, and when the buoy turns around ([50, 60] s) the encounter frequency becomes high due to head seas. The cable tension, seen in the (T ) subplot, increases with velocity when the frequency of encounter with waves is high, and when the rate of change of the azimuth angle is high. Finally, we note that the thrust and torques inputs (u (1,2,3,4)c ) are bounded and stable. Even though the DSVCS achieves good tracking of the desired surge and azimuth states, it does not necessarily translate into identical trajectories in the world frame. This is demonstrated in Fig. 7, which shows the planar trajectory of the buoy for cases C1 and C2. It is seen that the two trajectories do not coincide, even though the results of Fig. 6 showed good tracking with no steady-state errors. It is it also noticed that the final direction of the trajectory is not exactly at zero as seen in the 'xy' subplot of Fig. 7, even though the final azimuth angle is zero. This can be understood by observing the non-zero sway velocity, v b , in the 'Vel.' subplot of Fig. 7. We also notice that the sway velocity of the buoy, v b , can be as large as its surge velocity, u b , which results in a larger absolute surface velocity of the buoy, V . Finally, we note that the buoy is prone to having a nonzero side-slip angle, β b , as seen in Fig. 7, given that it does not have any constraint on its planar dynamics. At the start of the simulation (< 17s) when the cable is slack, the side-slip angle points to a backward motion (β b = 180°) since u b ≈ −1 m s −1 and v b ≈ 0 m s −1 due to the initial water current velocity.
As mentioned in Section 3, the above demonstration motivates the need for a higher level path planner in the event that the buoy is required to track trajectories in the inertial frame. That said, the obtained simulation results demonstrate the effectiveness of the DSVCS in controlling the buoy's surge velocity, u b , and the direction of the pulling force by manipulating the azimuth (ϕ) and elevation (α) angles. This provides the system with the ability to move forward and steer as a locomotive, which also lays the foundation for designing a path planner in the inertial frame.

Extended Simulations
To further validate the system performance in different conditions, we offer additional simulation scenarios by varying the water current direction, the wave components in the water environment, and well as changing the system parameters (buoy size and tether length). The physical system index is referred to as Si where i = {1, 2}, and the environment index is  referred to as Ej where j = {1, ..., 7}. Hence, we introduce the following system variations: (a) System size variants (buoy and cable): S1. Same system presented in Section 4, except for the cable length: l = 10 m (S1), and m b = 16 kg. S2. Same system as (S1) except for cable length:   8 shows the resulting trajectories of the simulation in environments E1 through E7, where system S1 is adopted with E1 through E5, and S2 is adopted with E6 and E7. From E1, we notice that without wave and current disturbances, the system can follow a circular trajectory by simply commanding a constant velocity and azimuth angle with a constant ramp. If water current exists, the buoy drifts towards the direction of the current as seen in E2, E4, and E6. When the system is deployed in wavy seas, additional drift may occur due to the waves' Stokes-drift effect, which increases with wave height, as seen by comparing the additional drift between E3 (lowest deviation with respect to the current-only case), E5, and E7 (highest deviation with respect to the current-only case), whereas the oscillating component of the waves' velocity tends to cancel out near the end. Finally, we see that the system is able to deal with relatively large waves and reject their respective disturbances, as observed in E7.
In all cases, tracking of the controlled states is achieved with small steady-state errors, which is practical for applications that do not require precision motion control. One room for improvement can be achieving better trajectory tracking in the world frame, which is sub-optimal due to the fact that the DSVCS has no direct effect on the buoy's sway velocity, and thus it cannot trace repetitive trajectories in the world frame. Again, if tracking trajectories in the world frame is desired for certain applications, that can be achieved by employing a high level path planner by means of commanding appropriate combinations of u b and ϕ, which is not within the scope of this current work.

Conclusion
This work extended the 2D marine locomotive UAV problem, realized as a tethered quadrotor UAV−buoy system, to the 3D space to enable the UAV to steer the buoy in the horizontal plane, and not being restricted to forward and backward motion only. The dynamics of the buoy and the UAV were given as 6-DOF rigid bodies with a proper definition of all effective forces on the system. A detailed derivation of the coupled system dynamics was included to serve as a reference for the rigorous mathematical formulation and modeling of the heterogeneous multi-physics robotic system under consideration, which features high order, nonlinear, and coupled dynamics. Finally, the control problem for applying a tensile force via the tether with specific magnitude and direction was solved through the proposed DSVCS, which was realized by controlling the buoy's surge velocity and the tether's orientation.
Future work naturally flows into designing a path planner to control the buoy's location in the world frame, and study the effect of the tether and the buoy's tilt dynamics on the system.
Similarly, its damping matrix in W, D b , is symmetric (i.e., D ij = D ji , for i = j), and has its non-zero elements explicitly given as: