Control structure design with constraints for a slung load quadrotor system

Abstract

We propose a control structure for a quadrotor carrying a slung load with swing-angle constraints. This quadrotor is supposed to pass through the waypoints at specified speeds. First, a cascaded PID autopilot is designed, which adaptively gives attention to position and speed requirements as a function of their errors. Its parameters are found from an optimization problem solved using the PSO algorithm. Second, this controller’s performance is improved by adding the Complementary Controller employing an ANN. 5. Training data for the ANN is created by solving optimal control problems. The ANN is activated when the swing angle constraint is about to be violated. It is trained using optimal control values corresponding to the cases where the swing angle falls in a particular band about the upper swing angle constraint. Simulations are performed in a MATLAB environment. Finally, some of the simulation results are validated on a physical system.

Introduction

In recent years, with the technological developments in unmanned aerial vehicles (UAV), studies on package delivery, ¹ agricultural applications, ² search and rescue in dangerous areas,^3,4 and information gathering from desired regions,^5–7 and military applications have gained momentum. ⁸ Load transportation is vital in all these applications. ⁸ Two different methods are being studied for carrying loads with rotary-wing UAVs. One can carry a load using a manipulator^9–12 integrated into the vehicle. The other is to connect the load to the vehicle with a cable or rope.^13–15 In vehicles with a manipulator, the vehicle’s inertia grows due to components associated with the manipulator – consequently, energy utilization of the vehicle with the manipulator increases. One of the biggest problems for UAVs is the short mission time. In addition, especially when it is required to perform agile movements, the growth in inertia will oppose the maneuvering. For these reasons, hanging loads seem more advantageous. However, using the suspended load method is also challenging because system dynamics change during the lift, hover, and flight. Even though we have not covered it in this study, we would like to emphasize that one of the important issues about load-bearing quadrotors is precise landing. This topic is covered in depth in Xuan-Mung et al. ¹⁶ and Xuan Mung et al. ¹⁷

As Khamseh et al. ⁴ highlighted, a complete mathematical model covering all phases of flight is the biggest challenge. Models are created assuming the cable or rope is constantly tense during flight. Flight profiles that violate these assumptions should be prevented during flight (e.g. the vehicle descends quickly). Usually, there are swings during the flight due to the quadrotor’s acceleration. These motions have distorting effects on the quadrotor’s dynamics.

A controller should be designed to make the vehicle follow the desired path and minimize these swinging effects. The controller designed in this study is carried out based on these considerations.

In the literature, extensive studies have been published to model the suspended load-bearing rotary-wing UAV systems. Modeling was carried out using the Udwadia-Kalaba Equations in Goodarzi et al. ¹⁸ and Almeida and Raffo, ¹⁹ the Euler-Lagrange approach in Lee et al. ²⁰ and Lee et al., ²¹ the Newton-Euler method in Jain, ²² and the Kane method in Klausen et al. ²³ In Pizetta et al., ²⁴ the system dynamics are defined as the second-order Euler-Lagrange model while ignoring the aerodynamic effects. In Barikbin and Fakharian, ²⁵ the Euler-Lagrange method was used, but the load’s motion was assumed to be planar rather than 3-dimensional. In this study, the controller is designed using the mathematical model detailed in Lee ²⁶ based on the Euler-Lagrange approach.

Measuring and estimating suspended load oscillations is critical for controlling the system. The oscillation frequency caused by the suspended load has been investigated by Geronel et al., ²⁷ in which the dynamic model is constructed by obtaining dimensionless equations of motion. A method is proposed to determine the vibration frequency of the payload during flight. In recent years, studies have been focused on the controller design for such systems to reduce the oscillations of the suspended load. In Sadr et al., ²⁸ a nonlinear control algorithm is developed that applies a pulse to the system to prevent oscillation, whose amplitude and time locations are calculated by estimating the system’s natural frequency and damping ratio. Due to the quadrotor’s movements, the slung load will start to swing while moving to the desired point. This swinging action affects the quadrotor’s traveling motion. ²⁹ The suspended load can create unstable oscillations due to its disturbing dynamics and reduce the performance of the UAV. They may even cause crashes in some cases. ³⁰ Therefore, it has been studied for a long time to solve the difficulties of UAVs in carrying suspended loads. Controller and path planning is the most critical parts of these studies. In Jain, ²² a controller is designed for a quadcopter carrying a point-mass load with a rigid cable. Predictive control is considered for carrying the payload from point to point with swing angle constraints along the desired trajectory.^23,24 Guerrero et al. ³¹ present a passivity based-controller design to minimize the load swing. However, a simplified 2D model is used in this study. In Kusznir and Smoczek, ³² sliding mode control is used for horizontal positioning and payload vibration damping, while a linearizing feedback controller is used for altitude and attitude control. In Yu et al., ³³ the backstepping procedure is applied to design a controller for the load to follow a pre-determined trajectory by defining a virtual thrust force. They find the actual forces using this virtual force. Another controller based on the backstepping method is proposed by Dai et al. ³⁴ to reduce the swing angle of the load. There, the swing angle of the suspended payload is estimated using the harmonically extended state observer. In Shi et al., ³⁵ an adaptive control scheme with unknown mass transportation is proposed. Lee et al. ³⁶ proposes an adaptive sliding mode controller to achieve the goal of altitude tracking control in the presence of considerable system parameter variations and the ground effect.

