**1. Introduction**

In recent years, systems for the transportation of suspended-payload using unmanned aerial vehicles (UAVs) have attracted research interest. Some important applications are described in [1]. Quadrotor vehicles exhibit complex dynamic behavior, and if a cable-suspended payload is added to a quadrotor, it increases the complexity of the system because additional degrees of freedom (DOFs) due to the payload oscillation are introduced. Moreover, if uncontrolled, a cable-suspended payload changes the dynamics of the flying vehicle and it can result in an unstable system.

Recently, some control methodologies have been developed to attenuate the payload swing and to solve the general problem of a quadrotor carrying a cable-suspended payload. For example, a control algorithm based on a backstepping strategy is obtained in [2], that ensures trajectory tracking of the quadrotor regardless of the payload swing. However, attitude control of the UAV is not considered. Additionally, a control design for a two-dimensional quadrotor with a cable-suspended payload that enables tracking of the vehicle rotation, the payload rotation, or the payload position is presented in [3], and it was extended to the three-dimensional case in [4].

In [5], the authors develop a nested saturation controller capable of driving the vehicle to a specified position while simultaneously limiting the payload dynamic effect. For this work, Nicotra et al. considered the design for the two-dimensional case only and the attitude control of the vehicle was not considered. Moreover, a feed-forward control algorithm for reducing or canceling the payload's oscillation is introduced in [6]. This controller was designed by implementing the input shaping theory. In [7], the authors present a geometric controller to exponentially stabilize the aerial robot position while aligning the links vertically below the vehicle. Similarly, a tracking control law for a UAV with a load attached by a cable represented as successively-attached inflexible links was designed in [8]. A fixed-gain geometric PD control strategy is developed to reach the desired goal for a nominal load. An adaptive control law for an aerial robot carrying a payload attached using a cable was presented in [9]. In [10], an algorithm for parameterizing aerial vehicles transporting a payload employing a complementary constraint is presented. A nonlinear attitude controller is developed in [11] to stabilize the altitude, and the translational dynamics control law is introduced by converting the vehicle velocity and position error into rotation control. Rego and Raffo [12] address trajectory tracking for a two-dimensional aerial robot transporting a payload. A discrete-time mixed *H*2/*H* ∞ linear control strategy is presented. In [13], an active-model-based linear controller is designed for a UAV transporting a payload. A linear model is obtained considering the vehicle in hover flight mode. In [14], a path tracking controller is developed based on existing Lyapunov-based path tracking control laws for free-flying aerial vehicles, which are further backstepped through the vehicle rotation dynamics.

A passivity-based control technique is used in [15] to control the UAV such that cable-suspended payload swing reduction is achieved for the planar case. Here, the attitude control of the vehicle is not considered. Also, interconnection and damping assignment passivity-based control (IDA-PBC) without total energy-shaping for a UAV transporting a cable-suspended payload for planar maneuvers is developed in [16,17]. Two control laws with total energy-shaping are presented in [18], where the closed-loop inertia matrices are modified. These works compute PDEs for synthesizing the control law. For this reason, the control strategy is only designed within the longitudinal plane. The control design for the three-dimensional case yields complex partial differential equations (PDEs).

In the literature, unmanned aerial vehicles have been controlled using energy-based controllers. In [19], the design of two nonlinear controllers to stabilize an aerial vehicle characterized with quaternions and their axis-angle depiction is studied. Also, [20] introduces a PBC for a vertical take-off and landing (VToL) vehicle. An estimator of unmodelled dynamics and external wrench acting on the UAV and based on the momentum of the system is used to compensate such disturbance effects. Moreover, [21] develops an IDA-PBC methodology that is able to change the apparent vehicle dynamical parameters, while [22] proposes a robust control of underactuated aerial manipulators via IDA-PBC.

On the other hand, some works apply control strategies based on linear matrix inequalities to UAVs. In [23], a method for using LMIs to synthesize controller gains for a UAV system is presented. In [24], a nonlinear state feedback controller based on LMIs, and a technique with pole placement constraint (PDC) is synthesized. The requirements of stability and pole placement in the linear matrix inequality (LMI) region are formulated based on the Lyapunov direct method.

In this work, the control approach is based on a hierarchical scheme considering the well-known time-scale separation between rotational and translational dynamics of the quadrotor. On one hand, the objective of this paper was to design an energy-based control law for the outer-loop (i.e., for the underactuated dynamics of the system). This control law is proposed for the three-dimensional case, and it is based on the translation dynamics, which is able to lead the vehicle to a desired position while simultaneously reducing the payload swing. Compared with similar works that present control laws based on passivity and energy, particularly for underactuated systems, the proposed controller avoids solving complex PDEs to obtain the control law. On the other hand, a feedback controller based on an LMI for the inner loop which is fully actuated is presented. The controller based on an LMI for the rotational dynamics results in a control algorithm with relative simplicity and with an easy analysis to demonstrate its stability.

The contribution of this paper is the synthesis of a new controller for a class of Lipschitz nonlinear systems. An important limitation of the classic results for Lipschitz nonlinear systems is that they perform well only for small values of the Lipschitz constant. In the case when the Lipschitz constant is large, most of the existing controller design approaches fail to contribute a solution to the LMI. This article introduces an algorithm that operates for larger Lipschitz constants compared with classical results in the literature.

This paper is organized as follows. Section 2 describes the dynamical model for a three-dimensional aerial vehicle carrying a payload. Section 3 presents an approach for LMI-based Lipschitz nonlinear systems. This theoretical approach is applied to stabilize the quadrotor rotational dynamics. Section 4 proposes an energy-based control to stabilize the vehicle translational dynamics and to attenuate the payload swing. Section 5 presents numerical simulations and results. Finally, Section 6 gives conclusions and perspectives.

#### **2. Dynamic Model**

In this section, we present the mathematical model of a quadrotor transporting a payload connected by a cable. The aim is to present a dynamic model that mathematically describes the relationship between the quadrotor, the cable, and the payload. For this purpose, consider a rigid body with mass *m*, being transported by a quadrotor as shown in Figure 1. Note that in addition to the six DOFs of the UAV, the payload adds another two DOFs, resulting in a system with eight DOFs and four inputs.

**Figure 1.** Three-dimensional quadrotor with a cable-suspended payload.

The following assumptions were made for modeling the quadrotor with a cable-suspended payload:


As shown in Figure 1, the body-fixed frame is defined by *B* = {*<sup>e</sup>*1,*e*2,*e*3} and the inertial frame by *O* = *ex*,*ey*,*ez* . The location of the mass center of the vehicle relative to *O* is represented by *ξ* = *xyz T*, the attitude of the quadrotor is denoted by *η* = *ψθφ T* (i.e., yaw, pitch, and roll, respectively). *α* defines the payload swing angle in space and has two components *αx* and *<sup>α</sup>y*, *μ* = *αx <sup>α</sup>y T*. *αx* is the swing angle projected on the *XZ* plane and *<sup>α</sup>y* is the swing angle projected on the *YZ* plane, *β* represents the angle of the *X* axis and the projected line of the cable to the *Y* plane. Thus, the state vector is denoted by *q* = *xyz ψθφαx <sup>α</sup>y T* ∈ R8. The control input is represented by *u* = *f τ T* ∈ R4, where *f* = *f*1 + *f*2 + *f*3 + *f*4 defines the total thrust magnitude and *τ* = *τψ τθ τφ T* denotes the input torques.

The mass of the quadrotor and the payload are defined by *M* and *m*, respectively. The length of the cable is represented by *l*, the gravitational acceleration by *g*. Finally, the distance between the motors and the gravity center is equal to *d*.

In this work, the system is modeled using the Euler–Lagrange formulation.

#### *2.1. Euler–Lagrange Formulation*

Let *ξ p* = *xp yp zp T* ∈ R<sup>3</sup> be the payload position in the inertial frame. Thus, the quadrotor and payload positions are related by

$$
\xi\_p = \xi + lp\_r
$$

where

$$p = [\sin(\mathfrak{a}\_x)\cos(\mathfrak{a}\_y) \quad \sin(\mathfrak{a}\_y) \quad -\cos(\mathfrak{a}\_x)\cos(\mathfrak{a}\_y) \ ]^T.$$

