A Predictive Control Strategy for Aerial Payload Transportation with an Unmanned Aerial Vehicle

Abstract

This paper presents the results of a model-based predictive control (MPC) design for a quadrotor aerial vehicle with a suspended load. Unlike previous works, the controller takes into account the hanging payload dynamics, the dynamics in three-dimensional space, and the vehicle rotation, achieving a good balance between fast stabilization times and small swing angles. The mathematical model is based on the Euler–Lagrange formulation and considers the dynamics of the vehicle, the cable, and the load. Then, the mathematical model is represented as an input-affine system to obtain the linear model for the control design. A constrained MPC strategy was designed and compared with an unconstrained MPC and an algorithm from the literature for the case of study. The constraints to be considered include the limits on the swing angles and the quadrotor position. The constrained control algorithm was constructed to stabilize the aerial vehicle. It aims to track a trajectory reference while attenuating the load swing, considering a maximum swing range of $\pm {10}^{\circ }$. Numerical simulations were carried out to validate the control strategy.

1. Introduction

In recent years, the development of aerial robotics has been rapid via both scientific and commercial research, since the use of these autonomous systems is now not restricted to the military field. The versatility of unmanned aerial vehicles (UAVs), together with recent technological advances, has allowed the development of new devices that have a higher degree of autonomy; hence, they considered for different types of applications, mainly in those where the user’s safety may be compromised or where the areas of interest are difficult to access. The transport of a load suspended by a cable is one of the most common activities with unmanned aerial vehicles. The modeling and design of control strategies for these systems is a topic that has attracted the attention of researchers in recent decades; this is because there is a variety of tasks where fast and efficient handling of a load is required, which leads to aggressive movements that generate oscillations in the load that when not controlled can cause accidents. The works focused on the use of UAVs for payload transport are relatively recent [1].

There are several devices reported in the literature that are used to attach a load to such a vehicle. For example, in [2,3] a clamp-type mechanical device with only one degree of freedom that keeps the load attached to the structure of the UAV was used. In [4] the use of a manipulator with a robotic arm was mentioned. In these configurations, the load is closer to the vehicle’s center of gravity, which causes increases in weight and rotational inertia that affect the maneuverability of the vehicle in flight. Another approach, which is widely used in transportation, consists of suspending the payload through a rope or cable connected directly to the vehicle’s structure, allowing fast and agile movements. This approach adds degrees of freedom to the system, and a swinging motion of the payload that acts as a pendulum that can transmit oscillations through the cable and negatively affect the operations and stability of the vehicle in flight [5].

In the literature, different control strategies have been proposed to address the general problem of transporting a payload suspended by a cable from a UAV. For example, in [6] the effect of dynamic load disturbances introduced by instantaneously increased load mass and how those affect UAVs under PID flight control was studied. Ref. [7] introduced the use of dynamic programming to generate an optimal trajectory and adaptive control for maneuvering without swinging effects. Ref. [8] presented a nonlinear dynamic model and two IDA-PBC control algorithms for a quadrotor with a cable-suspended payload that explicitly incorporate total energy-shaping; however, the vehicle travels only in the x–z plane and then in the y–z plane to achieve control in three-dimensional space. Additionally, ref. [9] proposed a mathematical model of the interconnected multi-body system using Kane’s equations and a tracking controller based on the backstepping technique. Ref. [10] considered a hierarchical scheme combining energy-based control (translational dynamics) and a nonlinear state feedback controller based on a linear matrix inequality (rotational dynamics) to solve the problem of transporting a cable-suspended payload by a UAV. Meanwhile, in [11] two novel nonlinear control algorithms for controlling a quadcopter’s position and attenuating the swing angle of the payload were designed. The controllers take advantage of the natural coupling existing between the horizontal quadcopter movement and the payload oscillation; however, the vehicle only travels in one plane at a time. Ref. [12] shows a nonlinear control approach with an elaborately constructed integral term for an aerial transportation systems. This method achieves satisfactory anti-swing and positioning performance and reduces steady errors during real flight. A model of a quadrotor with a suspended payload was derived in [13]. The hierarchical control structure of inner-loop attitude and outer-loop position using PID was adopted in the design of trajectory tracking of the quadrotor transporting a payload. Ref. [14] presents a quadrotor with the payload connected to the vehicle through a cable. The payload swing is rejected by using a nonlinear controller based on a system model derived using the Euler–Lagrange equations. However, the payload is considered a disturbance. Ref. [15] presents a model and a proportional-derivative control of a quadrotor with payload uncertainties; the payload is also considered a disturbance. A mathematical model of the payload and quadrotor was presented in [16] through the Euler–Lagrange formulation. In addition, a robust sliding mode control was designed for the transitional movements to cope with both disturbances such as wind, and payload swings. Ref. [17] proposes an adaptive control technique to control the velocity of a volumetric payload with unknown mass, which is transported by a quadrotor in the presence of rotor down-wash.

Optimal control techniques have also been applied. For example, ref. [18] developed an iterative linear quadratic regulator (iLQR) optimal strategy of a mathematical model for a quadrotor with a cable-suspended heavy rigid body. The algorithm considers an anti-swinging load in a transporting task; however, the stabilization times are very long compared to the results of works found in the literature and to those obtained in the present article. Reference [19] explores a trajectory planning method based on predictive control for a quadrotor with a suspended payload. The authors contemplated the anti-swing and obstacle avoidance in the flight path. However, this research did not consider the rotational dynamics of the vehicle and treated the quadrotor as a three-dimensional overhead crane. Moreover, ref. [20] presented a linear quadratic Gaussian control method for a quadrotor transporting an unknown suspended load: the load parameters were unknown, and the load state was not available. This paper also did not consider the rotational dynamics of the vehicle. Reference [21] presents a cascade control strategy with anti-disturbance functionality and a predictive optimal function to control the rotation of a quadrotor carrying unknown payloads. However, the translational control was not realized and the payload was considered a disturbance.

