Particle Swarm Optimization and Cuckoo Search-Based Approaches for Quadrotor Control and Trajectory Tracking

Abstract

This paper explores the full control of a quadrotor Unmanned Aerial Vehicles (UAVs) by exploiting the nature-inspired algorithms of Particle Swarm Optimization (PSO), Cuckoo Search (CS), and the cooperative Particle Swarm Optimization-Cuckoo Search (PSO-CS). The proposed PSO-CS algorithm combines the ability of social thinking in PSO with the local search capability in CS, which helps to overcome the problem of low convergence speed of CS. First, the quadrotor dynamic modeling is defined using Newton-Euler formalism. Second, PID (Proportional, Integral, and Derivative) controllers are optimized by using the intelligent proposed approaches and the classical method of Reference Model (RM) for quadrotor full control. Finally, simulation results prove that PSO and PSO-CS are more efficient in tuning of optimal parameters for the quadrotor control. Indeed, the ability of PSO and PSO-CS to track the imposed trajectories is well seen from 3D path tracking simulations and even in presence of wind disturbances.

1. Introduction

Over the past few years, Unmanned Aerial Vehicles (UAVs) of type quadrotors or quadcopters have seen an increasing interest since their wide range of civilian and military applications. These quadrotors are considered the most complex flying machines due to many physical effects influencing their dynamics including aerodynamic effects, gravity, gyroscopic effects, friction, and inertia. However, they have advantages over conventional helicopters. Given that the left and the right motors rotate clockwise, and the front and the rear motors rotate counterclockwise, gyroscopic effects and aerodynamic torques tend to cancel in trimmed flight.

Many researchers modeled the quadrotor using the formalisms of Newton-Euler [1,2] and Euler-Lagrange [2,3,4]. For both methods, the quadrotor modeling is regarded as a delicate task and the model obtained using these approaches is strongly nonlinear, fully coupled, under-actuated (Six Degrees of Freedom System (6-DOF) with only four actuators) and dynamically unstable with complex behavior. In this work, we adopt the most used formalism of Newton-Euler to define the dynamics of the quadrotor.

Regarding the quadrotor’s control, there are several types of research done both in linear and nonlinear control methods [5,6,7]. For practicality, the linear controllers, especially PID (Proportional, Integral, and Derivative)/PD (Proportional and Derivative)/ PI (Proportional and Integral), are easy to design and simple to tune, which perform better in practical implementation [8]. In literature, PIDs are the most used for quadrotor control [2,9,10]. However, tuning the parameters of these controllers is very important to get best performances. In this scope, classical methods can be used. They include the methods of Ziegler-Nichols [11], Graham-Lathrop [12], Naslin [13], Cohen-Coon, Reference Model (RM) [14,15], and many other techniques.

Nature-inspired optimization algorithms can effectively resolve complex problems compared to classical and statistical methods. Some of these algorithms proposed in the literature are Genetic Algorithm (GA) [16,17], Particle Swarm Optimization (PSO) [18], Ant Colony Optimization (ACO) [19], Artificial Bee Colony (ABC) [20], and Cuckoo Search (CS) [21,22]. GA is inspired by the crossover, mutation, and recombination principles of genetics. PSO algorithm is inspired by social behavior patterns of organisms that live and interact within large groups. It incorporates swarming behaviors observed in flocks of birds and schools of fish. ACO is inspired by ants’ behavior in seeking food sources. The ABC algorithm is based on the natural behavior of honeybee swarm, and CS algorithm is motivated by the aggressive breeding of a bird called ‘cuckoo’. The PSO method has shown superior performances. It has features such as a straightforward algorithm, rapid convergence, easy implementation, and is computationally efficient when applied to a diverse set of optimization problems [15]. In CS, immigration and environmental specifications have the advantage to help cuckoos’ groups to converge and reach the best places for breeding and egg laying [21]. Another advantage of CS compared to PSO and GA is that it uses a smaller number of parameters to be tuned, which makes it more adaptable [21]. Due to these advantages, PSO and CS are still successfully applied in recent works for control problems [23,24,25], and PSO is especially used in tuning PID controllers for quadrotor’s control [26,27]. Therefore, regarding PSO and CS programs, there are many common points. Both use the same initialization matrices of particles in n-dimensional research space for PSO and nests with D eggs (problem’s dimension) for CS. The fitness function evaluates the quality of particles and nests in the same way. Then, the velocities of displacement of particles and nests are calculated in different ways for PSO and CS, which generates differences in the convergence speeds of both algorithms. For these reasons, it would be very interesting to analyze the social thinking and local search abilities of PSO and CS in quadrotor control.

The main contribution of our paper lies on using PSO and CS to improve results concerning the quadrotor full control. However, the CS algorithm suffers from a low convergence speed, since it uses a fixed step size over generations. To overcome this problem, Naik et al. [28] proposed an adaptive cuckoo search algorithm (ACS) with adaptive step size. The ACS converges to a near optimum solution faster than CS algorithm. However, the problem of ACS algorithm is that it does not achieve better solutions compared to CS in most cases. In our work, we propose the cooperative PSO-CS algorithm that combines the ability of social thinking in PSO with the local search capability of CS. Such a cooperative process is likely to offer proper guidance for cuckoos to the global best positions and to ensure a balance between exploitation and exploration of the search space. Consequently, we are interested in examining if it is possible for us to improve results of quadrotor control by effectively integrating these two algorithms with complementary strengths. Our work treats position and attitude control in order to cover all topics in quadrotor control. For this objective, the control design includes six PID controllers: Three for angles (ϕ, θ, ψ) and three for positions (x, y, z). These controllers are optimally tuned using the heuristics PSO and CS, the cooperative PSO-CS and the classical method of Reference Model (RM). A comparative study is done to highlight the efficiency of the proposed intelligent controllers of PSO and PSO-CS. Indeed, their robustness is demonstrated in presence of wind disturbance.

The rest of this paper is organized as follows. In Section 2, Newton-Euler formalism is used to establish the quadrotor dynamic model. In Section 3, the proposed algorithms of PSO, CS, the cooperative PSO-CS and RM are presented. In Section 4, the quadrotor full control is treated using PID controllers optimally tuned with the proposed methods. In the last section, our conclusions are given.