After choosing the controller type, its parameters should be determined to optimize the performance. Usually, evolutionary algorithms are preferred for this purpose; PSO (The PSO technique is used to design the PID controller to stabilize the humanoid robot,^37,38) Ant Colony Optimization (ACO) technique, ³⁹ Genetic Algorithm (GA), ⁴⁰ Differential Evolution (DE), ⁴¹ (the pitch angle controller of a rocket system is optimized algorithm). Comparisons made by El-Ghandour and Elbeltagi ⁴² have demonstrated through in-depth studies that PSO is more efficient and has lower computational complexity simultaneously.

There are many reasons why Artificial Neural Networks are used in control applications. First, the self-learning ability of neural networks eliminates the use of complex and challenging mathematical analysis that predominates in many adaptive and optimal control methods. Second, it provides a nonlinear mapping capability to solve highly nonlinear control problems by incorporating the Activation Function into hidden neurons. Third, unmodelled dynamics and perturbations in the system are expected to be dealt with using an ANN. Thus, neural controllers can be applied over a broader range of uncertainty.

As ANN has proven to have excellent approximation and learning capabilities, they have become popular tools in control applications of nonlinear systems.^43,44 A unified framework for describing and controlling nonlinear dynamical systems is proposed in Polycarpou. ⁴⁵ The adaptive nonlinear control and the parametric method of adaptive linear control theory can be applied to perform stability analysis. Typical stable neural network approach control schemes based on Lyapunov training designs are given in Mughees and Mohsin. ³⁹ Training neural networks using a backpropagation algorithm (BP) has dramatically increased the development of these structures in control applications. ⁴⁶ Many research studies have proven the approximation capabilities of neural networks. ⁴⁷ In this context, many neural control approaches have been developed with BP neural networks.^45,48,49 Moreover, the ANN structure has been used in different control applications; Guo et al. ⁵⁰ provides a decentralized robust optimal monitoring control scheme for the class of nonlinear systems interconnected with disturbances and actuator attacks, for which the optimal control problem is solved by NN-based adaptive dynamic programing.

We aim to make a quadrotor carry a suspended load following a given path without the load swinging too much. Here, the slung-load system is supposed to follow a path using the generalized waypoint guidance algorithm (waypoints and desired speeds are specified). The original contributions of this study are described below:

1. An autopilot is designed consisting of cascaded PID controllers for inner and outer loops. This autopilot is supposed to control the position and velocity of the quadrotor and prevent high swing angles of the suspended load.

2. The primary controller adaptively changes (the percentage of) its attention to speed or position depending on its errors. While the PID coefficients are found, the controller gains associated with positions and velocities are adaptively tuned by gain adjustment (see equations (8)–(11)).

3. Finally, the whole structure is fine-tuned using delta force and torques produced by a Complementary Controller with an ANN to attenuate the oscillations if there is too much swinging.

4. Training data for the ANN is created by solving optimal control problems to calculate the attenuating delta force and torque inputs to the plant when the swing angle is above a certain threshold. These calculations are performed while the vehicle tries to follow specially chosen paths.

To the best of our knowledge, there is no such study for a slung-load quadrotor system.

The rest of the paper is organized as follows; Firstly, the dynamic model of the system is given in Section II. Secondly, the structure and design procedure of the controller are described in Section III. Section IV contains the experimental results. Then, this study is concluded in Section V.

Dynamic model

This section will describe the 3D dynamic model for the slung-load quadrotor system. It starts with the coordinate system definitions in the mathematical model, as shown in Figure 1. The mathematical model used in this study is borrowed from Lee. ²⁶ The detailed derivation of this model can be found there.

Figure 1.

Slung-load quadrotor system coordinate frames.