The equations of motion are obtained using the Euler–Lagrange method. The kinetic energy of the payload-quadrotor system is given by

$$\mathcal{K}\_{Q-P} = \frac{1}{2} M^{\dot{\zeta}^T \dot{\zeta}} \dot{\zeta} + \frac{1}{2} m \dot{\zeta}^T\_{\mathcal{P}} \dot{\zeta}\_{\mathcal{P}} + \frac{1}{2} \dot{\eta}^T I \dot{\eta},\tag{1}$$

where *J* = *diag*[*Ixx*, *Iyy*, *Izz*]<sup>3</sup>×<sup>3</sup> is a symmetric positive definite constant inertia matrix of the quadrotor with respect to *B*.

The total potential energy is defined by

$$V\_{\mathbb{Q}-P} = (M+m)\lg z - mgl\cos(\boldsymbol{a}\_x)\cos(\boldsymbol{a}\_y). \tag{2}$$

From (1) and (2), the Lagrangian is given by

$$L = \frac{1}{2} M \mathfrak{J}^T \mathfrak{J} + \frac{1}{2} m \mathfrak{J}\_p^T \mathfrak{J}\_p + \frac{1}{2} \dot{\eta}^T l \dot{\eta} - (M + m) \mathfrak{J} z + mgl \cos(a\_x) \cos(a\_y). \tag{3}$$

Using the Lagrangian and the derived formula for the equations of motion:

$$
\mathcal{M}(q)\ddot{q} + \mathbb{C}(q,\dot{q})\dot{q} + \mathbb{G}(q) = \ddot{B}u,\tag{4}
$$

where M(*q*) ∈ R8×<sup>8</sup> denotes the inertia matrix, which is symmetric and positive definite, *<sup>C</sup>*(*q*, *q*˙) ∈ R8×<sup>8</sup> the Coriolis and centrifugal matrix, *<sup>G</sup>*(*q*) ∈ R<sup>8</sup> the gravitational vector, and the matrix *B*¯ ∈ R8×<sup>4</sup> is determined by the manner in which the control *u* ∈ R<sup>4</sup> is the input of the system. These matrices are defined as

$$
\mathcal{M}(q) = \begin{bmatrix}
 m+M & 0 & 0 & 0 & 0 & 0 & mlc\_{a\_x}c\_{a\_y} & -mls\_{a\_x}s\_{ay} \\
 0 & m+M & 0 & 0 & 0 & 0 & 0 & mlc\_{a\_y} \\
 0 & 0 & m+M & 0 & 0 & 0 & mls\_{a\_x}c\_{a\_y} & mlc\_{a\_x}s\_{ay} \\
 0 & 0 & 0 & l\_{xx} & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & l\_{yy} & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & l\_{zz} & 0 & 0 \\
 mlc\_{a\_x}c\_{a\_y} & 0 & mls\_{a\_x}c\_{a\_y} & 0 & 0 & 0 & ml^2c\_{a\_y}^2 & 0 \\
 -mls\_{a\_x}s\_{ay} & mlc\_{a\_y} & mlc\_{a\_x}s\_{ay} & 0 & 0 & 0 & 0 & -ml^2
\end{bmatrix},
$$

$$\mathbf{C}(q) = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & (I\_{yy} - I\_{zz})\theta \\ 0 & 0 & 0 & 0 & 0 & I\_r\Omega\_r + (I\_{zz} - I\_{xx})\dot{\psi} \\ 0 & 0 & 0 & 0 & (I\_{xx} - I\_{yy})\dot{\psi} - f\_r\Omega\_r & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ \end{bmatrix}$$

$$\begin{aligned} -ml(s\_{ax}c\_{ay}\dot{\kappa}\_x + c\_{ax}s\_{ay}\dot{\upsilon}\_y) & -ml(s\_{ax}c\_{ay}\dot{\upsilon}\_x + c\_{ax}s\_{ay}\dot{\upsilon}\_y) \\ 0 & -mls\_{ay}\dot{\upsilon}\_y \\ ml(-s\_{a}s\_{a}\dot{\upsilon}\_x + c\_{a}c\_{a}\dot{\upsilon}\_y) & ml(-s\_{a}s\_{a}\dot{\upsilon}\_x + c\_{a}c\_{a}\dot{\upsilon}\_y) \end{bmatrix} \end{aligned}$$

$$\begin{array}{cccc} ml(-s\_{a\_x}s\_{a\_y}\dot{a}\_y + c\_{a\_x}c\_{a\_y}\dot{a}\_x) & ml(-s\_{a\_x}s\_{a\_y}\dot{a}\_x + c\_{a\_y}c\_{a\_x}\dot{a}\_x) \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ -ml^2s\_{a\_y}c\_{a\_y}\dot{a}\_y & -ml^2s\_{a\_y}c\_{a\_y}\dot{a}\_x \\ ml^2s\_{a\_y}c\_{a\_y}\dot{a}\_x & 0 \end{array}$$

⎤

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

$$\begin{aligned} G(q) &= \begin{bmatrix} 0 & 0 & (M+m)\mathbf{g} & 0 & 0 & 0 & mgl s\_{\mathbf{a}\_x} c\_{\mathbf{a}\_y} & mgl c\_{\mathbf{a}\_x} s\_{\mathbf{a}\_y} \end{bmatrix}^T \\\\ B(q) &= \begin{bmatrix} s\_{\Phi}s\_{\Psi} + c\_{\Phi}c\_{\Psi}s\_{\theta} & c\_{\Phi}s\_{\theta}s\_{\Psi} - c\_{\Psi}s\_{\Phi} & c\_{\theta}c\_{\Phi} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & l & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & l & 0 & 0 \end{bmatrix}^T \end{aligned}$$

where *sθ* = sin(*θ*), *cθ* = cos(*θ*).

We can see in the above expressions that the cable-suspended payload affects the translational motion, but the rotational motion is not affected. Then, the vehicle attitude dynamics is decoupled from the payload–quadrotor translational dynamics. Thus, the quadrotor with a cable-suspended payload model can be divided into payload–quadrotor translational and quadrotor rotational dynamics. The following subsections show the translational and rotational dynamics.

#### *2.2. Translational Dynamics*

From (1) the model considering the translational motion is

$$
\mathcal{M}\left(\vec{q}\right)\ddot{\vec{q}} + \vec{\mathcal{C}}\left(\vec{q},\dot{\vec{q}}\right)\dot{\vec{q}} + \vec{\mathcal{G}}\left(\vec{q}\right) = \vec{B}\vec{u},\tag{5}
$$

$$\text{where } \overline{\eta} = \begin{bmatrix} \int \xi & \mu \end{bmatrix}^T, \overline{\mu} = \begin{bmatrix} F\_x & F\_y & F\_z \end{bmatrix}^T = \begin{bmatrix} f\left(\mathbf{s}\_\theta \mathbf{s}\_\theta + \mathbf{c}\_\theta \mathbf{c}\_\theta \mathbf{s}\_\theta\right) & f\left(\mathbf{c}\_\theta \mathbf{s}\_\theta \mathbf{s}\_\theta - \mathbf{c}\_\theta \mathbf{s}\_\theta\right) & f\mathbf{c}\_\theta \mathbf{c}\_\phi \end{bmatrix}^T,$$

the matrices are given by

$$
\tilde{\mathcal{M}}(\tilde{q}) \quad = \begin{bmatrix}
 m+M & 0 & 0 & ml c\_{\mathfrak{a}\_{x}} c\_{\mathfrak{a}\_{y}} & -ml s\_{\mathfrak{a}\_{x}} s\_{\mathfrak{a}\_{y}} \\
 0 & m+M & 0 & 0 & ml c\_{\mathfrak{a}\_{y}} \\
 0 & 0 & m+M & ml s\_{\mathfrak{a}\_{x}} c\_{\mathfrak{a}\_{y}} & ml c\_{\mathfrak{a}\_{x}} s\_{\mathfrak{a}\_{y}} \\
 ml c\_{\mathfrak{a}\_{x}} c\_{\mathfrak{a}\_{y}} & 0 & ml s\_{\mathfrak{a}\_{x}} c\_{\mathfrak{a}\_{y}} & ml^{2} c\_{\mathfrak{a}\_{y}} & 0 \\
 -ml s\_{\mathfrak{a}\_{x}} s\_{\mathfrak{a}\_{y}} & ml c\_{\mathfrak{a}\_{x}} s\_{\mathfrak{a}\_{y}} & 0 & ml^{2}