2. Quadrotor Dynamic Modeling

The quadrotor is a complex flying system, strongly nonlinear, multivariable, fully coupled, and under-actuated (6-DOF and only four control inputs). Therefore, its modeling is a delicate task. For better understanding of its dynamic modeling, different working hypotheses were assumed: A rigid and a symmetrical quadrotor structure (diagonal matrix of inertia); a rigid propellers (negligible effect of deformation during rotation); the lift and the drag forces are proportional to the square of the rotational speed of the rotors (a very close approximation of the aerodynamic behavior); and the center of mass is exactly the origin of reference related to the structure.

From the configuration shown in Figure 1, the quadrotor has two pairs of rotating rotors attached to the end of a cross and the control electronics is situated in the center of the cross. The front and the rear propellers spin counterclockwise, while the left and the right ones spin clockwise, which effectively neutralizes the unwanted reactive torque and allows the vehicle to fly without overturning.

Notice that the absolute position of the mass center is described by the three coordinates (x, y, z) and its attitude by the three Euler’s angles (ϕ, θ, ψ). Using Newton-Euler formalism [2], the dynamic model of the quadrotor can be expressed as:

$$\begin{array}{l}\ddot{\phi }=\frac{l{U}_2}{I_x}+\frac{\left(I_y-I_z\right)}{I_x}\left(\dot{\psi }cos\phi cos\theta -\dot{\theta }sin\phi \right)\left(\dot{\theta }cos\phi +\dot{\psi }sin\phi cos\theta \right)\\ -\frac{J_r\left({\omega }_1-{\omega }_2+{\omega }_3-{\omega }_4\right)}{I_x}\left(\dot{\psi }sin\phi cos\theta +\dot{\theta }cos\phi \right)\\ -\frac{K_{fax}}{I_x}\left({\dot{\phi }}^2-2\dot{\phi }\dot{\psi }sin{\theta }^2\right)\end{array}$$

(1)

$$\begin{array}{l}\ddot{\theta }=\frac{l{U}_3}{I_y}+\frac{\left(I_z-I_x\right)}{I_y}\left(\dot{\psi }cos\phi cos\theta -\dot{\theta }sin\phi \right)\left(\dot{\phi }-\dot{\psi }sin\theta \right)\\ -\frac{J_r\left({\omega }_1-{\omega }_2+{\omega }_3-{\omega }_4\right)}{I_y}\left(\dot{\psi }sin\theta -\dot{\varphi }\right)\\ -\frac{K_{fay}}{I_y}\left({\dot{\theta }}^{2}cos{\phi }^2+2\dot{\phi }\dot{\psi }sin\phi cos\phi cos\theta +{\dot{\psi }}^{2}sin{\phi }^{2}cos{\theta }^2\right)\end{array}$$

(2)

$$\begin{array}{l}\ddot{\psi }=\frac{U_4}{I_z}+\frac{\left(I_x-I_y\right)}{I_z}\left(\dot{\psi }sin\phi cos\theta +\dot{\theta }cos\phi \right)\left(\dot{\phi }-\dot{\psi }sin\theta \right)\\ -\frac{J_r\left({\omega }_1-{\omega }_2+{\omega }_3-{\omega }_4\right)}{I_y}\left(\dot{\psi }sin\theta -\dot{\phi }\right)-({\dot{\theta }}^{2}sin{\phi }^2\\ -2\dot{\phi }\dot{\psi }sin\phi cos\phi cos\theta )+\frac{K_{faz}}{I_z}{\dot{\psi }}^{2}cos{\phi }^{2}cos{\theta }^2\end{array}$$

(3)

$$\ddot{x}=\frac{\cos \phi \cos \phi \sin \theta +\sin \phi \sin \psi }{m}{U}_1-\frac{K_{ftx}}{m}\dot{x}$$

(4)

$$\ddot{y}=\frac{\sin \psi \sin \theta \cos \phi -\cos \psi \sin \phi }{m}{U}_1-\frac{K_{fty}}{m}\dot{y}$$

(5)

$$\ddot{z}=-g+\frac{\cos \phi \cos \theta }{m}{U}_1-\frac{K_{ftz}}{m}\dot{z}$$

(6)

with