The primary mission of aerial transportation is to maintain the hanging load within the desired operating and safety limits. Therefore, the prime aim is to design control algorithms that ensure minimal swing angles without sacrificing response times. However, some works listed above presented swing angles between $\pm {12}^{\circ }$ and $\pm {35}^{\circ }$, or in order to achieve small angles, they employed very conservative control laws, resulting in long settling times—e.g., [18,22]. Thus, motivated by the previous works—where the authors dealt with the hanging payload transportation problem by (i) considering the payload as a disturbance; (ii) including the payload dynamics, but considering only a longitudinal translation of the UAV; (iii) considering three-dimensional dynamics but neglecting the rotation of the vehicle; or (iv) reducing the payload swing by sacrificing performance—in this paper, we offer an alternative approach to stabilize a quadrotor transporting a payload by using a predictive control strategy. Here, the quadrotor-payload system model includes the hanging payload dynamics, and the predictive controller considers both the system’s translation dynamics in three-dimensional space and the vehicle’s rotation. Constraints on the load swing angles are to keep the vehicle transporting a load within the desired operating and safety limits. Considering constraints on the quadrotor positioning in x and y allows one to avoid overshoot in the transient responses of this positioning and thus to evade further oscillations of the load angles. Predictive control was chosen for its inherent ability to handle constraints on the system’s output, which are applied to restrict the maximum swing angle of the load, guaranteeing the safe flight of the aircraft. An additional contribution of this work is the linear model representation of the quadrotor-payload system used in the design of the model-based predictive control (MPC) strategy; although a modest contribution on its own, it is valuable within the controller design methodology. Moreover, as evidenced in the results section, this strategy achieves a good balance between fast stabilization times and small swing angles.

The article is organized as follows: Section 2 describes the nonlinear model of the quadrotor with the suspended load; in the same section, the linearized model used for the controller design is presented. Section 3 introduces the theoretical basis of the predictive controller. Section 4 shows the numerical results obtained. Section 5 presents the conclusions and future work. Finally, due to the many of variables used, the nomenclature is summarized in Appendix A.

2. Dynamic Model

This section aims to present the dynamic model that describes the complete system, including the coupling problem between the load suspended by means of a cable and the aerial vehicle as an integral problem and not as a disturbance. For this purpose, the following assumptions are made:

(a) 

The multi-rotor fuselage is considered rigid and symmetrical.

(b) 

The cable that connects to the load is attached to the center of mass of the air vehicle.

(c) 

The cable is considered rigid, inelastic, and massless. Its length is constant and is known.

(d) 

The payload is considered to be a point-mass.

(e) 

Aerodynamic effects on the load are neglected.

Consider a quadcopter type multirotor vehicle evolving in three dimensions, capable of carrying a suspended payload by means of a cable as shown in Figure 1. The inertial coordinate system fixed to ground is represented by $\left\{\mathit{O}\right\}$ and the body-fixed frame is defined by $\left\{\mathit{B}\right\}$.

The generalized coordinates of the model can be written as:

$$\begin{array}{c}\hfill q={\left[\begin{array}{ccc}\xi & \eta & \mu \end{array}\right]}^T\in {\mathbb{R}}^8\end{array}$$

(1)