\end{bmatrix} \tag{6}
$$

*C* ˜(*q*˜, ˙*q*˜) = ⎡ ⎢⎢⎢⎢⎢⎣ 000 <sup>−</sup>*ml*(*<sup>s</sup>αx <sup>c</sup>αyα*˙*<sup>x</sup>* + *<sup>c</sup>αx <sup>s</sup>αyα*˙*y*) <sup>−</sup>*ml*(*<sup>s</sup>αx <sup>c</sup>αyα*˙*<sup>y</sup>* + *<sup>c</sup>αx <sup>s</sup>αyα*˙*x*) 000 0 <sup>−</sup>*mls<sup>α</sup>yα*˙*<sup>y</sup>* 000 *ml*(−*sαx <sup>s</sup>αyα*˙*<sup>y</sup>* + *<sup>c</sup>αx <sup>c</sup>αyα*˙*x*) *ml*(−*sαx <sup>s</sup>αyα*˙*<sup>x</sup>* + *<sup>c</sup>αy <sup>c</sup>αxα*˙*y*) 000 <sup>−</sup>*ml*2*sαy <sup>c</sup>αyα*˙*<sup>y</sup>* <sup>−</sup>*ml*2*sαy <sup>c</sup>αyα*˙*<sup>x</sup>* 000 *ml*2*sαy <sup>c</sup>αyα*˙*<sup>x</sup>* 0 ⎤ ⎥⎥⎥⎥⎥⎦ (7)

$$\begin{array}{rcl} \mathcal{G}(\vec{q}) &=& \begin{bmatrix} 0 & 0 & (M+m)\mathbf{g} & m\mathbf{g}\mathbf{l}\mathbf{s}\_{\mathfrak{a}\_{\mathfrak{a}}}\mathbf{c}\_{\mathfrak{a}\_{\mathfrak{g}}} & m\mathbf{g}\mathbf{l}\mathbf{c}\_{\mathfrak{a}\_{\mathfrak{a}}}\mathbf{s}\_{\mathfrak{a}\_{\mathfrak{g}}} \end{bmatrix}^{T} \end{array} \tag{8}$$

$$
\mathcal{B}\overline{u} = \begin{bmatrix} F\_{\overline{x}} & F\_{\overline{y}} & F\_{\overline{z}} & 0 & 0 \end{bmatrix}^T \tag{9}
$$

,

where *sθ* = sin(*θ*), *cθ* = cos(*θ*).

#### *2.3. Rotational Dynamics*

From (1) the model considering the rotational motion is

$$\mathcal{M}^\*\left(\eta\right)\ddot{\eta} + \mathcal{C}^\*\left(\eta,\dot{\eta}\right)\dot{\eta} = B^\*\tau\_\prime$$

where