This model is differentially flat; the flat outputs are the load’s position and the quadrotor’s yaw angle. Unlike Lee, ²⁶ the 3D dynamic model is preferred over the flat model since we imposed constraints on the quadrotor’s position and attitude. Additionally, different from previous studies,^46,51,52 we do trajectory planning for both the load and the quadrotor to avoid obstacles. The state vector is, $\mathit{x}\left(t\right)={\left[{\mathit{x}}_{\mathit{L}}^{\mathit{T}}{\mathit{p}}^{\mathit{T}}\mathit{R}{\stackrel{·}{\mathit{x}}}_{\mathit{L}}^{\mathit{T}}{\stackrel{·}{\mathit{p}}}^{\mathit{T}}{\mathit{\Omega }}^{\mathit{T}}\right]}^{\mathit{T}}\in {\mathbb{R}}^{18}$ where ${\mathit{x}}_{\mathit{L}}$ is the position vector of the load, and $\mathit{p}$ is the unit vector from the quadrotor to the load in the inertial frame. The input vector is,

$\mathit{v}={\left[f{M}_{\varphi }{M}_{\theta }{M}_{\psi }\right]}^T$ where $f=\mathrm{T}1+\mathrm{T}2+\mathrm{T}3+\mathrm{T}4$ is the lifting force acting on the system $\mathbf{M}={\left[M_{\varphi }{M}_{\theta }{M}_{\psi }\right]}^T$ are roll, pitch, and yaw moments (see Figure 1).

The dynamic model for the slung-load quadrotor system is valid under the following assumptions:

• The payload is connected with a cable at the quadrotor’s center of gravity.

• The cable is rigid, massless, and inelastic.

• The quadrotor has no external disturbances or interactions, such as wind, except for the load.

The input vector $\mathit{v}$ is found as;

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{v}=\left[\begin{array}{c}\begin{array}{c}f\hfill \\ \multicolumn{1}{c}{M_{\varphi }}\\ \multicolumn{1}{c}{M_{\theta }}\end{array}\hfill \\ \multicolumn{1}{c}{M_{\psi }}\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{K}_f\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{K_f}\end{array}\hfill \\ \multicolumn{1}{c}{K_M}\end{array}\begin{array}{c}\begin{array}{c}{K}_f\hfill \\ \multicolumn{1}{c}{-K_f}\\ \multicolumn{1}{c}{0}\end{array}\hfill \\ \multicolumn{1}{c}{-K_M}\end{array}\begin{array}{c}\begin{array}{c}{K}_f\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-K_f}\end{array}\hfill \\ \multicolumn{1}{c}{K_M}\end{array}\begin{array}{c}\begin{array}{c}{K}_f\hfill \\ \multicolumn{1}{c}{K_f}\\ \multicolumn{1}{c}{0}\end{array}\hfill \\ \multicolumn{1}{c}{-K_M}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\omega }_1^2\hfill \\ \multicolumn{1}{c}{{\omega }_2^2}\\ \multicolumn{1}{c}{{\omega }_3^2}\end{array}\hfill \\ \multicolumn{1}{c}{{\omega }_4^2}\end{array}\right]}\end{array}$$

(1)

where ${\omega }_j,(j=1,2,3,4)$ is the angular velocity of rotor $j$ , $K_f$ aerodynamic force constant, $K_M$ aerodynamic moment constant. The mathematical model equations are given below:

$$\frac{d}{\mathrm{dt}}{\mathit{x}}_{\mathit{L}}={\stackrel{·}{\mathit{x}}}_{\mathit{L}}$$

(2)

$$\frac{d}{\mathrm{dt}}{\stackrel{·}{x}}_L=\frac{1}{m_q+m_l}(\mathit{p}.f\mathit{R}{\mathit{e}}_{\mathit{z}}-m_{q}l\left|\left|\stackrel{·}{\mathit{p}}\right|\right|)\mathbf{p}-g{\mathbf{e}}_{\mathbf{z}}$$

(3)

$$\frac{d}{\mathrm{dt}}\mathbf{p}=\stackrel{·}{\mathit{p}}$$

(4)