where $\xi ={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T\in {\mathbb{R}}^3$ represents linear positions and $\eta ={\left[\begin{array}{ccc}\psi & \theta & \varphi \end{array}\right]}^T\in {\mathbb{R}}^3$ represents angular positions of the vehicle. $\mu ={\left[\begin{array}{cc}\alpha & \beta \end{array}\right]}^T\in {\mathbb{R}}^2$ represents the oscillation angles of the suspended load. l is the length of the cable, d is the distance between the motors and the center of mass, the gravitational acceleration constant is represented by g, the mass of the quadrotor is defined by M, and the payload is m. The vector of control inputs is represented by $u={\left[\begin{array}{cc}{u}_1& \tau \end{array}\right]}^T\in {\mathbb{R}}^4$, where $u_1=f_{P1}+f_{P2}+f_{P3}+f_{P4}$ is the total propulsive force and $\tau ={\left[\begin{array}{ccc}{\tau }_{\psi }& {\tau }_{\theta }& {\tau }_{\varphi }\end{array}\right]}^T$ are the input torques.

2.1. Euler–Lagrange Methodology

The Euler–Lagrange equations of motion represent an adequate and convenient analytical method to obtain the dynamic model of the system. Their mathematical structure allows one to analyze and study in detail all the physical phenomena of the plant quadrotor-payload, since the dynamic nature of the system is highly nonlinear, underactuated, and multivariable, and the subsystems are strongly coupled.

The Euler–Lagrange formulation is based on the energy of the system; therefore, the kinetic energy has to be calculated as the potential energy to obtain the Lagrangian of the system. The total kinetic energy $T_{QSP}\left(\dot{q}\right)$ of the plant is represented by the following expression:

$$\begin{array}{c}\hfill T_{QSP}\left(\dot{q}\right)=\underset{Quadrotor}{\underbrace{\frac{1}{2}M{\dot{\xi }}^T\dot{\xi }+\frac{1}{2}{\dot{\eta }}^{T}J\dot{\eta }}}+\underset{Payload}{\underbrace{\frac{1}{2}m{\dot{\xi }}_p^T{\dot{\xi }}_p+\frac{1}{2}{I}_p({\dot{\alpha }}^2+{\dot{\beta }}^2)}}\end{array}$$

(2)

where J is the inertia matrix and $I_p$ contains the mass moments of inertia of the payload. The positions relative to the center of the load in Cartesian coordinates are described by:

$$\begin{array}{c}\hfill {\xi }_p=\xi +lr,\end{array}$$

(3)

where ${\xi }_p={\left[\begin{array}{ccc}{x}_p& y_p& z_p\end{array}\right]}^T$ and $r={\left[\begin{array}{ccc}{s}_{\alpha }{c}_{\beta }& s_{\alpha }{s}_{\beta }& -c_{\alpha }\end{array}\right]}^T$, to simplify the presentation, we have used the short notation $s_{\alpha }=sin\left(\alpha \right)$ and $c_{\alpha }=cos\left(\alpha \right)$. The potential energy is due to the position or the configuration of the body with respect to a reference frame, so the total potential energy is described by:

$$\begin{array}{c}\hfill V\left(q\right)=Mgz+mg(z-l{c}_{\alpha })\end{array}$$

(4)

The Lagrangian L, which is what we are interested in calculating, is obtained by subtracting the kinetic energy $T_{QSP}$ from the potential energy $V\left(q\right)$, as shown in the following expression:

$$\begin{array}{c}\hfill L=\frac{1}{2}M{\dot{\xi }}^T\dot{\xi }+\frac{1}{2}{\dot{\eta }}^{T}J\dot{\eta }+\frac{1}{2}m{\dot{\xi }}_p^T{\dot{\xi }}_p+\frac{1}{2}{I}_p({\dot{\alpha }}^2+{\dot{\beta }}^2)-Mgz-mg(z-l{c}_{\alpha })\end{array}$$

(5)

By applying the Euler–Lagrange formulation,

$$\begin{array}{c}\hfill \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}}\right)-\left(\frac{\partial L}{\partial q}\right)=Bu,\end{array}$$

(6)

we obtain the non-linear model of the quadrotor system with a suspended load by means of a cable represented by the following equations:

$$\begin{array}{cc}\hfill \ddot{x}=& u_1((s_{\varphi }{s}_{\psi }+c_{\varphi }{c}_{\psi }{s}_{\theta })-(ml{c}_{\alpha }\ddot{\alpha }+lm{s}_{\alpha }{\dot{\alpha }}^2))/(m+M);\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \ddot{y}=& u_1((c_{\varphi }{s}_{\theta }{s}_{\psi }-c_{\psi }{s}_{\varphi })-(ml{c}_{\beta }\ddot{\beta }+lm{s}_{\beta }{\dot{\beta }}^2))/(m+M);\hfill \end{array}$$

(8)

$$\ddot{z}=(u_1\left(c_{\theta }{c}_{\varphi }\right)-(m+M)g-ml(s_{\alpha }{c}_{\beta }\ddot{\alpha }-ml(s_{\beta }{c}_{\alpha }\ddot{\beta }-ml{\dot{\beta }}^2{c}_{\alpha }-ml{\dot{\beta }}^2{c}_{\beta }+2Ml{s}_{\beta }{s}_{\alpha }\dot{\alpha }\dot{\beta })))/(m+M);$$

(9)