$$\mathcal{M}^\*(\eta) = \text{diag}\left[I\_{\text{xx}}, I\_{yy}, I\_{zz}\right]\_{3\times 3'} \mathcal{B}^\* = I\_{3\times 3}$$

$$\mathcal{C}^\*(\eta) = \begin{bmatrix} 0 & 0 & (I\_{yy} - I\_{zz})\theta \\ 0 & 0 & I\_r \Omega\_r + (I\_{zz} - I\_{xx})\dot{\psi} \\ 0 & (I\_{xx} - I\_{yy})\psi - I\_r \Omega\_r & 0 \end{bmatrix}$$

where *Jr* is the rotational moment of inertia.

#### **3. LMI-Based Approach for Lipschitz Nonlinear Systems**

In this section, we propose a nonlinear state feedback controller based on a linear matrix inequality to stabilize the quadrotor rotational dynamics.

Taking the state vector of the attitude dynamics as *δ*(*t*) = *ψ ψ θ* ˙ ˙*θ φ φ*˙ *T*, one obtains

$$
\delta = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -a\_1 \Omega\_r \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & a\_2 \Omega\_r & 0 & 0 \end{bmatrix} \delta + \begin{bmatrix} 0 & 0 & 0 \\ b\_1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & b\_2 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & b\_3 \end{bmatrix} \tau + \begin{bmatrix} 0 \\ I\_1 \theta \dot{\phi} \\ 0 \\ 0 \\ I\_2 \psi \phi \\ 0 \\ I\_3 \dot{\psi} \theta \end{bmatrix}. \tag{10}
$$

Therefore, we can represent the rotational dynamics as a particular nonlinear system of the form

$$
\dot{\delta} = A\delta + B\tau + \varrho(\delta, t),
$$

where *A* and *B* are constant matrices, *B* is chosen such that (*<sup>A</sup>*, *B*) is controllable, and *ϕ*(*<sup>δ</sup>*, *t*) denotes the nonlinearities of the system. The parameters are given by:

$$I\_1 = \frac{I\_{zz} - I\_{yy}}{I\_{xx}}, \ I\_2 = \frac{I\_{xx} - I\_{zz}}{I\_{yy}}, \ I\_3 = \frac{I\_{yy} - I\_{xx}}{I\_{zz}}, \ a\_1 = \frac{I\_r}{I\_{yy}}, \ a\_2 = \frac{I\_r}{I\_{zz}}, \ b\_1 = \frac{1}{I\_{xx}}, \ b\_2 = \frac{d}{I\_{yy}}, \ b\_3 = \frac{d}{I\_{zz}}.$$

Now, we consider that the following assumptions and lemmas are accomplished.

**Assumption 1.** *The function ϕ*(*<sup>δ</sup>*, *t*) *is Lipschitz w.r.t. δ, with a Lipschitz constant γ, if:*

$$\|\|\boldsymbol{\varrho}(\delta,t) - \boldsymbol{\varrho}(\bar{\delta},t)\|\| \le \gamma \|\|\delta - \bar{\delta}\|. \tag{11}$$

**Lemma 1.** *Let ϕ* : [*a*¯, ¯ *b*] ×D→ R*n be continuous for some domain* D ⊂ R*n. Suppose that ∂ϕ∂δ*(*δ*, *t*) *exists and is continuous on* [*a*¯, ¯ *b*] × D*. If, for a convex subset W* ⊂ D*, there is a constant γ* - 0 *such that*

$$\left\| \frac{\partial \boldsymbol{\varrho}}{\partial \boldsymbol{\delta}}(t, \boldsymbol{\delta}) \right\| \leqslant \gamma$$

*on* [*a*¯, ¯ *b*] × *W, then*

*ϕ*(*<sup>δ</sup>*, *t*) − *<sup>ϕ</sup>*(˜*δ*, *t*) ≤ *γ<sup>δ</sup>* − ˜*δ*

*for all t* ∈ [*a*¯, ¯ *b*], *δ* ∈ *W, and* ˜*δ* ∈ *W.*

**Lemma 2.** *If ϕ*(*<sup>δ</sup>*, *t*) *and ∂ϕ∂δ*(*δ*, *t*) *are continuous on* [*a*¯, ¯*b*] × R*n, then ϕ is globally Lipschitz in δ on* [*a*¯, ¯ *b*] × R*n if and only if ∂ϕ∂δ* (*δ*, *t*) *is uniformly bounded on* [*a*¯, ¯*b*] × R*n.*

**Lemma 3.** *If ϕ*(*<sup>δ</sup>*, *t*) *and ∂ϕ∂δ* (*δ*, *t*) *are continuous on* [*a*¯, ¯*b*] × D*, for some domain* D ⊂ R*n, then ϕ is locally Lipschitz in δ on* [*a*¯, ¯ *b*] × D [25]*.*

**Assumption 2.** *For a and b* ∈ R*n and ε* > 0 *we have*

$$
\varepsilon \mathcal{L} a^T b \le \varepsilon^{-1} a^T a + \varepsilon b^T b. \tag{12}
$$

#### *3.1. Classical LMIs for the Quadrotor's Orientation*

We propose a controller *u* = −*Kδ* based on an LMI. Firstly, one introduces a classical LMI approach, which states the following lemma.

**Lemma 4.** *Consider system (10) and Assumption 1. Assume also that the nonlinearity ϕ is locally Lipschitz (Lemma 3) with Lipschitz constant γ (Lemma 1). Then, there exist a constant ϑ* > 0*, matrices X* = *X<sup>T</sup>* > 0*, and W with appropriate dimensions such that the following LMI is satisfied:*

$$
\begin{bmatrix}
A X + X A^T - B \mathcal{W}^T - \mathcal{W} B^T + 2(\theta - \gamma) X & \gamma I\_n + X \\
\gamma I\_n + X & -I\_n
\end{bmatrix} < 0,\tag{13}
$$

.

*with the feedback gain K* = *WTX*−1*. The system (10) is exponentially stable when the control input is u* = <sup>−</sup>*Kδ.*

We now need to obtain the value of the Lipschitz constant *γ* of the nonlinearities of the rotation system function *ϕ*(*<sup>δ</sup>*, *t*). From the rotational dynamics model (10), the function *ϕ*(*<sup>δ</sup>*, *t*) is given by

$$
\begin{bmatrix}
\dot{\varphi}(\delta,t)
\end{bmatrix} = \begin{bmatrix}
0 & I\_1\dot{\theta}\dot{\phi} & 0 & I\_2\dot{\psi}\dot{\phi} & 0 & I\_3\dot{\psi}\dot{\theta}
\end{bmatrix}^T.
$$

We are interested in calculating a Lipschitz constant for *ϕ*(*<sup>δ</sup>*, *t*) over the convex set

$$\mathcal{W} = \left\{ (\psi, \dot{\psi}, \theta, \dot{\theta}, \dot{\phi}, \dot{\phi}) \mid (|\psi| \le l\_{\theta'} |\dot{\psi}| \le l\_{\dot{\psi}'} |\theta| \le l\_{\theta'} |\dot{\theta}| \le l\_{\theta'} |\phi| \le l\_{\theta'} |\dot{\phi}| \le l\_{\dot{\phi}} \right\}.$$

Applying Lemma 1, one gets:

$$\begin{split} \sup\_{\delta \in \mathcal{W}} \left| \frac{\partial \varphi(\delta, t)}{\partial \delta} \right| &= \sup\_{\mathbf{x} \in \mathcal{W}} \left| \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & I\_1 \dot{\Phi} & 0 & I\_1 \dot{\theta} \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & I\_2 \dot{\Psi} & 0 & 0 & 0 & I\_2 \dot{\Phi} \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & I\_3 \dot{\theta} & 0 & I\_3 \dot{\Psi} & 0 & 0 \end{bmatrix} \right| \leqslant \gamma \\ &= \max \left\{ \left| \begin{aligned} |I\_1 \dot{\Phi}| + |I\_1 \dot{\theta}|, |I\_2 \dot{\Psi}| + |I\_2 \dot{\Phi}|, |I\_3 \dot{\theta}| + |I\_3 \dot{\Psi}| \right) \right\} \\ &= I\_2 \left\{ |\dot{\psi}| + |\dot{\phi}| \right\} = I\_3 \left\{ |\dot{\theta}| + |\dot{\psi}| \right\} \\ &\leqslant \quad I\_2 \left( l\_\psi + l\_\phi \right) = I\_3 \left( l\_\theta + l\_\phi \right) . \end{split} \right\}. \end{split}$$

Then,

$$\gamma = I\_2 \left( l\_{\psi} + l\_{\phi} \right) = I\_3 \left( l\_{\phi} + l\_{\psi} \right) \dots$$

The velocities are bounded as *lψ*˙ = *l* ˙*θ* = *lφ*˙ = 13*π*9 rad/s because the rotors are driven by DC permanent magne<sup>t</sup> motors, which support a maximum voltage of 9 V. This implies that the torque input vector *τ* is bounded and that the rotation speed capability of the motors has a maximum. When 9 V is applied to the motor the angular speed reaches *lψ*˙ = *l* ˙*θ* = *lφ*˙ = 13*π*9 rad/s. Considering these bounded values and the values of the Table 1, we can compute the Lipschitz constant as *γ* = 26*π*9 . Then, *ϕ* is locally Lipschitz.

**Table 1.** Parameters.


With *γ* = 26*π*9 , we try to solve the LMI (13) using the LMI Toolbox-R in MATLAB-R software. However, the LMI is infeasible. It may happen that the classical LMI controller can deal with the problem only for smaller values of the Lipschitz constant. Thus, we propose an LMI following some ideas from [26]. The LMI presented in [26] is conceived for observer design. The new LMI design techniques are significantly less conservative than the classical LMI design technique.

#### *3.2. Rotational Subsystem Control*

Based on the linear-state-feedback approach, the control law for the attitude dynamics is *u*(*t*) = −*K* ˜ *<sup>δ</sup>*(*t*). Therefore, the orientation dynamics can be constructed as follows:

˙

$$\begin{split} \tilde{\delta} &= A\tilde{\delta}(t) + Bu(t) + q(\tilde{\delta}(t), t) \\ &= A\tilde{\delta}(t) - BK\tilde{\delta}(t) + q(\tilde{\delta}(t), t). \end{split} \tag{14}$$

The goal is to find a suitable *B* and *K*, such that *δ* → ˜ *δ*. Thus, the attitude error dynamics *e*(*t*) is represented as:

$$\begin{aligned} \dot{\varepsilon}(t) &= \delta - \dot{\delta} \\ &= (A - BK)\varepsilon(t) + q(\delta(t)) - q(\delta(t)). \end{aligned} \tag{15}$$

The following assumptions are needed for the derivation of the control law.

**Assumption 3.** *There exists a matrix G with appropriate dimensions such that:*

$$\|\|\boldsymbol{\varrho}(\delta, t) - \boldsymbol{\varrho}(\bar{\delta}, t)\|\| \le \|\|\boldsymbol{G}(\delta - \bar{\delta})\|\|. \tag{16}$$

*This matrix G is a sparsely populated matrix. G*(*δ* − ˜*δ*) *can be much smaller than the constant γ*(*<sup>δ</sup>* − ˜*δ*) *used earlier in Equation (11) for the same nonlinear function.*

Let us now consider a larger Lipschitz constant of the nonlinear system. We can achieve a state feedback controller that is able to bring the state of the nonlinear system *δ*(*t*) to the desired state ˜*<sup>δ</sup>*(*t*). This controller is given in the following statement:

**Theorem 1.** *For attitude error dynamics (15), assume that Assumption 1, Lemmas 1 and 3 are satisfied and there exist a constant ε* > 0*, matrices X* = *X<sup>T</sup>* > 0*, and W with suitable dimensions, such that the following LMI holds:*

$$
\begin{bmatrix}
A\ X + XA^T - B\mathcal{W}^T - \mathcal{W}B^T + 2\theta X + XG^T + GX & \varepsilon I\_{\mathbb{H}} + XG^T \\
\varepsilon I\_{\mathbb{H}} + GX & -\varepsilon I\_{\mathbb{H}}
\end{bmatrix} < 0,\tag{17}
$$

*where In denotes an identity matrix with appropriate dimensions, ϑ* > 0 *is a constant, and W* = (*KTX*)−1*, the matrix K is a suitable feedback gain. Then, system (15) is exponentially stable, implying that the systems (10) and (14) are exponentially stable, then δ*(*t*) → ˜*<sup>δ</sup>*(*t*)*.*

**Proof.** Define the Lyapunov function *V* = *<sup>e</sup>TPe*. From the trajectory error (15), one gets:

$$\begin{split} \dot{V} &= \boldsymbol{\varepsilon}^{T} P \dot{\boldsymbol{\varepsilon}} + \dot{\boldsymbol{\varepsilon}}^{T} P \boldsymbol{\varepsilon} \\ &= \boldsymbol{\varepsilon}^{T} \left( P (A - BK) + (A - BK)^{T} P \right) \boldsymbol{\varepsilon} + 2 \boldsymbol{\varepsilon}^{T} P \left( \boldsymbol{\varrho} (\delta(t)) - \boldsymbol{\varrho} (\delta(t)) \right) . \end{split} \tag{18}$$

From Assumption 3, one gets:

$$\begin{split} \| 2\boldsymbol{\varepsilon}^{T}(\boldsymbol{t})P(\boldsymbol{\varrho}(\boldsymbol{u}(t),\boldsymbol{\delta}(\boldsymbol{t})) - \boldsymbol{\varrho}(\boldsymbol{u}(t),\boldsymbol{\delta}(\boldsymbol{t}))) \| &\leq 2\| \boldsymbol{\varepsilon}^{T}(\boldsymbol{t})P\| \| \| \boldsymbol{\varrho}(\boldsymbol{\delta}(\boldsymbol{t})) - \boldsymbol{\varrho}(\boldsymbol{\delta}(\boldsymbol{t})) \| \\ &\leq 2\| \boldsymbol{\varepsilon}^{T}(\boldsymbol{t})P\| \| \| \boldsymbol{\mathcal{G}} \boldsymbol{\varepsilon} \|. \end{split} \tag{19}$$

According to Assumption 2, *a* = *Pe*(*t*) and *b* = *Ge*, one can rewrite (19) as follows:

$$\left\| \left| 2\mathbf{e}^{T}(t)P(\boldsymbol{\varphi}(\delta(t)) - \boldsymbol{\varphi}(\tilde{\delta}(t))) \right| \right\| \leq \boldsymbol{\varepsilon}\_{1}^{-1} \mathbf{e}^{T} P \mathbf{P} \mathbf{e} + \boldsymbol{\varepsilon}\_{1} \mathbf{e}^{T} \mathbf{G}^{T} \mathbf{G} \mathbf{e}.\tag{20}$$

Now, replacing (20) into (18), one gets:

$$\dot{V} \le e^T(t) \left( P(A - BK) + (A - BK)^T P + \mathfrak{e}\_1^{-1} P P + \mathfrak{e}\_1 G^T G \right) e(t). \tag{21}$$

If *V* ˙ ≤ −2*ϑe<sup>T</sup>*(*t*)*Pe*(*t*) < 0, where *ϑ* > 0, one can rewrite (21) as:

$$\dot{V} \le e^T(t) \left( P(A - BK) + (A - BK)^T P + \varepsilon\_1^{-1} P P + \varepsilon\_1 G^T G + 2\theta P \right) e(t). \tag{22}$$

Indeed, the attitude error dynamics (15) is exponentially stable, and hence the two coupled systems (10) and (14) are exponentially stable. Using the Schur complement, Equation (22) can be easily represented in an LMI as:

$$
\begin{bmatrix}
P(A - BK) + (A - BK)^T P + 2\theta P + G^T P + PG & P + \varepsilon\_1 G^T \\
P + \varepsilon\_1 G^T & -\varepsilon\_1 I\_{\text{fl}}
\end{bmatrix} < 0. \tag{23}
$$

Multiplying the above inequality by *P*−<sup>1</sup> 0 0 <sup>−</sup><sup>1</sup> 1 *In* from the left-hand and right-hand sides, respectively, and letting *X* = *<sup>P</sup>*−1, *ε* = <sup>−</sup><sup>1</sup> 1 , and *W* = (*KP*−<sup>1</sup>)*<sup>T</sup>*, then the above inequality is further transformed into the following LMI:

$$
\begin{bmatrix}
A \ X + XA^T - B\mathcal{W}^T - \mathcal{W}B^T + 2\theta X + XG^T + GX & \varepsilon I\_{\mathbb{H}} + XG^T \\
\varepsilon I\_{\mathbb{H}} + GX & -\varepsilon I\_{\mathbb{H}}
\end{bmatrix} < 0. \tag{24}
$$

If suitable *X* > 0 matrix and *W* are selected such that the LMI (17) is satisfied, then the attitude error dynamics (15) with the feedback gain *K* = *WTX*−<sup>1</sup> is exponentially stable, implying that the coupled systems (10) and (14) are exponentially synchronized.

Now, we apply the main result of this paper (Theorem 1) to system (10). Then, we compute a controller for the rotational dynamics of the vehicle with guaranteed stability. Using Theorem 1, the LMI is then solved to obtain the control gain matrix *K* with *ε* = 1, *ϑ* = 25, and *γ* = 26*π*9 . Therefore, one easily obtains *K* from (24) by using the MATLAB-R LMI Toolbox-R :

$$K = \begin{bmatrix} 3.281 \times 10^8 & 2.235 \times 10^6 & 0 & 0 & 0 & 0\\ 0 & 0 & 38.75 & 1.16 & -0.0114 & -2.27 \times 10^{-4} \\ 0 & 0 & 0.0114 & 2.27 \times 10^{-4} & 38.75 & 1.16 \end{bmatrix} . \tag{25}$$

The results of the state feedback controller with the gain matrix (25) are shown in Figures 2 and 3. Figure 2 displays the quadrotor attitude (roll, pitch, and yaw) with *δ*0 = 0◦ 0◦ 1◦ 0◦ 2◦ 0◦ . Note that the stabilization time is about 0.2 s, thus the state feedback controller with the gain matrix *K* calculated by LMI provides good transient performance. Figure 3 shows the input torques. Here we can observe that they are smooth.

**Figure 2.** Quadrotor attitude.

**Figure 3.** Control inputs.

#### **4. Energy-Based Control**

In this section, we propose an energy-based strategy to control the quadrotor's translation movements and to attenuate the payload oscillation.

#### *4.1. Planar Case*

In order to show the energy-based control law design in a simple manner, let us consider the system as a planar system operating in the *XZ* axis, as illustrated in Figure 4. The controller synthesis for the three-dimensional case is presented in Section 4.2.

In this case, *y* = *ψ* = *φ* = *<sup>α</sup>y* = 0, then *q*˜ = *x z α T*, *u*˜ = *Fx Fz T* = *f* sin *θ f* cos *θ T* and the matrices of the model (5) for the translational dynamics are defined as

$$
\tilde{\mathcal{M}}(\vec{q}) = \begin{bmatrix} M+m & 0 & mlc\_{\mathfrak{a}} \\ 0 & M+m & mls\_{\mathfrak{a}} \\ mlc\_{\mathfrak{a}} & mls\_{\mathfrak{a}} & ml^2 \end{bmatrix}, \quad \tilde{\mathcal{C}}(\vec{q}, \dot{\vec{q}}) = \begin{bmatrix} 0 & 0 & -mls\_{\mathfrak{a}}\dot{\mathfrak{a}} \\ 0 & 0 & mlc\_{\mathfrak{a}}\dot{\mathfrak{a}} \\ 0 & 0 & 0 \end{bmatrix} \tag{26}
$$

$$\mathcal{G}(\vec{q}) = \begin{bmatrix} 0 \\ (M+m)\mathcal{g} \\ m\mathcal{l}\mathcal{g}s\_n \end{bmatrix}, \quad \mathcal{B}\vec{u} = \begin{bmatrix} F\_x \\ F\_z \\ 0 \end{bmatrix}. \tag{27}$$

**Figure 4.** Two-dimensional quadrotor with a cable-suspended payload.

The total energy of the translational dynamics can be described by

$$\begin{aligned} H(\vec{q}, \dot{\vec{q}}) &= \frac{1}{2} \dot{\vec{q}}^T \mathcal{M} \dot{\vec{q}} + V(\vec{q}) \\ &= \frac{1}{2} \dot{\vec{q}}^T \mathcal{M} \dot{\vec{q}} + (M + m) gz - mgl \cos(a). \end{aligned} \tag{28}$$

Differentiating (28) along the trajectories of the system, we obtain

$$
\hat{H}(\vec{q}, \dot{\vec{q}}) = \dot{\vec{q}}^T \mathcal{N}(\vec{q}) \ddot{\vec{q}} + \frac{1}{2} \dot{\vec{q}}^T \dot{\mathcal{N}}(\vec{q}) \dot{\vec{q}} + \dot{\vec{q}}^T \mathcal{G}(\vec{q}).
$$

Substituting (5) into the above yields

$$\dot{H}(\vec{q}, \dot{\vec{q}}) = \dot{\vec{q}}^T \left( \vec{n} - \vec{\mathcal{C}}(\vec{q}, \dot{\vec{q}}) \dot{\vec{q}} - \vec{\mathcal{G}}(\vec{q}) \right) + \frac{1}{2} \dot{\vec{q}}^T \dot{\mathcal{N}}(\vec{q}) \dot{\vec{q}} + \dot{\vec{q}}^T \vec{\mathcal{G}}(\vec{q}) \dots$$

Taking into account that the skew-symmetric relationship ˙*q*˜*<sup>T</sup>* 12 ˙ M˜ (*q*˜) ˙*q*˜ − *<sup>C</sup>*˜(*q*˜, ˙*q*˜)!˙*q*˜ = 0 is satisfied, we obtain

$$\begin{aligned} \dot{H}(\bar{q}, \dot{\bar{q}}) &= \quad \dot{\bar{q}}^T \bar{u} \\ &= \quad \dot{\mathfrak{x}} F\_{\mathfrak{x}} + \dot{\mathfrak{z}} F\_{\mathfrak{z}}. \end{aligned}$$

Considering *x*¯ = *x* − *x*˜ and *z*¯ = *z* − *z*˜, the total energy in terms of the error is given by

$$
\hat{H} = \dot{\mathfrak{X}} F\_{\mathfrak{x}} + \dot{\mathfrak{z}} F\_{\mathfrak{z}}.\tag{29}
$$

We propose the following Lyapunov candidate function:

$$E = \frac{1}{2}\ddot{H}^2 + \frac{k\_{vx}}{2}\hat{\mathbf{x}}^2 + \frac{k\_{vx}}{2}\hat{\mathbf{z}}^2 + \frac{k\_{px}}{2}\hat{\mathbf{x}}^2 + \frac{k\_{pz}}{2}\hat{\mathbf{z}}^2,\tag{30}$$

where *kpx*, *kpz* are proportional constant gains and the *kvx*, *kvz* constants inject damping into the system. Differentiating (30) with respect to time, we have

$$\begin{split} \dot{E} &= \quad \dot{H}\dot{H} + k\_{\upsilon x}\hat{x}\tilde{x} + k\_{\upsilon z}\hat{z}\tilde{z} + k\_{px}\bar{x}\hat{x} + k\_{pz}\bar{z}\hat{z} \\ &= \quad \pounds \left( \bar{H}F\_{\text{x}} + k\_{\upsilon x}\mathbb{1} + k\_{px}\mathfrak{x} \right) + \hat{z} \left( \bar{H}F\_{z} + k\_{\upsilon z}\sharp + k\_{\upsilon z}\sharp \right). \end{split} \tag{31}$$

We can obtain *x*¨¯ and *z*¨¯ from (5), (26), and (27):

$$\begin{array}{rcl} \tilde{x} &=& \frac{M\left(lm\sin\pi\hbar\dot{\alpha}^2 + F\_x\right) + F\_x m\cos^2\pi + F\_z m\cos\pi\sin\pi}{M^2 + Mm},\\ \tilde{z} &=& -g + \frac{F\_z m + M\left(F\_z - lm\hbar^2\cos\delta\right) - F\_z m\cos^2\delta + F\_x m\cos\delta\sin\delta}{M^2 + Mm}.\end{array}$$

.

Introducing the above into (31), we ge<sup>t</sup>

$$\begin{split} E^{\ddagger} &= -\hbar \left[ H F\_{\ddagger} + \frac{k\_{\square \Omega} M \left( \left| m \sin \hbar \hat{\mathbf{n}} \right|^{2} + F\_{\ddagger} \right) + k\_{\square \Omega} m F\_{\ddagger} \cos^{2} \tilde{\mathbf{n}} + k\_{\square \Omega} m F\_{\ddagger} \cos \hbar \sin \tilde{\mathbf{n}}}{M(M+m)} + k\_{\square \Omega} \tilde{\mathbf{g}} \right] + \\ &\quad \hbar \left[ H F\_{\ddagger} + \frac{k\_{\square \Omega} M \left( -\left| m \cos \hbar \hat{\mathbf{n}} \right|^{2} + F\_{\ddagger} \right) + k\_{\square \Omega} m F\_{\ddagger} - k\_{\square \Omega} m F\_{\ddagger} \cos^{2} \tilde{\mathbf{n}} + k\_{\square \Omega} m F\_{\ddagger} \cos \hbar \sin \tilde{\mathbf{n}}}{M(M+m)} + k\_{\square \Omega} \tilde{\mathbf{z}} - k\_{\square \Omega} \tilde{\mathbf{g}} \right]. \end{split}$$

We propose a control law such that

$$
\begin{bmatrix}
E + \frac{k\_{\text{rx}}}{M+m} + \frac{k\_{\text{rx}}m\cos^2 a}{M(M+m)} & \frac{k\_{\text{rx}}m\cos a \sin a}{M(M+m)} \\
\frac{k\_{\text{rx}}m\cos a \sin a}{M(M+m)} & E + \frac{k\_{\text{rx}}}{M} - \frac{k\_{\text{rx}}m\cos^2 a}{M(M+m)}
\end{bmatrix}
\begin{bmatrix}
F\_x \\
F\_z
\end{bmatrix} = 
\begin{bmatrix}
\frac{k\_{\text{rx}}lm\cos a\dot{\mathbf{z}}^2}{M+m} - k\_{\text{rx}}z + k\_{\text{rx}}g - k\_{\text{rx}}\dot{z}
\end{bmatrix}.
\tag{32}
$$

The matrix that multiplies the vector *Fx Fz T* is a nonsingular matrix. Which leads to

$$
\dot{E} = -k\_{ix}\dot{x}^2 - k\_{iz}\dot{z}^2.
$$

From (32) we can obtain

$$F\_x = \frac{-k\_{ix}\dot{x} - k\_{xx}\dot{x} - \frac{k\_{xx}l\cos\Delta t}{M+m} - \frac{k\_{ix}m\cos\Delta t\sin\left(\frac{k\_{xx}l\cos\Delta t\dot{h}}{M+m} - k\_{ix}z + k\_{ix}g - k\_{ix}\frac{\dot{x}}{\dot{x}}\right)}{M+\frac{k\_{xx}}{M+m} + \frac{k\_{xx}m\cos^2\Delta t}{M(M+m)} - \frac{k\_{ix}k\_{ix}m^2\cos^2\Delta t\sin^2\Delta t}{M(M+m)\left(M(M+m)\dot{H} + k\_{ix}(x+M) - k\_{ix}m\cos^2\Delta t\right)}},$$

$$F\_z = \frac{-k\_{ix}\dot{z} - k\_{xz}\ddot{z} + \frac{k\_{ux}l\cos\Delta t\dot{h}}{M+m} + k\_{vz}g + \frac{k\_{vx}m\cos\Delta t\sin\left(\frac{k\_{vx}l\sin\Delta t\dot{h}}{M+m} + k\_{vx}x + k\_{ix}x\right)}{M+\frac{k\_{ux}}{M} - \frac{k\_{ux}l\cos\Delta^2t\dot{h}}{M(M+m)} - \frac{k\_{ux}k\_{vx}m^2\cos^2\Delta t\sin^2\Delta t}{M(M+m)(M+m)\left(M(M+m)\dot{H} + k\_{vx}M + k\_{vx}m\cos^2\Delta t\right)}}.$$

#### *4.2. Three-Dimensional Case*

The mathematical model for the three-dimensional case is presented in (5)–(9). The total energy of the system is

$$\begin{aligned} H(\vec{q}, \dot{\vec{q}}) &= \begin{aligned} \frac{1}{2} \dot{\vec{q}}^T \tilde{\mathcal{M}} \dot{\vec{q}} + V(\vec{q}) \\ &= \frac{1}{2} \dot{\vec{q}}^T \tilde{\mathcal{M}} \dot{\vec{q}} + (\mathcal{M} + m) \text{g}z - m \text{gl} \cos(a\_x) \cos(a\_y) . \end{aligned} \tag{33}$$

In a similar way to the planar case, differentiating (33) along the trajectories of the system, we obtain

$$\begin{aligned} \dot{H}(\vec{q}, \dot{\vec{q}}) &= \dot{\vec{q}}^T \vec{u} \\ &= \dot{\vec{x}} F\_x + \dot{y} F\_y + \dot{z} F\_z. \end{aligned} \tag{34}$$

The total energy in terms of the error can be rewritten as

$$
\dot{H} = \dot{\mathfrak{X}}F\_{\mathfrak{x}} + \dot{\mathfrak{Y}}F\_{\mathfrak{Y}} + \dot{\mathfrak{z}}F\_{\mathfrak{z}}.\tag{35}
$$

Consider the following Lyapunov candidate function:

$$E(\bar{\eta}, \hat{\eta}) = \frac{1}{2}\bar{H}^2 + \frac{k\_{vx}}{2}\hat{\mathbf{x}}^2 + \frac{k\_{vy}}{2}\hat{y}^2 + \frac{k\_{vz}}{2}\hat{z}^2 + \frac{k\_{px}}{2}\bar{\mathbf{x}}^2 + \frac{k\_{py}}{2}\bar{y}^2 + \frac{k\_{pz}}{2}\bar{z}^2. \tag{36}$$

Differentiating (36) along the trajectories of the system, it follows that

˙

$$\begin{split} \dot{E}\_{\perp} &= \; \dot{H}\dot{H} + k\_{\text{vx}}\dot{x}\ddot{x} + k\_{\text{vy}}\dot{y}\ddot{y} + k\_{\text{vx}}\dot{z}\dddot{x} + k\_{\text{px}}x\dddot{x} + k\_{\text{py}}y\dot{y} + k\_{\text{pz}}z\dot{z} \\ &= \; \dot{\mathfrak{x}}\left(\ddot{H}F\_{\text{x}} + k\_{\text{vx}}\ddot{x} + k\_{\text{px}}\ddot{x}\right) + \dot{\mathfrak{y}}\left(\ddot{H}F\_{\text{y}} + k\_{\text{vy}}\ddot{y} + k\_{\text{py}}\ddot{y}\right) + \dot{\mathfrak{z}}\left(\ddot{H}F\_{\text{z}} + k\_{\text{vz}}\ddot{z} + k\_{\text{pz}}\ddot{z}\right). \end{split} \tag{37}$$

From (5)–(9) we can obtain *x*¨¯, *y*¨¯ and *z*¨¯. These expressions are defined as:

*x* ¨ ¯ = *Fx* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *y c*4*α*¯ *x* + *c*2*α*¯ *x c*2*α*¯ *y*! + *Mc*4*α*¯ *y s*2*α*¯ *x*! − *Fymc*<sup>2</sup>*α*¯ *x c*2*α*¯ *y sα*¯ *x sα*¯ *y* + *Fzmsα*¯ *x cα*¯ *x c*3*α*¯ *y* + *<sup>λ</sup>*(*μ*, *μ*˙) (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! , *y* ¨ ¯ = *Fy* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *x c*4*α*¯ *y* + *c*2*α*¯ *x c*2*α*¯ *y*! + *Mc*4*α*¯ *x s*2*α*¯ *y*! − *Fxmc*<sup>2</sup>*α*¯ *x c*2*α*¯ *y sα*¯ *x sα*¯ *y* + *Fzmsα*¯ *y cα*¯ *y c*3*α*¯ *x* + *<sup>κ</sup>*(*μ*, *μ*˙) (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! , *z* ¨ ¯ = *Fz* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *y c*4*α*¯ *x* − *s*2*α*¯ *x c*4*α*¯ *y*! − *Mc*2*α*¯ *x c*2*α*¯ *y*! − *Fxmc*<sup>3</sup>*α*¯ *y cα*¯ *x sα*¯ *x* − *Fymsα*¯ *y cα*¯ *y c*3*α*¯ *x* + *M*<sup>2</sup> + *Mm go*(*μ*) + *<sup>χ</sup>*(*μ*, *μ*˙) (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! ,

where

*<sup>λ</sup>*(*α*¯ *x*, *α*¯ *y*) = *Mlm <sup>s</sup>α*¯ *x c*2*α*¯ *x c*4*α*¯ *y α*¯˙ *x*2 + *α*¯˙ *y*2! + *c*2*α*¯ *x c*2*α*¯ *y sα*¯ *x s*2*α*¯ *yα*¯˙ *y*2 + 2*cα*¯ *x c*3*α*¯ *y s*2*α*¯ *x sα*¯ *yα*¯˙ *xα*¯˙ *y* + *c*4*α*¯ *y s*3*α*¯ *xα*¯˙ *x*2! , *κ*(*μ*¯, *μ*¯˙) = *Mlm <sup>s</sup>α*¯ *y c*2*α*¯ *y c*4*α*¯ *x α*¯˙ *x*2 + *α*¯˙ *y*2! + *c*2*α*¯ *x c*2*α*¯ *y sα*¯ *y s*2*α*¯ *xα*¯˙ *x*2 + 2*cα*¯ *y c*3*α*¯ *x s*2*α*¯ *y sα*¯ *xα*¯˙ *xα*¯˙ *y* + *c*4*α*¯ *x s*3*α*¯ *yα*¯˙ *y*2! , *χ*(*μ*¯, *μ*¯˙) = *Mlm <sup>c</sup>*3*α*¯ *y c*3*α*¯ *x α*¯˙ *x*2 + *α*¯˙ *y*2! + *cα*¯ *x c*3*α*¯ *y s*2*α*¯ *xα*¯˙ *x*2 − 2*sα*¯ *y c*2*α*¯ *x c*2*α*¯ *y sα*¯ *xα*¯˙ *xα*¯˙ *y* − *c*3*α*¯ *x s*2*α*¯ *y cα*¯ *yα*¯˙ *y*2! , *o*(*μ*¯) <sup>=</sup> *<sup>c</sup>*4*α*¯ *y s*2*α*¯ *x* + *c*2*α*¯ *x c*2*α*¯ *y* − *s*2*α*¯ *y c*4*α*¯ *x*! .

Substituting *x*¨¯, *y*¨¯, and *z*¨¯ into (37) yields

*E* ˙ =*x* ¯ ˙ ⎡⎣*HF*¯ *x*+*kvx Fx* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *y c*4*α*¯ *x* + *c*2*α*¯ *x c*2*α*¯ *y*! *Mc*4*α*¯ *y s*2*α*¯ *x*! − *Fymc*<sup>2</sup>*α*¯ *x c*2*α*¯ *y sα*¯ *x sα*¯ *y* + *Fzmsα*¯ *x cα*¯ *x c*3*α*¯ *y* + *<sup>λ</sup>*(*μ*¯, *μ*¯˙) (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! +*kpx x*¯⎤⎦+ *y* ¯ ˙ ⎡⎣*HF*¯ *y*<sup>+</sup>*kvy Fy* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *x c*4*α*¯ *y* + *c*2*α*¯ *x c*2*α*¯ *y*! + *Mc*4*α*¯ *x s*2*α*¯ *y*! − *Fxmc*<sup>2</sup>*α*¯ *x c*2*α*¯ *y sα*¯ *x sα*¯ *y* + *Fzmsα*¯ *y cα*¯ *y c*3*α*¯ *x* + *κ*(*μ*¯, *μ*¯˙) (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! +*kpyy*¯⎤⎦+ *z* ¯ ˙ ⎡⎣*HF*¯ *z* + *kvz Fz* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *y c*4*α*¯ *x* − *s*2*α*¯ *x c*4*α*¯ *y*! − *Mc*2*α*¯ *x c*2*α*¯ *y*! − *Fxmc*<sup>3</sup>*α*¯ *y cα*¯ *x sα*¯ *x* − *Fymsα*¯ *y cα*¯ *y c*3*α*¯ *x* (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! +*kvz M*<sup>2</sup> + *Mm go*(*μ*¯) + *χ*(*μ*¯, *μ*¯˙) (*M*<sup>2</sup> + *Mm*) *<sup>s</sup>*4*α*¯ *x s*2*α*¯ *y* + *s*2*α*¯ *x s*4*α*¯ *y* − <sup>3</sup>*s*2*α*¯ *x s*2*α*¯ *y* + 1! + *kpzz*¯⎤⎦

.

We propose a control law such that

$$\dot{E} = -k\_{ix}\dot{x}^2 - k\_{iy}\dot{y}^2 - k\_{iz}\dot{z}^2 \dots$$

Finally, we solve the following system of equations to obtain *Fx*, *Fy*, and *Fz*.

*Fx* ⎛⎝*H*¯ + *kvx* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *y c*4*α*¯ *x* + *c*2*α*¯ *x c*2*α*¯ *y*! + *Mc*4*α*¯ *y s*2*α*¯ *x*! *<sup>δ</sup>*(*μ*¯) ⎞⎠ − *Fy kvxmc*<sup>2</sup>*α*¯ *x c*2*α*¯ *y sα*¯ *x sα*¯ *y <sup>δ</sup>*(*μ*¯) + *Fz kvxmsα*¯ *x cα*¯ *x c*3*α*¯ *y <sup>δ</sup>*(*μ*¯) + *<sup>λ</sup>*(*μ*¯, *μ*¯˙) *<sup>δ</sup>*(*μ*¯) + *kpx x*¯ = −*kix x*¯˙, (38) *Fy* ⎛⎝*H*¯ + *kvy* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *x c*4*α*¯ *y* + *c*2*α*¯ *x c*2*α*¯ *y*! + *Mc*4*α*¯ *x s*2*α*¯ *y*! *<sup>δ</sup>*(*μ*¯) ⎞⎠ − *Fx kvymc*<sup>2</sup>*α*¯ *x c*2*α*¯ *y sα*¯ *x sα*¯ *y <sup>δ</sup>*(*μ*¯) + *Fz kvymsα*¯ *y cα*¯ *y c*3*α*¯ *x <sup>δ</sup>*(*μ*¯) + *κ*(*μ*¯, *μ*¯˙) *<sup>δ</sup>*(*μ*¯) + *kpyy*¯ = −*kiyy*¯˙, (39) *Fz* ⎛⎝*H*¯ + *kvz* (*M* + *<sup>m</sup>*) *<sup>s</sup>*2*α*¯ *y c*4*α*¯ *x* − *s*2*α*¯ *x c*4*α*¯ *y*! − *Mc*2*α*¯ *x c*2*α*¯ *y*! *<sup>δ</sup>*(*μ*¯) ⎞⎠ − *Fx kvzmc*<sup>3</sup>*α*¯ *y cα*¯ *x sα*¯ *x <sup>δ</sup>*(*μ*¯) − *Fy kvzmsα*¯ *y cα*¯ *y c*3*α*¯ *x <sup>δ</sup>*(*μ*¯) +*M*<sup>2</sup> + *Mm go*(*μ*¯) *<sup>δ</sup>*(*μ*¯) + *χ*(*μ*¯, *μ*¯˙) *<sup>δ</sup>*(*μ*¯) + *kpzz*¯ = −*kizz*¯˙, (40)

where

$$\mathcal{S}(\boldsymbol{\mu}) = \left(\boldsymbol{M}^2 + \boldsymbol{M}m\right) \left(\mathbf{s}\_{\bar{\alpha}\_x}^4 \mathbf{s}\_{\bar{\alpha}\_y}^2 + \mathbf{s}\_{\bar{\alpha}\_x}^2 \mathbf{s}\_{\bar{\alpha}\_y}^4 - \mathbf{3} \mathbf{s}\_{\bar{\alpha}\_x}^2 \mathbf{s}\_{\bar{\alpha}\_y}^2 + 1\right).$$

#### **5. Numerical Simulations and Results**

In order to check the performance of the designed control scheme, some simulations were carried out. The objective was to move the vehicle transporting a payload to the desired position of a square of 1 m length at 1 m height. The desired trajectory is then defined by,

$$\xi\_d = \begin{cases} \begin{bmatrix} 0 & 0 & 1 \end{bmatrix}^T & \text{,} & t & < 2 \\ \begin{bmatrix} 0 & 1 & 1 \end{bmatrix}^T & \text{,} & 2 << & t & < 14 \\ \begin{bmatrix} 1 & 1 & 1 \end{bmatrix}^T & \text{,} & 14 << & t & < 26 \\ \begin{bmatrix} 1 & 0 & 1 \end{bmatrix}^T & \text{,} & 26 << & t & < 38 \\ \begin{bmatrix} 0 & 0 & 1 \end{bmatrix}^T & \text{,} & 38 << & t & < = 50 \end{cases}$$

The values of the model parameters used for the simulation were the following: *M* = 0.4 kg, *m* = 0.03 kg, *l* = 0.35 m, *g* = 9.81 m/s2.

These parameters are close to real aerial platforms. The corresponding simulation results are presented in the following figures.

Figure 5 illustrates the *x*, *y*, and *z* positions of the vehicle during the validation. We can see that the position for each axis was stabilized according to the desired reference points. The position *z* was regulated in less than 2 s, the performance of the *x* and *y* positions dynamics were similar and were regulated in less than 5 s.

Figure 6 displays the payload swing angles *α* and *β*. It is clear that the proposed control law exhibited good performance, since the payload swing angles were regulated to 0◦ at around 5 s and the maximum overshoot was ±2.2◦.

The simulation results of the feedback controller based on LMI (inner-loop of the system) show the quadrotor's orientation dynamics in Figure 7. We can see that the attitude converged to the desired

points with a null steady-state error. The evolution of the control inputs *f* , *τψ*, *τθ*, and *τφ* are presented in Figure 8. Finally, a three-dimensional view of the path followed by the vehicle is depicted in Figure 9.

**Figure 5.** Quadrotor position.

**Figure 7.** Quadrotor attitude.

**Figure 9.** Three-dimensional trajectory.

One more numerical experiment was carried out. Figures 10–13 illustrate the tracking of an ascending circular trajectory to prove the efficiency of the proposed controller in a scenario that involves simultaneous variations of both *α* and *β* angles. These figures show that the proposed control strategy was capable of achieving accurate trajectory tracking since the aircraft converged to the reference trajectory while attenuating the swing angles of the payload.

In summary, these numerical experiments show that the proposed control scheme presented a satisfactory performance in position control and the attenuation of cable-suspended payload swing. It succeeded in transporting the payload to a desired position with attenuation of the oscillation angles. In contrast, the algorithm of [18] involves solving complicated partial differential equations for obtaining the control law, and they cannot be solved for the 3D case.

**Figure 10.** Quadrotor position.

**Figure 12.** Quadrotor attitude.

**Figure 13.** Three-dimensional trajectory.

## **6. Conclusions**

This work presents an energy-based control strategy and a nonlinear state feedback controller based on a linear matrix inequality to solve the problem of transporting a cable-suspended payload by an unmanned aerial vehicle. On one hand, a new methodology based on an LMI for stabilizing the orientation dynamics is proposed. The main contribution is that we can employ the proposed methodology for Lipschitz nonlinear systems with larger Lipschitz constants than other classical techniques based on LMIs. Moreover, the LMI-based controller results in a control algorithm with relative simplicity and guaranteed stability. In addition, the design of the LMI-based control takes into account physical limits of the system such as the maximum motor voltage or its rotation speed capability through the velocity and torque bounds which are used to calculate the Lipschitz constant, while with the energy-based control these limits are not part of the controller design. On the other hand, an energy-based control to stabilize the quadrotor's translational dynamics and to attenuate the cable-suspended payload swing is designed in this work. This strategy is based on an energy approach and the passivity properties of the translational dynamics. Passivity-based control is employed, as this part of the system is underactuated. The main contribution is that the computation of excessive and complex partial differential equations is not needed to obtain the control law. The results showed an excellent performance of the proposed control scheme. Thus, the new approach achieves precise payload positioning with rapid oscillation attenuation.

Future work will extend the energy-based control method in order to consider variations in the payload. In addition, the methodology will be extended so that the linear matrix inequality can be replaced by an algebraic Riccati equation.

**Author Contributions:** Conceptualization, M.-E.G.-S. and O.H.-G.; methodology, M.-E.G.-S., O.H.-G., R.L., C.-D.G.-B., G.V.-P., and F.-R.L.-E.; formal analysis, M.-E.G.-S. and O.H.-G.; writing—original draft preparation, M.-E.G.-S. and O.H.-G.; writing—review and editing, G.V.-P. and F.-R.L.-E.; supervision, R.L. and C.-D.G.-B.

**Funding:** This work has been partially supported by PRODEP Mexican Program gran<sup>t</sup> number 511-6/18-10333.

**Conflicts of Interest:** The authors declare no conflict of interest.