$$\frac{d}{\mathrm{dt}}\stackrel{·}{\mathit{p}}=\frac{1}{m_{q}l}(\mathbf{p}\times (\mathit{p}\times f\mathit{R}{\mathit{e}}_{\mathit{z}})-\left|\left|\stackrel{·}{\mathit{p}}\right|\right|\mathbf{p}$$

(5)

$$\frac{d}{dt}\mathbf{R}=\mathbf{R}\widehat{\mathbf{\Omega }}$$

(6)

$$\frac{d}{\mathrm{dt}}\mathbf{\Omega }={\mathbf{I}}^{-1}(\mathbf{M}-\mathbf{\Omega }\times \mathbf{I}\mathbf{\Omega })$$

(7)

where R is the rotation matrix from Body Frame to Inertial Frame, ${\mathbf{e}}_{\mathbf{z}}$ is the unit vector along the reference z-direction. $\mathbf{\Omega }$ is the angular velocity of the quadrotor. I is the inertia matrix, and $m_q$ and $m_l$ are the mass of the quadrotor and the load, respectively. $\hat{\mathbf{\Omega }}$ is the skew-symmetric matrix described in Lee. ⁵³

Controller design

Our study gives the waypoints the unmanned aerial vehicle should visit and the desired speed values at these points. When the system reaches a certain proximity to a desired waypoint, this point is considered already visited, and a command is produced to go to the next waypoint. Figure 2 shows this structure for two consecutive waypoints. In our study, as in the trajectory tracking approach, the entire path is not given, but only specific waypoints (the dashed lines between the waypoints are drawn to show the approach easily).

Figure 2.

Waypoint geometry used in controller design.

The slung-load quadrotor system is an under-actuated system having 8 degrees of freedom (DOF) and 4 control surfaces. It is a challenging task to control this system with velocity constraints. The controller structure shown in Figure 3 has been created to overcome this challenge. In this figure, each controller block consists of a PID controller. The vehicle’s roll and pitch angles must be changed for position and speed control. Since we want to perform speed and position control, another block for gain adjustment has been added. The gain adjustment block adjusts effectiveness by changing the roll or pitch angles according to speed and position errors. This effect is achieved by calculating the $g1$ and $g2$ parameters at each time step. They are set between 0.2 and 0.8 according to the error in speed and the angle between the current position and the target waypoint adaptively (equations (8)–(11)). Figure 3 shows the controller structure used for position and speed control. This structure does not have a complementary controller, and the swing angle constraint is not expected to be satisfied in scenarios not used to determine PID gains. The gain adjustment parameters $g1$ and $g2$ are calculated as follows:

$${\mathrm{\Delta }}_{\mathrm{direction}}=\mathrm{abs}\left(\frac{\arctan \left(\frac{w{p}_y-c{p}_y}{w{p}_x-c{p}_x}\right)}{\pi }\right)$$

(8)

$${\mathrm{\Delta }}_{\mathrm{speed}}=\frac{(V_{\mathrm{desired}}-V_{\mathrm{measured}})}{V_{\max }}$$

(9)

$$g1=\min (max(\frac{{\mathrm{\Delta }}_{\mathrm{direction}}}{{\mathrm{\Delta }}_{\mathrm{direction}}+{\mathrm{\Delta }}_{\mathrm{speed}}},0.001),1)$$

(10)

$$g2=1-g1$$

(11)

where $V_{\max }$ is the maximum allowable quadrotor’s speed, $w{p}_x\mathrm{and}w{p}_y$ are the $x\mathrm{and}y$ coordinate of the target waypoint, $c{p}_x\mathrm{and}c{p}_y$ are the $x\mathrm{and}y$ coordinate of the quadrotor’s current position. ${\mathrm{\Delta }}_{\mathrm{direction}}$ is the normalized angle between the quadrotor’s current position and the target waypoint. ${\mathrm{\Delta }}_{\mathrm{speed}}$ is the normalized speed error. ${\mathrm{\Delta }}_{\mathrm{direction}}$ and ${\mathrm{\Delta }}_{\mathrm{speed}}$ are saturated between 0.001 and 1 to ensure not to get zero division error and to compensate for round-off errors.

Figure 3.

Controller structure without the complementary controller.

Tuning the PID coefficients has been proposed as an optimization problem solved using the PSO algorithm. The problem definition is given below in detail.

Definition of the optimization problem

The control objective is to ensure that all waypoints are visited in a certain order and at the desired speed at the waypoint in minimum time. The system must visit the waypoints sequentially. A proximity radius is defined for waypoints considered visited when the distance to a waypoint less than the proximity radius is. The proximity radius may change according to the task performed, which is considered one of the problem’s input parameters.

$$\begin{array}{l}\underset{K_1,\dots ,K_{27}}{\min }J\left(K_1,\dots ,K_{27}\right)=\\ ({\alpha }_1\sum _j^{N}perro{r}_j^2+{\alpha }_2\sum _j^{N}verro{r}_j^2)+{\alpha }_3{T}_{end}\\ s.t.\\ \left|{\beta }_i\right|\le {\beta }_{max}i=0,\mathrm{\Delta }t,2\mathrm{\Delta }t\dots T_{end}\\ Eq.2to7\end{array}$$

(12)

where $K_1,\dots ,K_{27}$ are the PID coefficients. For ease of expression, they are enumerated starting from the proportional gain of the PID1. (i.e. $K_1$ is the proportional gain of the PID1, $K_2$ is the integral gain of the PID1, and $K_{27}$ is the derivative gain of the PID9). ${\alpha }_1,{\alpha }_2,\mathrm{and}{\alpha }_3$ values are appropriately selected, they create a convex sum. ${\beta }_i$ is the swing angle at ith time-step, and ${\beta }_{\max }$ is the maximum allowable swing angle. $T_{\mathrm{end}}$ is the duration of the mission. $\mathrm{\Delta }t$ is the simulation time step ( $\mathrm{\Delta }t=1\mathrm{msec}).$ The connection between optimization problem (also Algorithm-1) and the model is that the system’s dynamical model is taken as the constraint of the optimization problem given in equation (12).

Details of the objective function calculation are given below.

$$\begin{array}{c}\mathrm{perro}{r}_j=\{\begin{array}{c}0j=\mathrm{next}\mathrm{wp}\mathrm{AND}{d}_j\le r_{\mathrm{proximity}}\hfill \\ \multicolumn{1}{c}{d_{j}j=\mathrm{next}\mathrm{wp}\mathrm{AND}{d}_j>r_{\mathrm{proximity}}}\\ \multicolumn{1}{c}{d_{\max }j\ne \mathrm{next}\mathrm{wp}}\end{array}\hfill \\ j=1,2,\dots N\hfill \end{array}$$

(13)

where, $r_{\mathrm{proximity}}$ is the distance required for the waypoint to be considered visited. It can be determined according to the mission to be performed. As stated before, the system must visit the waypoints in order. Considering this constraint, $\mathrm{perro}{r}_j$ is calculated. This calculation includes the system’s distance to the waypoints at the time of the visit. $d_j$ is the distance between the jth waypoint and the quadrotor’s position. $d_{\max }$ is the furthest distance that the system can travel. For waypoints that have not yet been visited and are not the next waypoint in the visiting queue, $\mathrm{perro}{r}_j$ is taken as $d_{\max }$ .

$$\begin{array}{c}\mathrm{verro}{r}_j=\{\begin{array}{c}\left||v_d-v_a|\right|d_j\le r_{\mathrm{proximity}}\hfill \\ \multicolumn{1}{c}{V_{\max }{d}_j>r_{\mathrm{proximity}}}\end{array}\hfill \end{array}$$

(14)

Since the speed objective is valid only within the radius of the proximity area, the function is written in parts, and the velocity error outside this region is considered the maximum allowable velocity value.

PID gains using particle swarm optimization

Our optimization problem (equations (12)–(14)) depends on the Simulink model; we prefer to use a derivative-free evolutionary algorithm. PSO is a population-based iterative algorithm.^49,54 Populations are structured with some topology comprising bidirectional edges connecting particle pairs.

The power of PSO comes from the interaction of individuals and movement to a better objective value that occurs when the individuals interact. All particles in the same neighborhood communicate so they move toward their new position, influenced by the best objective value obtained by the individuals in the topology. Topologies of the neighborhood are scrutinized in Poli et al. ⁴⁷

In PSO, many particles are positioned in the problem’s search space, each evaluating the objective function at its present location. Each particle in the swarm is composed of three vectors in ${\mathbb{R}}^{27}$ ; the current position ${\mathit{x}}_iϵ{\mathbb{R}}^{27}$ , the previous best position ${\mathit{p}}_{\mathit{i}}ϵ{\mathbb{R}}^{27}$ and the velocity ${\mathit{v}}_{\mathit{i}}ϵ{\mathbb{R}}^{27}$ . At each iteration of the algorithm, ${\mathit{x}}_i$ vectors are taken as the controller gains. If ${\mathit{x}}_i$ is better, then the PID gains are stored in the second vector, ${\mathit{p}}_{\mathit{i}}$ . The best value in the neighborhood is denoted by $\mathit{g}ϵ{\mathbb{R}}^{27}$

Data generation for ANN

The controller structure given in Figure 3 was tested in many different scenarios. During these trials, if the vehicle is supposed to make a sharp turn between consecutive waypoints, the swing angle of the load may violate the constraints. This problem has been solved by adding a complementary controller, which computes the required controller inputs using an ANN. This structure (see Figure 4) with ANN creates complementary control inputs when necessary. It accepts velocity and position error vectors as inputs and generates necessary force and moments as its outputs. If the swing angle of the load is within the specified range, this structure will not contribute to the overall control inputs.

Figure 4.

Overall controller structure with ANN.

The input data for training the ANN has been created by running the structure shown in Figure 3 for different and challenging scenarios. Instances that did not meet the angle constraint were found in the results of these runs. Figure 5 shows the plot of a swing angle obtained from a scenario and the points that do not satisfy the constraint value. A particular optimal control problem is solved for each point (indicated by red dots in Figure 5) in the zoomed image. In these optimal control problems, the aim is to obtain the value of the complementary control inputs that will bring the swing angle below the constraint value. A complementary input $\mathrm{\Delta }{\mathit{v}}_{\mathrm{ij}}\left(t\right)(t=t_{\mathrm{ij}}-\mathrm{\Delta }T,\dots ,t_{\mathrm{ij}})$ have been added to the model. Then, the optimal control problem is formulated as follows:

$$\underset{\mathrm{\Delta }{\mathit{v}}_{\mathrm{ij}}\left(\mathit{t}\right)}{min}{\int }_{t_{\mathrm{ij}}-\mathrm{\Delta }T}^{t_{\mathrm{ij}}}{(max(\beta \left(t\right)-\epsilon ))}^2\mathrm{dt}$$

(15)

Figure 5.

Points at which swing angle constraint is violated.

subject to

$$\begin{array}{c}\mathrm{Dynamic}\mathrm{constraints}\mathrm{of}\mathrm{quadrotor}\mathrm{slung}\mathrm{load}\mathrm{system}\hfill \\ (i.e.,\mathrm{Eq}.3-8)\mathrm{and}\mathrm{PID}\mathrm{controllers}.\hfill \end{array}$$

where $i\in \{1,\dots ,N_s\}\mathrm{and}j\in \{1,\dots ,N_i\}$ are scenario and instance indices, respectively. The total number of scenarios is $N_s$ , and the total number of instances in scenario $i$ is $N_i$ . So, we have obtained $N_s{N}_i$ points to train ANN. Eighty percent of these points were taken as training and 20% as validation sets. In equation (15), $\mathrm{\epsilon }$ is inserted into the equation to satisfy this strict inequality constraint. The objective of the optimal control problem (equation (15)) is to determine the complementary inputs that should additionally be applied to the PID outputs such that the swing angle constraint is satisfied. The $\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}\left(t\right)$ values, which should be applied $\mathrm{\Delta }\mathrm{T}$ steps before the time step when the constraints are violated, were considered the optimal control problem variables. Since applying complementary inputs when the swing angle exceeds the constraint value does not guarantee the satisfaction of the constraint, the complementary controller should start contributing some time before a violation occurs. We preferred to use an interval of $\mathrm{\Delta }\mathrm{T}=200$ msec (i.e. 200 time-step) by trial and error to take care of the system’s inertia. These complementary input values obtained were used as ANN outputs during training. In the training process, ANN’s inputs and outputs are:

Input: ${\mathbf{\chi }}_{\mathit{ij}}\left(t\right)=\left[\varphi \theta p_x{p}_y{p}_z{V}_x{V}_y{V}_{z}pqr{\Omega }_x{\Omega }_y{\Omega }_z\right]$ ,

Output: $\mathrm{\Delta }{\mathit{v}}_{\mathit{ij}}\left(t\right)$ (complementary inputs for force and moments).

Training of the ANN

The scaled conjugate gradient backpropagation algorithm Møller ⁵⁵ is used because it avoids the line search per learning iteration. Instead, this algorithm uses gradient calculations to update weights instead of the Jacobian calculations used by other methods such as Levenberg-Marquardt and Bayesian regularization. The feedforward calculation has been given as follows;

$$N{N}_{\mathrm{out}}=h({\sum }_{j=1}^{n_H}{w}_{\mathrm{kj}}h(\sum _{i=1}^m{w}_{\mathrm{ji}}{\chi }_i+w_{j0})+w_{k0})$$

(16)

where the subscript $i$ indexes units on the input layer, $j$ for the hidden; $w_{\mathrm{ji}}$ denotes the input-to-hidden layer weights at the hidden unit $j,$ the subscript $k$ indexes units in the output layer (three, in our case) and $n_H$ denotes the number of hidden units (20, in our case). The activation function, $h$ , is the hyperbolic tangent sigmoid activation function;

$$h\left(s\right)=\frac{2}{1+e^{(-2s)}}-1$$

(17)

Simulation results

Results without using the complementary controller

A sufficiently challenging scenario with a 30° swing angle constraint is suggested to determine the PID coefficients using Algorithm 1. The top view of the resultant paths for the quadrotor and slung load is given in Figure 6. As shown in Figure 7, where the swing angle of the load is plotted with respect to time, the maximum swing angle has been obtained as 9.83°. In Figure 8, the graph of the vehicle’s speed is presented, and the speeds at the waypoints are marked. It is observed that the speeds at the waypoints are close to the desired speeds. There are minor deflections from desired speeds due to the swinging of the load. A more challenging scenario with very sharp turning points is used in this new study (Scenario-2). In this scenario, first, PID coefficients are found using an angle constraint of 15°.

Algorithm 1:

PSO-based PID tuning

Input: $w1$ , $w2$ ▹Adjustment Weights Input : $N$ , $S$ ▹Swarm Size, Neighbor fraction Input : $W$ ▹Inertia Range Output: $\overrightarrow{x}$  = $[K_1,...,K_{27}]$  ▹ PID Gains $\mathrm{Stop}\leftarrow \mathrm{false}$ ; GENERATESWARM $(N,S)$ ; for Each Particle-i do   CHECKBOUNDS $({\overrightarrow{x}}_i)$ ;    $f_i$  ← QSLMODEL $({\overrightarrow{x}}_i)$    $\mathrm{pbes}{t}_i$ ← $f_i$ ▹ initial personel bests end $\overrightarrow{g}$ ← FINDBESTNEIGHBOR $(\mathbf{f})$ ;▹ g: neighbor pos. while $\mathrm{Stop}$ = $\mathrm{false}$ do   for Each Particle-i do      ${\overrightarrow{v}}_i$ ←W ${\overrightarrow{v}}_i$ + $w_1$ ${\overrightarrow{u}}_1$ ⊗ $({\overrightarrow{p}}_i-{\overrightarrow{x}}_i)$ + $w_2$ ${\overrightarrow{u}}_2$ ⊗ $(\overrightarrow{g}-{\overrightarrow{x}}_i)$      ${\overrightarrow{x}}_i$ ← ${\overrightarrow{x}}_i$ + ${\overrightarrow{v}}_i$     CHECKBOUNDS $({\overrightarrow{x}}_i)$ ;      $f_i$ ← QSLMODEL $({\overrightarrow{x}}_i)$     if $f_i$ < $\mathrm{pbes}{t}_i$ then        $\mathrm{pbes}{t}_i$ ← $f_i$        $p_i$ ← $x_i$     end   end    $\overrightarrow{g}$ ← FINDBESTNEIGHBOR $(\mathbf{f})$ ;   UPDATENEIGHBORS();   ifObjective Limit is reached then      $\mathrm{Stop}\leftarrow \mathrm{true}$   end end

Figure 6.

Top view of the path of quadrotor slung-load system (Scenario-1).

Figure 7.

The swing angle of the load versus time (Scenario-1).

Figure 8.

Quadrotor’s speed versus time (Scenario-1).

The resultant paths followed by the quadrotor and slung load are presented in 3D in Figure 9. The top view of these paths is shown in Figure 10. The swing angle of the load is presented in Figure 11. As observed in this figure, the maximum swing angle is 12.49°. In Figure 12, the vehicle’s speed with respect to time is given.

Figure 9.

3D view of the path of quadrotor slung-load system (Scenario-2_1).

Figure 10.

Top view of the path of quadrotor slung-load system (Scenario-2_1).

Figure 11.

The swing angle of the load versus time (Scenario-2_1).

Figure 12.

Quadrotor’s speed versus time (Scenario-2_1).

Afterward, we wondered what happens to PID coefficients and the system’s performance if the swing angle constraint is relaxed to 30° (using the same scenario with a different swing angle constraint, i.e. Scenario-2_2).

When we run Algorithm 1 with this constraint on this scenario, we obtain the paths given in Figure 13, and the corresponding swing angle is given in Figure 14. The maximum swing angle is 25.81° in this case. We want to point out that by relaxing the system’s swing angle constraint (i.e. Figures 11 and 14) from 15° to 30°, task duration is decreased from 35.619 s to 25.764 s. This (i.e. decrease in task duration) is an expected result because as the system is more relaxed, it can behave faster to minimize costs further.

Figure 13.

Top view of the path of quadrotor slung-load system (Scenario-2_2).

Figure 14.

The swing angle of the load versus time (Scenario-2_2).

The PID coefficients found in the scenario with 15° swing angle constraint are tested on a new scenario (Scenario-3) shown in Figure 15. The resultant maximum swing angle (14.78°) can be observed in Figure 16, which is not much different than the result in Figure 11.

Figure 15.

Top view of the path of quadrotor slung-load system (Scenario-3).

Figure 16.

The swing angle of the load versus time (Scenario-3).

To check the sensitivity of PID autopilot coefficients concerning the path in the chosen scenario and the constraints, we performed the final test on a new scenario (Scenario-4), shown in Figure 17. This scenario has many sharp turns that may cause swing-angle constraint violations. Indeed, when the controller structure in Figure 3 is used, the swing angle of the load is given in Figure 18. The swing angle constraint for this scenario is 30°, but the resultant swing angle is about 38°. Therefore, since optimizing the PID coefficients for every possible path, desired speeds, and swing angle constraint is impossible, we have concluded that the autopilot structure should be extended/complemented in Figure 4.

Figure 17.

Resultant paths using controller structure without the complementary controller (Scenario-4).

Figure 18.

Swing Angles using controller structure without the complementary controller (Scenario-4).

Results using the complementary controller

The training data for ANN in improved structure is obtained from Figure 18 when the swing angle constraint was 30°. This new structure (with ANN) is tested on the same scenario (Scenario-4) with a swing angle constraint of 20°.

The results are shown in Figure 19; the path with the same set of waypoints is given. In this case, the control structure given in Figure 4 is employed. The ANN part is activated automatically when the swing angle of the load is within 18°–22°. Otherwise, it does not generate any delta inputs. As shown in Figure 20, the constraint is violated around 40th seconds by 0.0077°.

Figure 19.

Resultant paths using controller structure with complementary controller (Scenario-4).

Figure 20.

Swing Angles using controller structure with using the complementary controller (ANN has been trained for $[{\mathit{\beta }}_{\mathit{lower}}{\mathit{\beta }}_{\mathit{upper}}]=[1822]$ ° degrees swing angle constraint) (Scenario-4_1).

Figures 21 to 23 give delta inputs for thrust, roll, and pitch moments. As can be seen from these figures, the ANN structure is activated when the swing angle is within the defined interval (i.e. between 18° and 22°). Also, delta inputs are saturated, as in the case of PSO. We now change the swing angle constraint from 20° to 17° to better observe the Complementary Controller’s performance. Data is collected from the band $\left[{\beta }_{\mathrm{lower}}{\beta }_{\mathrm{upper}}\right]=\left[1519\right]$ ° to train the ANN. One can observe the swing angle graph in Figure 24. Note that the maximum value of the swing angle is 18.79°, which is slightly above the constraint for a short duration. You can compare this performance with the case when the complementary controller is trained for the band $\left[{\beta }_{\mathrm{lower}}{\beta }_{\mathrm{upper}}\right]=\left[2832\right]$ °; naturally, the resultant controller satisfies the constraint better as the constraint relaxes (Figure 25).

Figure 21.

Thrust generated by complementary controller versus time (Scenario-4_1).

Figure 22.

Roll Moment. It was generated by complementary controller versus time (Scenario-4_1).

Figure 23.

Pitch Moment. It was generated by complementary controller versus time (Scenario-4_1).

Figure 24.

Swing Angles using controller structure with using the complementary controller (ANN has been trained for $[{\mathit{\beta }}_{\mathit{lower}}{\mathit{\beta }}_{\mathit{upper}}]=[1519]$ ° degrees swing angle constraint) (Scenario-4_2).

Figure 25.

Swing Angles using controller structure with using the complementary controller (ANN has been trained for $[{\mathit{\beta }}_{\mathit{lower}}{\mathit{\beta }}_{\mathit{upper}}]=[2832]$ ° degrees swing angle constraint) (Scenario-4_3).

The overall simulation structure, including the state flow diagram (yellow block), is given in Figure 26. In this model, the state flow diagram manages waypoint guidance.

Figure 26.

Overall Simulation structure constructed in MATLAB/Simulink.

A few of these simulation studies were realized physically in an experimental setup. We use a similar method in Li et.al. ⁵⁶ to measure swing angles. Firstly, swing angle measurements of slung load were tested in a laboratory medium using the setup in Figure 27. The measurement accuracy was verified using a Single Board Computer (SBC) and a camera at the top of the setup. Next, the angular displacements of the slung load measurements were obtained using the same (Inertial Measurement Unit) IMU in the quadrotor-slung load system. Later, this swing angle measurement setup is integrated into a quadrotor.

Figure 27.

The swing angle measurement setup to collect IMU data.

Conclusions

This study designed a controller for a quadrotor carrying a suspended load. The quadrotor is supposed to reach its destination by passing through waypoints with speed constraints. This problem is solved successfully using the proposed complex control structure in Figure 4. In this structure, the primary controller adaptively changes the gains for speed and position as a function of their errors. The Complementary Controller is activated whenever the swing angle constraint is violated. The resultant structure performs much better than the one without the complementary structure. Simulations have been performed in MATLAB/Simulink environment. Measurements of the swing angle are first obtained in a laboratory setup, then tested on a quadrotor in the open field.

Studies on this problem will be extended to the slung load case with multiple quadrotors.