$$\begin{array}{cc}\hfill \ddot{\psi }=& ({\tau }_{\psi }-\ddot{\theta }\left(I_y{I}_z{c}_{\theta }{s}_{\varphi }{c}_{\varphi }{c}_{\theta }{s}_{\varphi }{c}_{\varphi }\right)+I_x\ddot{\varphi }{s}_{\theta }-\dot{\psi }(I_x\dot{\theta }{s}_{\theta }{c}_{\theta }+I_y(-\dot{\theta }{s}_{\theta }{c}_{\theta }{s}_{\varphi }^2+\dot{\varphi }{c}_{\theta }^2{s}_{\varphi }{c}_{\varphi })-I_z(\dot{\theta }{s}_{\theta }{c}_{\theta }{c}_{\varphi }^2+\dot{\varphi }{s}_{\varphi }{c}_{\varphi }{c}_{\theta }^2)\hfill \\ & -\dot{\theta }(I_x\dot{\psi }{s}_{\theta }{c}_{\theta }-I_y(\dot{\theta }{s}_{\theta }{s}_{\varphi }{c}_{\varphi }+\dot{\varphi }{c}_{\theta }{s}_{\varphi }^2-\dot{\varphi }{c}_{\theta }{c}_{\varphi }^2+\dot{\psi }{s}_{\theta }{c}_{\theta }{s}_{\varphi }^2+I_z(\dot{\varphi }{c}_{\theta }{s}_{\varphi }^2-\dot{\varphi }{c}_{\theta }{c}_{\varphi }^2-\dot{\psi }{s}_{\theta }{c}_{\theta }{c}_{\varphi }^2+\dot{\theta }{s}_{\theta }{s}_{\varphi }{c}_{\varphi }))\hfill \\ & +\dot{\varphi }(I_x\dot{\theta }{c}_{\theta }-I_y\dot{\psi }{c}_{\theta }^2{s}_{\varphi }{c}_{\varphi }+I_z\dot{\psi }{c}_{\theta }^2{s}_{\varphi }{c}_{\varphi }))/(I_x{s}_{\theta }^2+I_y{c}_{\theta }^2{s}_{\varphi }^2+I_z{c}_{\theta }^2{c}_{\varphi }^2);\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \ddot{\theta }=& ({\tau }_{\theta }-\ddot{\psi }(I_y{c}_{\theta }{s}_{\varphi }{c}_{\varphi }-I_z{c}_{\theta }{s}_{\varphi }{c}_{\varphi })-\dot{\psi }(-I_x\dot{\psi }{s}_{\theta }{c}_{\theta }+I_y\dot{\psi }{s}_{\theta }{c}_{\theta }{s}_{\varphi }^2+I_z\dot{\psi }{s}_{\theta }{c}_{\theta }{\left(c_{\varphi }\right)}^2)+\dot{\theta }(I_y\dot{\varphi }{s}_{\varphi }{c}_{\varphi }-I_z\dot{\varphi }{s}_{\varphi }{c}_{\varphi })\hfill \\ & -\dot{\varphi }(I_x\dot{\psi }{c}_{\theta }+I_y(-\dot{\theta }{s}_{\varphi }{c}_{\varphi }+\dot{\psi }{c}_{\theta }{c}_{\varphi }^2-\dot{\psi }{c}_{\theta }{s}_{\varphi }^2)+I_z(\dot{\psi }{c}_{\theta }{s}_{\varphi }^2-\dot{\psi }{c}_{\theta }{c}_{\varphi }^2)+\dot{\theta }{s}_{\varphi }{c}_{\varphi })+\dot{\theta }(I_y\dot{\varphi }{s}_{\varphi }{c}_{\varphi }-I_z\dot{\varphi }{s}_{\varphi }{c}_{\varphi }))/(I_y{c}_{\varphi }^2+I_z{s}_{\varphi }^2);\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \ddot{\varphi }=& ({\tau }_{\varphi }+\ddot{\psi }{I}_x{s}_{\theta }+\dot{\psi }(I_y\dot{\psi }{c}_{\theta }^2{s}_{\varphi }{c}_{\varphi }-I_z\dot{\psi }{s}_{\varphi }{c}_{\varphi }{c}_{\theta }^2-\dot{\theta }(-I_x\dot{\psi }{c}_{\theta }+I_y(\dot{\theta }{s}_{\varphi }{c}_{\varphi }+\dot{\psi }{c}_{\theta }{s}_{\varphi }^2-\dot{\psi }{c}_{\theta }{c}_{\varphi }^2)-I_z(\dot{\psi }{c}_{\theta }{s}_{\varphi }^2-\dot{\psi }{c}_{\theta }{c}_{\varphi }^2+\dot{\theta }{s}_{\varphi }{c}_{\varphi })))/\left(I_x\right);\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \ddot{\alpha }=& (-ml{c}_{\alpha }\ddot{x}-ml{c}_{\beta }{s}_{\alpha }\ddot{z}-m{l}^2{c}_{\alpha }{c}_{\beta }{s}_{\alpha }{s}_{\beta }\ddot{\beta }-gml{s}_{\alpha }{c}_{\beta }+m{l}^2{\dot{\alpha }}^2{c}_{\alpha }{s}_{\alpha }-m{l}^2{c}_{\beta }^2{s}_{\alpha }{c}_{\alpha }{\dot{\alpha }}^2\hfill \\ & +2m{l}^2{s}_{\alpha }^2{c}_{\beta }{s}_{\beta }\dot{\alpha }\dot{\beta }-m{l}^2{c}_{\beta }^2{s}_{\alpha }{c}_{\alpha }{\dot{\beta }}^2)/(m{l}^2(c_{\alpha }^2+c_{\alpha }^2{c}_{\beta }^2)+I_p);\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \ddot{\beta }=& (-ml{c}_{\beta }\ddot{y}-ml{c}_{\alpha }{s}_{\beta }\ddot{z}-m{l}^2{c}_{\alpha }{c}_{\beta }{s}_{\alpha }{s}_{\beta }\ddot{\alpha }-gml{s}_{\beta }{c}_{\alpha }+m{l}^2{\dot{\beta }}^2{c}_{\beta }{s}_{\alpha }-m{l}^2{c}_{\alpha }^2{s}_{\beta }{c}_{\beta }{\dot{\beta }}^2\hfill \\ & +2m{l}^2{s}_{\beta }^2{c}_{\alpha }{s}_{\alpha }\dot{\alpha }\dot{\beta }-m{l}^2{c}_{\alpha }^2{s}_{\beta }{c}_{\beta }{\dot{\alpha }}^2)/(m{l}^2(c_{\beta }^2+c_{\beta }^2{c}_{\alpha }^2)+I_p).\hfill \end{array}$$

(14)

In the previous expressions, $I_x$, $I_y$, and $I_z$ represented the inertias of the quadrotor, and $I_p$ represented the inertia of the suspended load.

2.2. Linear Model

In this subsection, the equations of motion of the quadrotor with a suspended load are linearized around the equilibrium points considering the state vector $\chi ={\left[\begin{array}{cccccccccccccccc}x& \dot{x}& y& \dot{y}& z& \dot{z}& \psi & \dot{\psi }& \theta & \dot{\theta }& \varphi & \dot{\varphi }& \alpha & \dot{\alpha }& \beta & \dot{\beta }\end{array}\right]}^T$ and the output vector $y_c={\left[\begin{array}{cccccccc}x& y& z& \psi & \theta & \varphi & \alpha & \beta \end{array}\right]}^T$. In order to determine them, the system is considered to be fixed in an equilibrium coordinate; i.e., the UAV is in a hover position. There is also an assumption of smooth movements, so the rotation and swing angles are relatively small. The translational and rotational system of the quadrotor with suspended load is linearized around the equilibrium point:

$$\begin{array}{c}\hfill {\chi }_{eq}={\left[\begin{array}{cccccccccccccccc}{x}_{eq}& 0& y_{eq}& 0& z_{eq}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]}^T,\end{array}$$

(15)

and for control inputs:

$$\begin{array}{c}{u}_{eq}={\left[\begin{array}{cccc}(m+M)g& 0& 0& 0\end{array}\right]}^T.\end{array}$$

(16)

Then, equations can be represented in state space as follows:

$$\begin{array}{cc}\hfill \dot{\chi }& =A_c\chi +B_{c}u\hfill \\ \hfill y_c& =C_c\chi \hfill \end{array}$$

(17)

By means of the Taylor series expansion and evaluating the equilibrium point, the linear model (17) is obtained with:

$$\begin{array}{cc}\hfill A_c& =\left[\begin{array}{cccccccccccccccc}0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& a_{2,9}& 0& 0& 0& a_{2,13}& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{4,11}& 0& 0& 0& a_{4,15}& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& a_{14,9}& 0& 0& 0& a_{14,13}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{16,11}& 0& 0& 0& a_{16,15}& 0\end{array}\right];\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill B_c& ={\left[\begin{array}{cccccccccccccccc}0& 0& 0& 0& 0& \frac{1}{M+m}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& \frac{1}{I_z}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{I_y}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \frac{1}{I_x}& 0& 0& 0& 0\end{array}\right]}^T;\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill C_c& =\left[\begin{array}{cccccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\end{array}\right];\hfill \end{array}$$

(20)

where $a_{2,9}=g(2l^2{m}^2+M{l}^{2}m+I_{p}m+I_{p}M)/(M+m)(m{l}^2+I_p)$; $a_{4,11}=-g(2l^2{m}^2+M{l}^{2}m+I_{p}m+I_{p}M/(M+m)(m{l}^2+I_p)$, $a_{2,13}=a_{4,15}=g{l}^2{m}^2(2l^2{m}^2+M{l}^{2}m+I_{p}m+I_{p}M)/{(M+m)}^2/{(m{l}^2+I_p)}^2$; $a_{14,9}=-glm(2l^2{m}^2+M{l}^{2}m+I_{p}m+I_{p}M)/(M+m){(m{l}^2+I_p)}^2=a_{14,13}=a_{16,11}$; $a_{16,15}=-glm(2l^2{m}^2+M{l}^{2}m+I_{p}m+I_{p}M)/(M+m){(m{l}^2+I_p)}^2$.

3. Model-Based Predictive Control

Predictive control makes explicit use of a model of the plant to predict the output of the process in future instants, and based on the minimization of a cost function, a strategy is proposed to obtain the optimal control signal. The minimization of the cost function can be obtained explicitly if the model is linear and there are no constraints; otherwise, an iterative optimization method must be used.

3.1. State Space Model and Input Increments

The behavior of the plant considering the system in discrete-time without disturbances and without measurement errors with n states, l outputs, and m inputs can be described by the following model in state space:

$$\begin{array}{cc}\hfill x_{k+1}& =A_d{x}_k+B_d{u}_k\hfill \\ \hfill y_k& =C_d{x}_k\hfill \end{array}$$

(21)

where $x_{k+1}\equiv x(k+1)$, $x_k$∈${\mathbb{R}}^n$ represents the sate vector, $y_k$∈${\mathbb{R}}^l$ the output vector, $u_k$∈${\mathbb{R}}^m$ the input vector, and k the time-step. For convenience, the increment in the control signal is considered an input variable as follows:

$$\begin{array}{c}\hfill \Delta u_k=u_k-u_{k-1}\end{array}$$

(22)

resulting in an augmented plant model:

$$\begin{array}{cc}\hfill {\zeta }_{k+1}& =A_{\zeta }{\zeta }_k+B_{\zeta }\Delta u_k\hfill \\ \hfill y_k& =C_{\zeta }{\zeta }_k\hfill \end{array}$$

(23)

where

$$\begin{array}{c}\hfill {\zeta }_{k+1}=\left[\begin{array}{c}{x}_{k+1}\\ u_k\end{array}\right];A_{\zeta }=\left[\begin{array}{cc}{A}_d& B_d\\ 0& I\end{array}\right];B_{\zeta }=\left[\begin{array}{c}{B}_d\\ I\end{array}\right];C_{\zeta }=\left[\begin{array}{cc}{C}_d& 0\end{array}\right].\end{array}$$

3.2. Predictions

The future evolution of the states can be calculated by iterating over the model (23):

$$\begin{array}{c}\hfill \underset{{\underrightarrow{\zeta }}_k}{\underbrace{\left[\begin{array}{c}{\zeta }_{k+1}\\ {\zeta }_{k+2}\\ {\zeta }_{k+3}\\ ⋮\\ {\zeta }_{k+n_y}\end{array}\right]}}=\underset{P_{\zeta \zeta }}{\underbrace{\left[\begin{array}{c}{A}_{\zeta }\\ A_{\zeta }^2\\ A_{\zeta }^3\\ ⋮\\ A_{\zeta }^{n_y}\end{array}\right]}}{\zeta }_k+\underset{P_{\zeta \Delta u}}{\underbrace{\left[\begin{array}{cccc}{B}_{\zeta }& 0& \cdots & 0\\ A_{\zeta }{B}_{\zeta }& B_{\zeta }& \cdots & 0\\ A_{\zeta }^2{B}_{\zeta }& A_{\zeta }{B}_{\zeta }& \cdots & B_{\zeta }\\ ⋮& ⋮& \ddots & ⋮\\ A_{\zeta }^{n_y-1}{B}_{\zeta }& A_{\zeta }^{n_y-2}{B}_{\zeta }& \cdots & A_{\zeta }^{n_y-n_u}{B}_{\zeta }\end{array}\right]}}\underset{{\underrightarrow{\Delta u}}_{k-1}}{\underbrace{\left[\begin{array}{c}\Delta u_k\\ \Delta u_{k+1}\\ ⋮\\ \Delta u_{k+n_u-1}\end{array}\right]}}\end{array}$$

(24)

where the arrow notation denotes prediction and is defined as ${\underrightarrow{x}}_k={\left[x_{k+1}^T{x}_{k+2}^T\dots \right]}^T$. $n_y$ is the prediction horizon and $n_u$ the control horizon—it is generally chosen by $n_u\le n_y$. The state predictions can be obtained by extracting them from the augmented states:

$$\begin{array}{cc}\hfill {\underrightarrow{x}}_k& =\underset{P_{x\zeta }}{\underbrace{diag\left(\left[\begin{array}{cc}I& 0\end{array}\right]\right)}}{\underrightarrow{\zeta }}_k.\hfill \end{array}$$

(25)

By substituting the state prediction into the output equation, the output predictions can be obtained:

$$\begin{array}{cc}\hfill {\underrightarrow{y}}_k& =P_{y\zeta }{\zeta }_k+P_{y\Delta u}{\underrightarrow{\Delta u}}_{k-1},\hfill \end{array}$$

(26)

where $P_{y\zeta }=diag\left(C\right)P_{\zeta \zeta }$ and $P_{y\Delta u}=diag\left(C\right)P_{\zeta \Delta u}$.

Finally, the input predictions are

$$\begin{array}{cc}\hfill {\underrightarrow{u}}_{k-1}& =\underset{P_{u\Delta u}}{\underbrace{\left[\begin{array}{cccc}I& 0& \cdots & 0\\ I& I& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ I& I& \cdots & I\end{array}\right]}}{\underrightarrow{\Delta u}}_{k-1}+\underset{P_{u\zeta }}{\underbrace{col\left(\left[\begin{array}{cc}0& I\end{array}\right]\right)}}{\zeta }_k.\hfill \end{array}$$

(27)

Our paper adopts an independent model approach to prediction [23,24]. Then, by defining ${\widehat{y}}_k$ as the output vector of the independent model (given by simulating model (21) in parallel with the plant), the model’s residual estimate is ${\widehat{d}}_k=y_k-{\widehat{y}}_k$ and it is assumed to be constant over the prediction horizon. This value is subtracted from the reference to compensate for disturbances or model mismatches.

3.3. Cost Function

The most common cost function for obtaining the control law consists of penalizing the squared errors multiplied by a weight, that is, penalizing the difference between the reference minus the squared output plus the control effort multiplied by a weight, along a prediction and control horizon. The general expression of this function is:

$$\begin{array}{c}{J}_{MPC}=\sum _{i=1}^{n_y}{\left(y_{k+i}-\widehat{r}\right)}^{T}Q\left(y_{k+i}-\widehat{r}\right)+\sum _{i=0}^{n_u-1}\Delta u_{k+i}^{T}R\Delta u_{k+i}\end{array}$$

(28)

Q and R are tuning parameters that impact future behavior and $\widehat{r}$ is the setpoint. By substituting the model of the predictions within the objective function by differentiating and setting it to zero, the optimal control sequence is obtained:

$$\begin{array}{c}\hfill {\underrightarrow{\Delta u}}_{k-1}^{*}=-{\left(P_{y\Delta u}^{T}diag\left(Q\right)P_{y\Delta u}+diag\left(R\right)\right)}^{-1}{P}_{y\Delta u}^{T}diag\left(Q\right)\left(P_{y\zeta }{\zeta }_k-\underrightarrow{\widehat{r}}\right)\end{array}$$

(29)

where * denotes optimality. Additionally, the control input $u_k$ is formed by adding to the previous input $u_{k-1}$, the first element of ${\underrightarrow{\Delta u}}_{k-1}^{*}$—that is, $\Delta u_{k-1}^{*}$.

3.4. The Constrained MPC Algorithm

Constraints help to keep the process within desired operating limits or conditions. Constraints have to be fulfilled for all the prediction horizon:

$$\begin{array}{cc}\hfill diag\left(A_x\right){\underrightarrow{x}}_k\le col\left(b_x\right);& diag\left(A_y\right){\underrightarrow{y}}_k\le col\left(b_y\right);\hfill \\ \hfill diag\left(A_u\right){\underrightarrow{u}}_{k-1}\le col\left(b_u\right);& diag\left(A_{\Delta u}\right){\underrightarrow{\Delta u}}_{k-1}\le col\left(b_{\Delta u}\right).\hfill \end{array}$$

(30)

where $A_x,A_y,A_u,A_{\Delta u}={[I,-I]}^T$ and $b_x,b_y,b_u,b_{\Delta u}$ are vectors that contain the maximum and minimum values allowed in the form ${\left[maxmin\right]}^T$.

The cost function (28) and constraints (30) can be expressed in terms of the desicion variable ${\underrightarrow{\Delta u}}_{k-1}$, and the following quadratic programming optimization problem is obtained:

$$\begin{array}{cc}\hfill {\underrightarrow{\Delta u}}_{k-1}^{*}=& \underset{{\underrightarrow{\Delta u}}_{k-1}}{\arg \min }\left\{\frac{1}{2}{\underrightarrow{\Delta u}}_{k-1}^T\left(P_{y\Delta u}^{T}diag\left(Q\right)P_{y\Delta u}+diag\left(R\right)\right){\underrightarrow{\Delta u}}_{k-1}+{\left(P_{y\Delta u}diag\left(Q\right)(P_{y\zeta }{\zeta }_k-{\underrightarrow{\widehat{r}}}_k)\right)}^T{\underrightarrow{\Delta u}}_{k-1}\right\}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & \mathrm{s}.\mathrm{t}.\underset{M_a}{\underbrace{\left[\begin{array}{c}diag\left(A_x\right)P_{x\zeta }{P}_{\zeta \Delta u}\\ diag\left(A_y\right)P_{\zeta \Delta u}\\ diag\left(A_u\right)P_{u\Delta u}\\ diag\left(A_{\Delta u}\right)\end{array}\right]}}{\underrightarrow{\Delta u}}_{k-1}\le \underset{q_a\left({\zeta }_k\right)}{\underbrace{\left[\begin{array}{c}col\left(b_x\right)\\ col\left(b_y\right)\\ col\left(b_u\right)\\ col\left(b_{\Delta u}\right)\end{array}\right]-\left[\begin{array}{c}diag\left(A_x\right)P_{x\zeta }{P}_{\zeta \zeta }\\ diag\left(A_y\right)P_{y\zeta }\\ diag\left(A_u\right)P_{u\zeta }\\ 0\end{array}\right]}}\hfill \end{array}$$

(32)

and the control law is:

$$\begin{array}{c}\hfill u_k=u_{k-1}+\Delta u_k^{*},\end{array}$$

(33)

where $\Delta u_k^{*}$ is the first element of ${\underrightarrow{\Delta u}}_{k-1}^{*}$.

3.5. MPC for a Quadrotor with a Suspended Load

This subsection presents the design and configuration with to implement the predictive control algorithm for a quadrotor carrying a suspended load. We start from a continuous state space model representation (17) to obtain a discretized and augmented version of the system (23) represented by the matrices $A_{\zeta }\in {\mathbb{R}}^{20x20}$, $B_{\zeta }\in {\mathbb{R}}^{20x4}$ and $C_{\zeta }\in {\mathbb{R}}^{8x20}$ that are used to calculate the prediction matrices $P_{\zeta \zeta },P_{\zeta \Delta u},P_{y\zeta },P_{y\Delta u},P_{x,\zeta },P_{u,\Delta u},P_{u,\zeta }$ to form the unconstrained control law (29) and the optimization problem (31) and (32) for the constrained case. To solve the optimization problem, Matlab’s $quadprog$ function was used. It is worth mentioning that for this particular process, only constraints on the positions of the quadrotor x, y, and z, and on the oscillation angles of the suspended load $\alpha$ and $\beta$ were considered.

The physical parameters of the suspended-load quadrotor system are presented in

Table 1. System parameters.

Parameter	Value	Units
M	$0.4$	m
m	$0.03$	m
d	$0.1$	m
l	$0.35$	m
g	$9.8$	m/s${}^2$
$I_x$	$1.77\times {10}^{-3}$	kgm${}^2$
$I_y$	$1.77\times {10}^{-3}$	kgm${}^2$
$I_z$	$3.54\times {10}^{-3}$	kgm${}^2$
$I_p$	$1.00\times {10}^{-6}$	kgm${}^2$

and were taken from [8]. Predictive control tuning parameters are the prediction $n_y$ and control $n_u$ horizons, and the weights Q and R. To select the horizons, the tuning guidelines provided by [25] were used. In general, for stable open-loop models, the literature suggests that $n_y$ should be larger than the settling time plus $n_u$. For unstable open-loop models, as this is the case, large values of $n_y$ are preferred, and then the value is detuned until the system achieves the desired robustness. Here, intermediate values of $n_y$ result in an adequate performance. The control horizon $n_u$ was chosen to be as small as possible to get the required behavior while minimizing computational load (although a higher value of $n_u$ should improve performance if the processing capabilities allow it). Finally, Q and R were heuristically adjusted—a compromise between response time and input rate (aggressiveness). Alternatively, a global search could be done to achieve better performance, as suggested by [26]. The tuning parameters of the predictive controller are presented in

Table 2. Controller parameters.

Parameter	Value	Units
$T_s$	$0.1$	s
$n_y$	10	–
$n_u$	3	–
Q	I	–
R	$0.00001\times I$	–

.

4. Numerical Simulations and Results

This section presents the numerical results of the simulations performed to validate the achievement of the designed MPC controller. The main goal was to bring the quadcopter to the desired reference position while the suspended load oscillation was restricted to a maximum swing of $\pm {10}^{\circ }=\pm 0.175 \mathrm{rad}$. For the model’s linearization, the values of $M,m,$ and g were constant and are defined in Table 1. The values of $x_{eq},y_{eq},z_{eq}$ were zero, which is an arbitrary point in three-dimensional space located well above the ground.

4.1. MPC Performance

In transport operations with multirotors, rapid handling of loads is required, which can be dangerous because the load can generate unwanted oscillations that can vary the flight conditions—the process must be kept within limits. In this particular process, the use of constraints for positions x and y of the vehicle and in the oscillation angles of the load were considered to avoid overshoots that can cause damage to the transported payload and guarantee safety in operations. The results that correspond to the constrained predictive controller are presented in Figure 2.

Figure 2a displays the positions x, y, and z of the quadrotor with their respective constraints $0\le x\le 1$, $0\le y\le 2$, and $0\le z\le 3$. The signals responded quickly to reach the reference in no longer than 3 s. The controller kept each constrained position x and y at its reference value without any overshoots in the transient responses to avoid further oscillations of the load angles.

In Figure 2b the attitude of the quadrotor is shown. No constraint was applied; however, lower values in $\theta$ and $\varphi$ angles occurred in comparison to the unconstrained predictive control.

Figure 2c presents the swing angles of the suspended load with their respective constraints, for $-0.1047\le \alpha \le 0.1047$ and $-0.1745\le \beta \le 0.1745$. As can be seen, the predictive controller kept the oscillation within the established values, demonstrating good performance.

Figure 2d shows the variations of the control inputs in the quadrotor, which ensured that the oscillations of the suspended load did not exceed the values established in the constraints and that the air vehicle could stabilize quickly. In the bottom plot of this figure, we can see the force $u_1$ in the z axis that stabilized at $(m+M)g$ maintain hovering, and in the top plot the control torques ${\tau }_{\psi }$, ${\tau }_{\theta }$, and ${\tau }_{\varphi }$.

In summary, numerical simulations have been presented that demonstrate good performance by the predictive controller when considering the constraints to attenuate the oscillations of the suspended load. It is possible to place the quadrotor at a desired position while the oscillations of the load do not exceed the proposed limits. Next, the results of not considering the constraints are shown, for comparative purposes.

4.2. Comparative Results

In Figure 3, a direct comparison of the performance of the proposed MPC algorithm with respect to a nonlinear controller based on the feedback linearization proposed in Equation (16) of [22] has been added. The control law in [22] was chosen because it manages to achieve small swing angles of the payload.

For a better understanding of the figures, know that the subscript “c” means the constrained MPC controller and the subscript “l” the algorithm in [22]. In order to contrast the results with a more aggressive strategy, we also added the unconstrained MPC control law, which is referred to with the subscript “u”. The performance responses of all controllers are displayed in Figure 3. In Figure 3c we can observe the swing angles. It is clear that [22] achieved smaller maximum angles, but the oscillation times for $\alpha$ and $\beta$ were longer than in the MPC strategy. Moreover, by analyzing Figure 3a, it is possible to appreciate better performance in the vehicle position response of the proposed constrained MPC controller with respect to [22], because the settling time of the vehicle position obtained from control law $u_l$ was larger than that obtained with the MPC law. In summary, ref. [22] achieveed small swing angles of the payload but at the cost of poor performance in the vehicle positioning and settling time of the payload swing. In contrast, this proposal produced shorter settling times for the payload and vehicle positioning while still meeting the constraints imposed on the payload swing angles.

In order to quantify the comparisons,

Table 3. Performances of the proposed control and the controller in [22].

Controller	Cost ${\mathit{J}}_{\mathit{MPC}}$	Swing Angle	Ts (s)	MS (Rad)
$u_c$	624.8	$\alpha$	3.5	0.1
		$\beta$	3.5	0.175
$u_u$	985.1	$\alpha$	3.5	0.28
		$\beta$	3.5	0.5
$u_l$	1997.3	$\alpha$	7.2	0.05
		$\beta$	7.2	0.15

summarizes the results of the proposed controller and those of the algorithm in [22] with respect to performance measures and oscillation angles. We computed for $\alpha$ and $\beta$, the setting time (Ts) and the maximum swing (MS). The superior performance of the proposed scheme over the one in [22] is evidenced by the running-time cost $J_{MPC}$.

5. Conclusions

This work presents the formulation of a model-based predictive controller meant to solve the problem of transporting a payload suspended by a cable from a multi-rotor-type unmanned aerial vehicle. This strategy considers the three-dimensional translation and rotation of the aerial vehicle, and the payload dynamics. The dynamic model was obtained by formulating a Euler–Lagrange in a global approach that mathematically describes the relationship of the vehicle, the load, and the cable. Next, an affine representation of the input of the model was derived to obtain a linear model for the MPC design. It is worth mentioning that the steps carried out to obtain the linear model of the quadrotor-payload are specific for this system, and the linear model was used as a prediction model for the control strategy. Finally, we evaluated the dynamic model of the UAV with a suspended load, and the proposed controller, using a numerical simulation software. The simulations showed excellent performance by the MPC in comparison to a nonlinear controller based on feedback linearization from the literature. Therefore, the proposed controller demonstrated its ability to precisely control the multi-rotor and achieve rapid attenuation of the oscillation of the suspended load while keeping it within the limits of the constraints. Future work will include considering a flexible cable in the mathematical model and an MPC with variable constraints.