$$\{\begin{array}{l}{U}_1=b\left({\omega }_1{}^2+{\omega }_2{}^2+{\omega }_3{}^2+{\omega }_4{}^2\right)\\ U_2=b\left({\omega }_4{}^2-{\omega }_2{}^2\right)\\ U_3=b\left({\omega }_3{}^2-{\omega }_1{}^2\right)\\ U_4=d\left({\omega }_1{}^2-{\omega }_2{}^2+{\omega }_3{}^2-{\omega }_4{}^2\right)\end{array}$$

(7)

where m is the total mass of the quadrotor; l the distance from the center of the body to the rotor; I_x, I_y and I_z are the inertias around x, y, and z axis; b and d thrust and drag coefficients respectively; K_fax, K_fay, and K_faz are frictions aerodynamics coefficients; K_ftx, K_fty, and K_ftz are translation drags coefficients; J_r is the rotor inertia; g represents gravity constant; ω_i is the angular speed of the motor i and finally ϕ, θ, and ψ are rotation around the roll, pitch, and yaw axis, respectively.

3. Proposed Methods

Intelligent and classical control methods of Particle Swarm Optimization (PSO), Cuckoo Search (CS), the cooperative PSO-CS, and the classical Reference Model (RM) are detailed in this section for further reason to be applied in quadrotor control.

3.1. PSO Algorithm

Particle Swarm Optimization (PSO) is a class of stochastic algorithms, which was developed by Eberhart and Kennedy in 1995 for hard optimization problems. The basic principle of PSO is inspired by the social behavior of animals moving in swarm as bird flocking. To search for food, each bird called particle, flies in the space of solutions, and determines its speed according to its personal experience and the information gained through interaction with other swarm members [18].

PSO is employed to intelligently select optimal parameters from N particles. The initialization matrix contains N particles dispersed in a D-dimensional search space. Like most optimization techniques, PSO requires a fitness function relevant to the particle’s position. Each particle i stores its best position Pb_i(t+1) and the best solution in its vicinity Pg(t+1), which is the position of the particle that has the smallest fitness value in the swarm as expressed in Equation (9). The mechanism of displacement of each particle is managed by three rules. First, the particle tends to follow the direction of its current velocity. Second, it wants to move toward its best position. Finally, it tends to move to the best position reached by its neighbors [18,29]. In fact, the new velocity matrix V_ij and position matrix X_ij of the particle i = {1, 2, …, N} in the search space of dimension D, with j = {1, 2, …, D} are calculated at iteration (t+1), according to Equations (10) and (11).

$$P{b}_i\left(t+1\right)=\{\begin{array}{ll}{X}_i\left(t+1\right)& iff\left(X_i\left(t+1\right)\right)<f\left(P{b}_i\left(t\right)\right)\\ P{b}_i\left(t\right)& else\end{array}$$

(8)

$$Pg=mi{n}_{i=1,2\dots N}f\left(P{b}_i\left(t+1\right)\right)$$

(9)

$$V_{ij}\left(t+1\right)=w{V}_{ij}\left(t\right)+R_1{C}_1\otimes \left(P{b}_{ij}\left(t\right)-X_{ij}\left(t\right)\right)+R_2{C}_2\otimes \left(P{g}_{ij}\left(t\right)-X_{ij}\left(t\right)\right)$$

(10)

$$X_{ij}\left(t+1\right)=X_{ij}\left(t\right)+V_{ij}\left(t+1\right)$$

(11)

where

• Pb_ij is the best position found by the particle i;
• Pg_ij is the best position found by the neighborhood;
• w, C₁, and C₂ are weighting coefficients;
• R₁ and R₂ are random variables generated from a uniform distribution in [0,1];
• ⊗ means element wise multiplications.

The weighting coefficients govern the process of finding the best solution and reflect the sociability of the particles. In order to move the particle towards the best solutions, selection of these coefficients must ensure a compromise between local and global exploration of the search space. In previous works [30,31], the constants C₁ and C₂ must verify the condition C₁ + C₂ ≤ 4. Also, experimental results in Reference [30] show that [1/2((C₁ + C₂) − 1)] < w ≤ 1. PSO search works in G iterations and the coordinates of the obtained global best position Pg are the parameters optimally tuned.

3.2. CS Algorithm

Cuckoo Search (CS) algorithm, proposed by Yang and Deb in 2009, is based on the life of the ”cuckoo” bird. The basic principle of this optimization algorithm is the specific breeding and egg laying of this bird. In the habitat of other host birds, adult cuckoos lay some eggs that grow and become mature cuckoos if they are not discovered and removed by host birds. The immigration of groups of cuckoos and environmental specifications hopefully lead them to converge and reach the best places for reproduction and breeding [21].

The primary population of CS is composed of N nests representing a set of solutions, where each nest has multiple eggs (D eggs). As in nature, each cuckoo dedicates from five to 20 eggs. These values are used as the limits of the dimension D. Each cuckoo lays one egg at a time and dumps it in a randomly chosen nest. The best nests with high-quality eggs (solutions) will carry over to the next generations, where the profit of a solution is obtained by evaluation of a fitness function F of the habitat, so the profit is an array of 1 × D. The number of available host nests is fixed, and a host can discover an alien egg with probability P_a from [0, 1]. In this case, the host bird can either throw the egg away or abandon the nest to build a completely new nest in a new location. For simplicity, this assumption can be approximated by a fraction P_a of the N nests being replaced by new nests, having new random solutions [21].

To replace solutions in the nests with new solutions, Lévy flights mechanism is used. A new solution X_i(t + 1) for cuckoo I is given by Equation (12), where ⊗ is entry wise product, similar to that used in PSO, α > 0 is the step size, and R is a random variable generated from a uniform distribution in the interval [0, 1] in order to provide a stochastic weight.

$$X_i\left(t+1\right)=X\left(i\right)+\alpha R\otimes Lévy\left(\lambda \right)$$

(12)

The random walk via Lévy flight is more efficient in exploring the search space, as its step length is drawn from a Lévy distribution and it is much longer in the long run. The Mantegna algorithm is used as Lévy stable distribution in most problems to decide the step length. Thus, the step length S from Mantegna algorithm can be written as represented by Equation (13), where u and v are samples drawn from Gaussian normal distribution given by Equation (14), and σ_u² is the variance of the distributions given by Equation (15), where Γ is the gamma function [22].

$$Lévy\left(\lambda \right)\approx S=u{\left|v\right|}^{-\frac{1}{\lambda }}\left(1<\lambda \le 3\right)$$

(13)

$$\{\begin{array}{c}u\approx Norm\left(0,{\sigma }_u^2\right)\\ v\approx Norm\left(0,1\right)\end{array}$$

(14)

$${\sigma }_u^2={\left\{\frac{\Gamma \left(1+\lambda \right)\sin \left(\pi \lambda /2\right)}{\Gamma \left[\left(1+\lambda \right)/2\right]\lambda 2^{\left(\lambda -1\right)/2}}\right\}}^{1/\lambda }$$

(15)

The fraction P_a of worse solutions are generated as given by Equation (16), where X_j(t) and X_k(t) are two random solutions chosen by random permutation, H is a Heaviside function, and r is a random number drawn from a normal distribution.

$$X_i\left(t+1\right)=X_i\left(t\right)+\left(P_a-r\right)\left(X_j\left(t\right)-X_k\left(t\right)\right)$$

(16)

3.3. Cooperative PSO-CS Algorithm

The initialization matrix in PSO contains N particles dispersed in a D-dimensional search space, when in CS, this matrix is composed of N nests and each nest has D eggs (problem’s dimension). In both PSO and CS, this primary population can be regarded of dimension D × N and the quality of these solutions is evaluated in the same way. The global best particle in PSO (or the best nest in CS) is the particle (or the nest) that has the smallest fitness value among all potential solutions. Then, the velocities of displacement of particles and nests are calculated in different ways in PSO and CS. However, due to the fast convergence of PSO algorithm and the low convergence speed of the CS algorithm, the cooperation between PSO and CS would be very interesting to combine the ability of social thinking in PSO with the local search capability of CS. This cooperative PSO-CS aims to modify the speed equation of displacement by combining the Lévy flights random walks of cuckoos and the movement of particles toward the global best solution in their vicinity Pg. New solutions X_i (t + 1) are given by:

$$X_i\left(t+1\right)=X_i\left(t\right)+\alpha R\otimes Lévy\left(\lambda \right)+C_2{R}_2\otimes \left(Pg\right(t)-X_i(t\left)\right)$$

(17)

The use of cooperative PSO-CS helps to improve searching at global and local scales, which allows a balance between exploration and exploitation of the search space. Both algorithms are combined and involved to produce results. This combination of PSO and CS capabilities contribute to increase particle diversification in the entire search space and ensure good coverage during iterations. Thanks to this cooperation, CS algorithm can be proposed to control the quadrotor. Indeed, to evaluate the efficiency of the proposed PSO-CS method, comparisons with the intelligent approaches of PSO and CS are presented in the next sections for quadrotor’s control and trajectory tracking.

3.4. RM Method

Reference Model (RM) technique is a linear method that approaches the behavior of any system of “n” order to the desired behavior of a reference system of first or second order. If the reference behavior is that of a first order, the desired settling time T_s is fixed and the value of τ is deduced, such that T_s = 3τ. The dominant pole is then to place at −1/τ and the (n−1) other poles on left of −1/τ in the complex plane [14,15]. If the reference behavior is that of a second order, with imposed values for damping ζ and natural pulsation ω_n, the two dominant poles are to place at −ζω_n ± jω_n√(1 − ζ²) and the (n−2) other poles on left of −ζω_n in the complex plane [14,15].

The closed loop transfer function H(p), of the system including the controller and its denominator D_H(p) are computed. By taking into account some desired specifications (T_s, τ or ζ, ω_n), the characteristic polynomial R_s(p) of the reference system is evaluated.

$${\mathrm{R}}_{\mathrm{s}}(\mathrm{P})=\{\begin{array}{ll}\left(p+\frac{1}{\tau }\right){\left(p+\frac{\alpha }{\tau }\right)}^{n-1}& ifafirstorderbehaviordesired\\ \left(p^2+2\zeta {\omega }_{n}p+{\omega }_n{}^2\right){\left(p+\alpha \zeta {\omega }_n\right)}^{n-2}& else\end{array}$$

(18)

The parameter “α” in (18) is made much greater than 1 to ensure the dominance of the poles −1/τ or −ζω_n ± jω_n√(1 − ζ²). By equating the coefficients of D_H(p) and D_s(p), the controller parameters can be extracted [14,15].

4. Quadrotor Control

4.1. Virtual Control

To manage the system dynamics, a simplified model is needed. By considering small rotations (cos x = 1 and sin x = x), the equations become linear and the PID controllers’ inputs can be clearly seen. Consequently, the control inputs for ϕ, θ, and ψ are U₂, U₃, and U₄, respectively. To deal with the under-actuated part, ϕ, θ, and U₁ are used to create three virtual control inputs so that every output will be controlled separately. Physically, these virtual controls mean that the translation motion along x, y, and z are controlled indirectly by the three common inputs (ϕ, θ, and U₁). These virtual controls are given by Equation (19) and the desired trajectories of roll (ϕ_d), pitch (θ_d) and U_1d are obtained by inverting Equation (19).

$$\{\begin{array}{l}{U}_x=(cos\phi cos\psi sin\theta +sin\phi sin\psi )U_1\\ U_y=(sin\psi sin\theta cos\phi -cos\psi sin\phi )U_1\\ U_z=\left(cos\phi cos\theta \right)U_1\end{array}$$

(19)

$$\{\begin{array}{l}{\phi }_d={\sin }^{-1}\left(\frac{U_{x}sin{\psi }_d-U_{y}cos{\psi }_d}{\sqrt{U_x^2+U_y^2+U_z^2}}\right)\\ {\theta }_d={\tan }^{-1}\left(\frac{U_{x}cos{\psi }_d+U_{y}sin{\psi }_d}{cos{\phi }_d}\right)\\ U_{1d}=\sqrt{U_x^2+U_y^2+U_z^2}\end{array}$$

(20)

4.2. Control Law Design

For an appropriate control of the UAV quadrotor system described by Equations (1)–(6), PID controllers are used. These controllers are simple to design and robust. The complete control scheme for the quadrotor is represented in Figure 2. It can be clearly seen from this figure that three PID controllers are used for positions (x, y, z) and three PIDs for attitudes (ϕ, θ, ψ). The inputs are divided to two parts, the desired and the sensor signals. Desired signals for positions and yaw angle (x_d, y_d, z_d, ψ_d) are fixed by the pilot or the autopilot program. Sensor signals for desired roll and pitch angles (ϕ_d, θ_d, and U_1d) are computed from the three virtual inputs U_x, U_y, and U_z.

The transfer functions C_s(p) of PID controllers are given in Equation (21), where s = {ϕ, θ, ψ, x, y, z}. Their gains (K_ps, K_is, K_ds) should be optimally tuned for the quadrotor in order to ensure stable, fast, precise, and robust responses. For this objective, PSO, CS, PSO-CS, and RM methods are applied to realize an optimal setting of these parameters.

$$C_s\left(p\right)=K_{ps}+\frac{K_{is}}{p}+K_{ds}p$$

(21)

4.3. Quadrotor Intelligent Control Using PSO, CS, and PSO-CS

To control the nonlinear and fully coupled quadrotor defined by Equations (1)–(6), it is suitable to use a method directly applicable to nonlinear systems as PSO, CS, and the cooperative PSO-CS. Instead of linear PID controllers, the intelligent and adaptive controllers can look for optimal PIDs parameters for the six quadrotors’ outputs [ϕ, θ, ψ, x, y, z].

The initial populations (N = 200) of PSO, CS, and PSO-CS are randomly dispersed in a search space that is chosen sufficient to contain all possible solutions [0, 200] and its dimension D is set to 18. To define the exploration capacity of each particle in order to improve the convergence of the process and to not exceed the velocity limits, we tested numerous values of w, and the compromise between local and global exploration in PSO was achieved for w = 0.8. The confidence coefficients C₁ and C₂ define the exploration capacity. They determine the degree of influence of the past positions of the particle itself and those of the other particles of the swarm in order to regulate the velocities relative to the best local and global positions, thus making it possible to determine the movements relative to these best positions. C₁R₁ and C₂R₂ are represented by matrices with dimension of 200 × swarm size. We considered that C₁R₁ takes random values in [0, 0.8] to avoid the problem of fast convergence, when C₂R₂ takes random values in [0, 1.2] to give more importance to the global best solution Pg (the interval of C₂R₂ is large to that of C₁R₁). In CS, the parameters used in experiments are as follows: Abandon probability P_a = 0.25, and the Lévy flights settings, α = 0.1 and λ = 1.5. In PSO-CS, the same settings of PSO and CS are conserved: P_a = 0.25, and the Lévy flights settings, α = 0.1, λ = 1.5 and C₂ with random values in [0, 1.2].

To evaluate the optimization performance of the quadrotor’s responses using the proposed approaches, the fitness function (profit) is defined in a similar way in order to minimize the differences between the desired and the controlled outputs responses of the quadrotor. This fitness function is defined in Equation (22) as the sum of the errors E_k that characterize the difference in behavior between the inputs and the outputs of the quadrotor based on the Integral Square Error (ISE), with l = 6 and k = {1, 2, 3, 4, 5, 6}.

$$F={{\displaystyle \sum }}_{k=1}^l{{\displaystyle \int }}_0^{\infty }{E}_k^2\left(t\right)dt$$

(22)

$$ISE={{\displaystyle \int }}_0^{\infty }{E}_k^2\left(t\right)dt$$

(23)

The performances found with ISE are compared to those established using three other performance indices: Integral Absolute Error (IAE), Integral Time Absolute Error (ITAE), and Integral Time Square Error (ITSE).

$$IAE={{\displaystyle \int }}_0^{\infty }\left|E_k\left(t\right)\right|dt$$

(24)

$$ITAE={{\displaystyle \int }}_0^{\infty }t\left|E_k\left(t\right)\right|dt$$

(25)

$$ITSE={{\displaystyle \int }}_0^{\infty }t{E}_k^2\left(t\right)dt$$

(26)

The maximum number of generations for the three programs (PSO, CS and PSO-CS) is fixed as the stop criterion and set to 20.

4.4. Quadrotor Classical Control Using RM

To use the linear RM technique for quadrotor control, the complete nonlinear Equations (1)–(6) should be simplified. The most used technique to simplify these nonlinear equations is to consider weak angles variations (reduced gyroscopic effects) and retain the linear portion, the Equations (1)–(6) become:

$$\{\begin{array}{l}{I}_x\ddot{\phi }=l{U}_2\\ I_y\ddot{\theta }=l{U}_3\\ I_z\ddot{\psi }=U_4\\ \ddot{x}=\frac{\cos \phi \cos \psi \sin \theta +\sin \phi \sin \psi }{m}{U}_1-\frac{K_{ftx}}{m}\stackrel{\cdot }{x}\\ \ddot{y}=\frac{\sin \psi \sin \theta \cos \phi -\cos \psi \sin \phi }{m}{U}_1-\frac{K_{fty}}{m}\stackrel{\cdot }{y}\\ \ddot{z}=-g+\frac{\cos \phi \cos \theta }{m}{U}_1-\frac{K_{ftz}}{m}\stackrel{\cdot }{z}\end{array}$$

(27)

The inputs U₂, U₃, U₄, U_x, U_y, and U_z affect respectively angles of roll, pitch, and yaw, and positions x, y, and z, the reason for their usage as the PID controllers’ inputs. By applying the Laplace transform, we obtain the following transfer functions for angles (ϕ, θ, ψ) and positions (x, y, z).

$$\{\begin{array}{l}{G}_{\phi }\left(p\right)=\frac{\phi \left(p\right)}{U_2\left(p\right)}=\frac{l}{I_x{p}^2}\\ G_{\theta }\left(p\right)=\frac{\theta \left(p\right)}{U_3\left(p\right)}=\frac{l}{I_y{p}^2}\\ G_{\psi }\left(p\right)=\frac{\psi \left(p\right)}{U_4\left(p\right)}=\frac{1}{I_z{p}^2}\\ G_x\left(p\right)=\frac{x\left(p\right)}{U_x\left(p\right)}=\frac{K_x}{p\left({\tau }_{x}p+1\right)}\\ G_y\left(p\right)=\frac{y\left(p\right)}{U_y\left(p\right)}=\frac{K_y}{p\left({\tau }_{y}p+1\right)}\\ G_z\left(p\right)=\frac{z\left(p\right)}{U_z\left(p\right)}=\frac{K_z}{p\left({\tau }_{z}p+1\right)}\end{array}$$

(28)

where K_j = 1/K_ftj, τ_j = m/K_ftj and j = {x, y, z}.

For quadrotor’s attitudes and positions control, the behavior of the linear equations defined in (28) is approached to a desired reference one of a first order system. The desired behaviors for angles and positions correspond to the characteristic polynomials R_s(p), where s = {ϕ, θ, ψ, x, y, z} and α, α₁, α₂ are taken equal to 30, 2 and 30, respectively.

$${\mathrm{R}}_s\left(p\right)=\left(p+\frac{1}{\tau }\right)\left(p+\frac{{\alpha }_1}{\tau }\right)\left(p+\frac{{\alpha }_2}{\tau }\right)$$

(29)

For all cases, we calculate the closed loop transfer functions Hs(p) of the six systems including the PID controllers, for which we desire search their parameters. The computed denominators D_Hs(p) are given in Equation (30).

$$\{\begin{array}{l}{D}_{H\phi }\left(p\right)=p^3+\frac{l{K}_{d\phi }}{I_x}{p}^2+\frac{l{K}_{p\phi }}{I_x}p+\frac{l{K}_{i\phi }}{I_x}\\ D_{H\theta }\left(p\right)=p^3+\frac{l{K}_{d\theta }}{I_x}{p}^2+\frac{l{K}_{p\theta }}{I_x}p+\frac{l{K}_{i\theta }}{I_x}\\ D_{H\psi }\left(p\right)=p^3+\frac{l{K}_{d\psi }}{I_x}{p}^2+\frac{l{K}_{p\psi }}{I_x}p+\frac{l{K}_{i\psi }}{I_x}\\ D_{Hx}\left(p\right)=p^3+\frac{K_{dx}+K_{ftx}}{m}{p}^2+\frac{K_{px}}{m}p+\frac{K_{ix}}{m}\\ D_{Hy}\left(p\right)=p^3+\frac{K_{dy}+K_{fty}}{m}{p}^2+\frac{K_{py}}{m}p+\frac{K_{iy}}{m}\\ D_{Hz}\left(p\right)=p^3+\frac{K_{dz}+K_{ftz}}{m}{p}^2+\frac{K_{pz}}{m}p+\frac{K_{iz}}{m}\end{array}$$

(30)

5. Results and Discussions

Simulations of the proposed intelligent and classical PID controllers are carried out on the quadrotor’s attitudes (ϕ, θ, ψ) and positions (x, y, z). These PIDs are tested firstly in the case of hovering to a determined position coordinates in order study the performances of the proposed intelligent PSO, CS and PSO-CS and the classical RM controllers. Then, for two different desired trajectories, the controllers’ gains are optimized through the intelligent proposed approaches. The main objective is to confirm the ability of these intelligent nature inspired methods to control the quadrotor and to efficiently realize different imposed trajectories even under the presence of disturbances.

Table 1. Quadrotor’s parameters.

Parameter (unit)	Value
m (Kg)	0.65
g (m/s2)	9.806
l (m)	0.4
b (N/rad/s)	2.9842 × 10−5
d (N.m/rad/s)	7.5 × 10−7
Ix, Iy (Kg.m2)	7.5 × 10−3
Iz (Kg.m2)	1.3 × 10−3
Jr (Kg.m2)	2.8385 × 10−5
Kfax, Kfay (N/rad/s)	5.567 × 10−4
Kfaz (N/rad/s)	6.354 × 10−4
Kftx, Kfty(N/rad/s)	5.567 × 10−4
Kftz (N/rad/s)	3.354 ×10−4

gives the parameters of the quadrotor and

Table 2. Trajectories’ conditions and dynamics.

	Initial Conditions	Trajectory Description	Attitude’s Saturation
Rectangular Trajectory	[ϕ, θ, ψ, x, y, z] = [0, 0, 0, 0, 0, 10]	$\begin{array}{l}{x}_d\left(t\right)=2\left[u\left(t-4\right)-u\left(t-12\right)\right]\\ y_d\left(t\right)=2\left[u\left(t-8\right)-u\left(t-16\right)\right]\\ z_d\left(t\right)=10u\left(t-20\right)\\ {\psi }_d\left(t\right)=0\end{array}$	−90° < ϕ, θ < 90
Helical Trajectory	[ϕ, θ, ψ, x, y, z] = [0, 0, 0, 0, 0, 2]	$\begin{array}{l}{x}_d\left(t\right)=2\sin \left(2\pi t\right)\\ y_d\left(t\right)=2\sin \left(2\pi t-\pi /2\right)\\ z_d\left(t\right)=tu\left(t\right)-\left(t-18\right)u\left(t-18\right)\\ {\psi }_d\left(t\right)=0\end{array}$

presents the suggested trajectories’ dynamics and their conditions.

5.1. Comparisons of Intelligent and Classical Control Performances

To get the quadrotor to hover at a determined position coordinates, desired signals for positions [x_d, y_d, z_d] are considered step signals of 1 m and desired yaw angle (ψ_d) is fixed to 1 degree.

Table 3. Optimal Proportional, Integral and Derivative (PID) parameters obtained with Particle Swarm Optimization (PSO), Cuckoo Search (CS), Particle Swarm Optimization-Cuckoo Search (PSO-CS), and Reference Model (RM) for the quadrotor.

Quadrotor’s Output	Optimal PID Gains
	PSO			CS
	Kp	Ki	Kd	Kp	Ki	Kd
x	42.829	0.167	12.366	27.345	0.865	14.023
y	62.881	0.226	16.065	36.714	0.505	14.056
z	34.545	14.508	9.498	14.756	8.647	8.077
ϕ	14.740	1.749	0.975	17.948	0.178	1.055
θ	13.872	0.809	0.555	19.430	0.867	1.806
ψ	15.299	0.702	0.994	7.967	0.198	0.278
	PSO-CS			RM
	Kp	Ki	Kd	Kp	Ki	Kd
x	69.822	−0.011	16.974	59.150	29.250	30.550
y	45.880	0.008	15.067	39.650	19.500	20.800
z	31.859	17.724	6.715	26.066	10.394	16.865
ϕ	206.617	3.941	1.493	2.739	2.608	0.783
θ	121.842	−1.648	2.217	3.293	1.630	1.695
ψ	38.469	0.469	1.466	1.092	1.040	0.31

summarizes the PID controllers gains computed online using PSO, CS, and PSO-CS, and offline with the conventional RM for the quadrotor’s outputs (ϕ, θ, ψ, x, y, z). Figure 3, Figure 4, Figure 5 and Figure 6 show the step responses for the corresponding outputs [x, y, z, ψ] obtained with these methods. Figure 7 and Figure 8 show the various obtained curves for the desired sensor signals as well as the controlled responses.

Table 4. Performances obtained with Particle Swarm Optimization (PSO), Cuckoo Search (CS), Particle Swarm Optimization-Cuckoo Search (PSO-CS), and Reference Model (RM) for the quadrotor.

Output	Performances	Method
		PSO	CS	PSO-CS	RM
x	Ts (s)	0.730	1.384	0.522	3.779
	Mp (%)	0.412	1.769	0.446	22.869
	ISE	0.227	0.346	0.1650	0.365
	IAE	0.324	0.492	0.234	0.714
	ITSE	0.034	0.083	0.018	0.176
	ITAE	0.088	0.262	0.038	0.847
y	Ts (s)	0.675	1.018	0.782	3.517
	Mp (%)	0.497	0.898	0.010	30.741
	ISE	0.237	0.296	0.226	0.441
	IAE	0.322	0.410	0.319	0.802
	ITSE	0.036	0.058	0.034	0.245
	ITAE	0.083	0.152	0.076	0.954
z	Ts (s)	1.070	1.821	0.813	1.855
	Mp (%)	0.107	2.463	0.002	3.467
	ISE	0.243	0.406	0.208	0.342
	IAE	0.397	0.587	0.308	0.617
	ITSE	0.050	0.129	0.031	0.122
	ITAE	0.134	0.312	0.081	0.521
ϕ	ISE	0.028	0.006	0.005	0.093
	IAE	0.197	0.053	0.032	0.431
	ITSE	0.016	0.001	0.000	0.076
	ITAE	0.388	0.036	0.013	0.741
θ	ISE	0.334	0.129	0.024	0.912
	IAE	0.886	0.271	0.068	1.509
	ITSE	0.386	0.026	0.001	1.238
	ITAE	2.006	0.134	0.023	2.822
ψ	Ts (s)	0.172	0.207	0.097	2.228
	Mp (%)	0.300	20.424	0.047	16.129
	ISE	0.036	0.037	0.022	0.138
	IAE	0.073	0.068	0.038	0.331
	ITSE	0.0012	0.012	0.000	0.033
	ITAE	0.038	0.011	0.006	0.244

gives the obtained quadrotor’s responses performances: Settling time (T_s), maximum overshoot (Mp), and the errors previously mentioned (ISE, IAE, ITAE, and ITSE).

For control systems, the RM method often respects the fixed set of specifications and gives acceptable results [15]. However, this method is linear, and it cannot be applied directly to quadrotor control (the nonlinear quadrotor model). It can be observed from the obtained responses in Figure 3 and Figure 4 and the performances in Table 4 that the search for the parameters of PIDs by RM is not optimal: Responses with overshoots, settling times, and some error values exceed the amplitude of the desired signal. These unsatisfactory performances of the linear RM method in search of PID gains can be explained by the fact that the linear transfer functions defining the system (26) already have integrations.

CS results present a considerable overshoot in yaw angle and the errors values of ISE and IAE are around 50% of the desired altitude z, which proves that ISE value weight more significantly the quality of solutions and more adequate to evaluate the fitness function. However, the PSO approach is suitable and results in better performances. Comparisons of PSO results to those of Genetic Algorithms presented in Reference [16] show that PSO reduces the settling time T_s from 2.5 s to 0.172 s for yaw angle. For altitude z, the value of T_s is increased to 1.07 s instead of 0.4 s for GA. These results illustrate that PSO and GA techniques are highly competitive to explore the optimal solution in quadrotor control. In the same project of quadrotor control, the work presented in Reference [32] suggests three intelligent controllers, multilayer neural networks controller based on reference model (RMNN), Neural Networks with PD’s controller behavior (CPDNN), and Fuzzy Logic to adjust PID controllers’ gains (Fuzzy-PID). The proposed Fuzzy-PID and CPDNN achieves better performances compared to RMNN. Thus, comparing PSO results to those of Fuzzy-PID and CPDNN confirm that PSO is highly competitive to soft computing methods.

Considering the Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 and Table 4, the controllers tuned using PSO-CS have good performances for positions and attitudes. Thus, PSO-CS results are seen concurrent to those obtained using PSO, where responses lead to fast settling times with no overshoot and reduced errors. The desired sensor signal curves in Figure 7 and Figure 8 (desired roll and pitch angles (ϕ_d and θ_d)) depend on the outputs of the controllers PID-x, PID-y, and PID-z. It is observed from this figure that these intelligent approaches, and especially PSO and PSO-CS, allow to follow perfectly the sensor signals for roll and pitch (ϕ_d and θ_d), and give controlled responses close to the desired curves.

These programs are implemented using Matlab R2016a and performed on a computer with a processor Intel^® Core™ i7-3770 CPU of 3.40 GHz and 8.0 GB of RAM. From the conducted experiments, we can get conclusions that for the same goal, the PSO algorithm spends less CPU time (10 min 45 s) than the PSO-CS algorithm (14 min 12 s) and the CS algorithm (12 min 15 s). During the calculation of the velocity update equation, PSO-CS consumes an additional CPU time as the algorithm combines the two matrices generating Lévy flights random walks of cuckoos and that generating the movement of particles toward the global best solution Pg. However, the PSO-CS and PSO algorithms get higher performances (reduced errors) than the CS algorithm and RM method. These results are achieved thanks to the good adjustment of PSO, CS and PSO-CS parameters: Position and velocity initializations, the appropriate choice of the fitness function, the good tuning of the weighting coefficients w, C₁, and C₂, the chosen probability P_a, and the Lévy flights settings for α and λ.

5.2. Rectangular Trajectory Tracking

This experimental simulation on the quadrotor engages to test the PIDs’ efficiency when the vehicle tracks a square trajectory without overturning. In this mode, inputs for desired positions (x_d, y_d, z_d) are considered pulse signals and desired yaw attitude is set to 0°. Figure 9 shows the tracking of each position coordinate (x, y, z) and the tracking trajectories in 3D plan obtained using the proposed methods. Responses for desired and controlled roll (ϕ) and pitch (θ) attitudes are illustrated in Figure 6.

Table 5. Errors (ISE-Integral Square Error) obtained with Particle Swarm Optimization (PSO), Cuckoo Search (CS), Particle Swarm Optimization-Cuckoo Search (PSO-CS), and Reference Model (RM) for the quadrotor (rectangular trajectory).

Method	Errors ISE for Output
	x	y	z	ϕ	θ	ψ
PSO	1.083	0.677	3.287	0.119	0.093	7.586 × 10−36
CS	1.447	1.327	21.46	0.152	0.318	1.001 × 10−35
PSO-CS	0.955	0.795	1.453	0.107	0.4022	1.651 × 10−36
RM	2.573	1.561	33.756	0.277	0.712	1.210 × 10−34

resumes the errors ISE used to compute fitness function.

It can be noticeably observed from Figure 9 and Table 5 that the adaptive autotuned controllers using PSO, CS, and PSO-CS have good performances (reduced errors) for positions (x, y, z) control. PSO-CS and PSO responses always lead to fast settling times with no overshoot for x position. RM presents overshoots for positions and CS shows few overshoots for y position and altitude z. The intelligent PSO and the cooperative PSO-CS drive the system toward perfect tracking of the prescribed square trajectory. Also, all these methods allowed maintaining the yaw angle in zero degrees: no overturning of the aircraft. Figure 10 and Figure 11 confirm that these intelligent approaches of PSO-CS nd PSO give controlled responses closed to the desired roll and pitch attitudes.

To test the robustness of these PID controllers, external disturbance was introduced. The most likely disturbance acting on the quadrotor is the wind disturbance along z axis. The external disturbance is modeled in this simulation by a random signal of variable amplitude inserted in the closed loop scheme of altitude z. It is noticed from Figure 12, Figure 13 and Figure 14 and the errors ISE presented in

Table 6. Errors (ISE-Integral Square Error) obtained Particle Swarm Optimization (PSO), Cuckoo Search (CS), Particle Swarm Optimization-Cuckoo Search (PSO-CS), and Reference Model (RM) for the quadrotor with disturbance (rectangular trajectory).

Method	Errors ISE for Output
	x	y	z	ϕ	θ	ψ
PSO	1.254	1.161	14.212	0.524	0.982	1.056 × 10−28
CS	1.498	1.334	21.612	0.147	2.441	7.378 × 10−35
PSO-CS	1.079	0.685	3.287	0.183	0.195	8.157 × 10−36
RM	2.555	1.573	33.776	0.272	0.704	5.467 × 10−35

that when random disturbance affects the quadrotor, PSO, CS, and PSO-CS controllers reject these undesirable effects and provide outputs responses close to the desired ones. Therefore, we can affirm that PSO, CS, and especially the cooperative PSO-CS, are able to tune optimally the six PID controllers and drive the quadrotor to the desired altitude z, positions x and y, and yaw angle.

5.3. Helical Trajectory Tracking

Through the last simulation, we will investigate more the efficiency of the proposed intelligent control methods (PSO, CS, and PSO-CS) to search for optimal PID autotuned controllers for the quadrotor, when the desired trajectory is a helix. Simulation results of each position’s tracking (x, y, z), and the tracking in 3D plan are presented in Figure 15. The tracking of roll and pitch attitudes (ϕ and θ) are presented in Figure 16 and Figure 17.

Table 7. Errors (ISE-Integral Square Error) obtained with Particle Swarm Optimization (PSO), Cucko Search (CS), Particle Swarm Optimization-Cuckoo Search (PSO-CS), and Reference Model (RM) for the quadrotor (helical trajectory).

Method	ISE Errors for Output
	x	y	z	ϕ	θ	ψ
PSO	0.006	0.013	0.069	0.010	0.069	7.540 × 10−38
CS	0.013	0.005	0.121	0.023	0.089	6.982 × 10−36
PSO-CS	0.011	0.010	0.007	0.001	0.033	7.495 × 10−35
RM	0.006	0.004	0.08	0.238	1.565	1.954 × 10−34

summarizes the obtained ISE errors.

Figure 15 clearly shows the possibility of the quadrotor system to follow a linear movement in altitude z in the first mode from 0 to 18 s and to keep perfectly a stable attitude z of 20 m in the second mode from 18 to 30 s. For translation motion, the quadrotor makes a complete circle of 2 m of radius in 30 s, where we notice that the errors in x, y, and z axis are very small (from Table 7). The same figure illustrates in 3D plan that the quadrotor tracks the helix trajectory and ends by a circle form. Also, the quadrotor avoids with high precision the overturning and maintains exactly 0° in yaw attitude.

Figure 16 and Figure 17 and Table 7 illustrate that PSO, CS, and PSO-CS help to follow perfectly the sensor signals (ϕ_d and θ_d). These results prove the possibility of our system to follow the imposed combined trajectory form with precision, for all variables (positions (x and y), altitude z, and attitudes (ϕ, θ, and ψ).

These controllers are tested in order to evaluate their performances under the same disturbance (a random input of variable amplitude) inserted in the altitude z closed loop scheme of the quadrotor. It can be observed from Figure 18, Figure 19 and Figure 20 and

Table 8. Errors (ISE-Integral Square Error) obtained with Particle Swarm Optimization (PSO), Cuckoo Search (CS), Particle Swarm Optimization-Cuckoo Search (PSO-CS), and Reference Model (RM) for the quadrotor with disturbance (helical trajectory).

Method	ISE Errors for Output
	x	y	z	ϕ	θ	ψ
PSO	0.072	0.048	0.163	1.079	2.606	2.001 × 10−34
CS	0.756	0.188	0.172	1.155	6.084	3.721 × 10−36
PSO-CS	0.009	0.009	0.007	0.010	0.023	1.128 × 10−35
RM	55.221	53.052	0.932	2.016	6.745	2.291 × 10−33

that PSO and PSO-CS controllers reject the disturbance effects and drive all variables to the desired ones. Therefore, we can affirm that despite the usage of a combined trajectory form and the presence of forcible disturbance, adaptive PSO-CS and PSO are proven as remarkably robust and efficient for quadrotor control.

6. Conclusions

In this paper, Particle Swarm Optimization (PSO), Cuckoo Search (CS), and the proposed cooperative PSO-CS were applied to control the quadrotor system. Results proved the efficiency of PSO approach and the cooperative PSO-CS to control optimally the quadrotor’s outputs (ϕ, θ, ψ, x, y, z) compared to CS and the classical Reference Model (RM) methods. Experimental simulations also confirm the effectiveness of the aircraft using the proposed intelligent PID controllers to track the imposed rectangular and helical trajectories. Indeed, the robustness and performances of the intelligent controllers are validated in presence of disturbance. The performances of PSO-CS showed not only the stability robustness against the external disturbance, but also the fast response and the excellent tracking capacity for the different trajectories. These achievements are the results of the good settings of PSO, CS, and PSO-CS parameters. The advantages of the proposed PSO-CS in seeking for the best solutions in control problem are the results of exploiting the local search capacity of CS and benefiting from the global intelligence offered by PSO.

Concerning the limitations of our PSO-CS, one more parameter is added to CS program, which is the acceleration constant C₂. A wrong setting of this parameter can cause premature convergence. C₂ makes the cuckoos move toward the global best solution and must ensure a compromise between local and global exploration of the search space. On the other hand, the PSO-CS program consumes additional CPU time to calculate the velocity update equation over generations in order to perform the search.