Model references frequency-domain design based fractional order PD controller for quadrotor control: attitude, altitude, and position studies

Abstract

This article proposes a control strategy based on a Fractional Order Proportional Derivative (FOPD) controller to enhance the performance and robustness of quadrotors operating in trajectory tracking mode. The tuning method proposed for the FOPD controller is based on an ideal Bode transfer function in a closed loop as a fractional reference model. The latter represents the quadrotor desired closed-loop response, and its specification is the key design decision. In this study, we first appropriately set model reference parameters to realistic desired time response characteristics (i.e. overshoot and time constant). Then, the parameters of the FOPD controller are computed using the direct synthesis design approach. The dynamic model of the quadrotor is developed using Lagrange formulation. The model is used to test the performance of the proposed FOPD-based control strategy for the attitude, position, and altitude of the flying robot. The results show that the proposed controller achieved accurately set points tracking with transient responses (i.e. overshoot and time constant) that satisfactorily meet those assigned for their respective reference models. These results validate the effectiveness of our proposed control strategy. Moreover, our comparative studies of the proposed FOPD controller with the conventional PD controller show that the former satisfies better-desired design requirements.

Keywords:

• Quadrotor–Lagrange approach
• fractional order PID controller
• ideal transfer function of Bode
• conventional PID controller

1. Introduction

Quadrotors, also known as Unmanned Aerial Vehicles (UAVs), are vehicles capable of flying missions without a human presence on board. Quadrotors can execute a broad range of applications without placing human life at risk and with low expenses. Nowadays, these vehicles are used in agricultural services, wildfire surveillance, search and rescue missions (Kashino, Nejat, & Benhabib, 2019), power plants inspection and object detection (Vaddi, 2019), mapping and photographing in indoor and outdoor environments (Behrooz, Mariun, Marhaban, Mohd Radzi, & Ramli, 2018; Falanga, Kleber, Mintchev, Floreano, & Scaramuzza, 2019), marine operations, battle damage assessment, border interdiction prevention, and many other applications. Quadrotor’s merits are its design simplicity, small size, good maneuverability, high agility, vertical take-off landing ability, and an increased payload. The quadrotor has simple mechanical control linkages because it relies on fixed-pitch rotors and for vehicle control, it uses the variation in motor speed (Idrissi & Annaz, 2020). However, these merits come at an expense as controlling a quadrotor is not easy because of the coupled dynamics and its commonly underactuated design configuration (Idrissi, Annaz, & Salami, 2021). In addition, the dynamics of the quadrotor are highly non-linear and several uncertainties are encountered during its missions (Idrissi, Salami, & Annaz, 2021). It has been revealed that controlling quadrotors is still a major challenge to be considered to further improve their performance despite unmodelled dynamic factors and external disturbances. This led to the motivation of researching designing control algorithms for these flying machines to be able to maneuver in harsh environments and cope with difficult tasks.

Several popular control strategies and their combination are explored in the literature (Idrissi, Salami, & Annaz, 2022; Mustafa Abro, Ali, Zulkifli, & Asirvadam, 2021; Roy et al., 2021; Zulu & John, 2014). These include classical PID, Linear quadratic regulator (LQR), Sliding mode control (SMC), Backstepping control, Adaptive control, and artificial intelligence algorithm as neural and fuzzy control. In Idrissi, Salami, et al. (2021) and Khalid, Zeb, and Haider (2019), a conventional PID controller was applied to control the attitude of a quadrotor as well as in Benbouali, Taleb, Chabni, and Bessaad (2016) and Jiao, Liu, Zhang, and Lian (2018) where the same controller is extended to control a position and orientation. These studies revealed that applying a conventional PID controller to the quadrotor limits its performance and the tuning of the PID controller, as a conventional type, could cause some challenges in achieving trajectory tracking according to desired specifications and requirements. In Bouabdallah, Noth, and Siegwart (2004) and Khatoon, Gupta, and Das (2014), a comparative study between PID and LQR control is made and the results concluded that both strategies provide average results. In Everett (2015), a linear quadratic regulator (LQR) combined with Integral Feedback was applied to a Parrot Minidrone. They assure that this technique is unsuitable for indoor navigation which requires consistent minimum steady-state error. A sliding mode control (SMC) is also applied, in Khalid et al. (2019), by dividing the quadrotor system into the full-actuated and under-actuated systems this study gets good results of stability and robustness but the chattering effect of SMC was present. Good tracking results for both position and attitude control were achieved when the authors of Karabacak and Eskikurt (2012) used a backstepping strategy to control our engine; However, it was difficult to implement and it wasn’t robust against the external perturbation. To increase the robustness, several algorithms of robust control were designed. In Liu, Bai, Lu, and Zhong (2013), for example, attitude control is designed based on linear control and robust compensation. The results present good robustness but a poor tracking trajectory.

On another side, artificial intelligent control systems are also used to control aerial robots. The fuzzy controller was used by (Benbouali et al., July 2021; Behrooz et al., 2018), which gives a good response, but the tuning approach of control variables presents the limitation of the method. The neural network is also presented by the authors of Razmi and Afshinfar (2019) and Jin, He, Wu, Wang, and Chi (2020), they give good results in tracking and robustness properties; However, these results were achieved with expensive control law. Among all the control strategies proposed above, none of them presents the best of the desired features. However, control strategies based on PID are by the far the most superior in terms of design simplicity and computation time, it is indeed the most applied controller in the industry (Hui-Fang, Qiu-Sheng, Zhi-Cheng, & Jing-Gang, 2015).

PID controller, however, presents some shortcomings such as sensitivity to system parameter changes, performance deterioration as the order of the system to be controlled increases, and poor performance with overshooting occurrence in nonlinear systems. To overcome these shortcomings, a generalized version of the conventional PID controllers, known as fractional-order PID (FO-PID) has recently been proposed (Chopade, Khubalkar, Junghare, Aware, & Das, 2018). FO-PID controllers have been shown to provide better performance and robustness than conventional PID controllers (Li, Ding, Wu, & Lin, 2021). The most important features of FO-PID controllers are more flexibility, better performance with high-order systems, and better response in the case of nonminimum phase systems (Leena, Sundaravadivu, Monisha, & Rajinikanth, 2018; Liu, Zhang, Pan, & Zhang, 2022; Tepljakov et al., 2021; Valluru, Kumar, & Ayush, 2019; Yongqiang, 2017; Zheng, Chen, Wang, Lin, & Guo, 2023; Zheng, Luo, Chen, & Wang, 2021; Zhuo-Yun, Yi-Min, Qing-Guo, Rui-Juan, & Lei-Jun, 2020). These advantages attracted the attention of researchers to apply fractional calculus in the field of control design.

In this work, the focus is on the use of a fractional-order PID controller (FOPID) and its tuning for the control of quadrotor motion. Several tuning methods for setting the parameters of FO-PID controllers have been developed (Chopade et al., 2018; Leena et al., 2018). It is claimed, that FO-PID is expected to be the best solution to overcome the control problems for attitude and position tracking and stabilization of the quadrotor. In this direction, the objective of this work is to investigate the performance of a control scheme based on a FOPD controller for the position, attitude, and altitude control of the quadrotor. The tuning method proposed for FOPD parameters is based on an ideal transfer function of Bode in a closed-loop as a reference model (Zheng et al., 2021; Zhuo-Yun et al., 2020) and is to be discussed in Sec. 3. The simplicity and effectiveness of the tuning method motivated us to investigate its performance to tune a FOPD controller for the control of the quadrotor. Our contribution is to demonstrate, via extensive simulation, that the proposed tuning method of the FOPD controller gives better performance than a pole placement tuned PD controller in terms of tracking accuracy, transient response, and robustness for controlling a quadrotor. The results obtained using the proposed method are very satisfactory and meet the desired specifications. Moreover, a comparative study is performed between the FO-PD and PD performances under the same operating conditions. It is shown that the FOPD is able to meet desired specifications better than PD controller.

This paper is organized as follows: Description and mathematical modeling of the quadrotor are addressed in Sec. 2. The configuration of the fractional PID controller and the proposed tuning method are detailed in Sec. 3. Simulation results for attitude and position control along with comparative studies with conventional PID controllers are presented in Sec. 4. Finally, the conclusion is provided in Sec. 5.

2. Dynamic modeling of the quadrotor

Several models of the quadrotor are developed in the literature to design a controller for quadrotors. In Idrissi, Annaz, et al. (2021), Idrissi, Salami, et al. (2021) and Idrissi and Annaz (2020) the authors derived a model based on selected physical parameters. A study in (Yongqiang, 2017) preferred to use a model obtained using experimental data and system identification techniques. The dynamic model proposed in this work is based on Lagrange formulation (Benbouali et al., 2016, 2021; Halcı, Gazi, & Cihan, 2019). A detailed description of the dynamic model developed in this work can be found in the reference (Benbouali et al., 2016). As shown in Figure 1, a quadrotor is an aerial robot with four rotors placed at the ends of a cross-form rigid body. The inertial and body-fixed frames of the quadrotor are defined as ($\overrightarrow{W_1};\overrightarrow{W_2};\overrightarrow{W_3}$) and ($\overrightarrow{B_1};\overrightarrow{B_2};\overrightarrow{B_3}$) respectively.

Figure 1. Quadrotor axis system.

The angular positions (attitude): the yaw angle $\varphi ,$ pitch angle $\theta ,$ and roll angle $\psi$ are defined respectively as follows:

1. Roll angle $\varphi$ determines the rotation of the quadrotor around $\overrightarrow{B_1}$

2. Pitch angle $\theta$ determines the rotation of the quadrotor around $\overrightarrow{B_2}$

3. Yaw angle $\psi$ determines the rotation of the quadrotor around $\overrightarrow{B_2}$

The operation of a quadrotor is quite particular. By ingeniously varying the power of the rotors’ motors, it is possible to raise/lower, to tilt it according to its three angular positions (roll, pitch, and angle yaw) and thus control its position (x, y, and z).

The dynamic model is derived under the following assumptions:

1. The structure is rigid and symmetrical.

2. The center of mass of the vehicle and the body-fixed frame origin is assumed to coincide.

3. The propellers are rigid in a plane.

Based on the above assumptions, the model of the quadrotor is developed using the Lagrange approach (Benbouali et al., 2016, 2021). (1) $$\{\begin{array}{l}{\Gamma }_i=\frac{d}{dt}\left(\frac{\partial L}{\partial q_i}\right)-\frac{\partial L}{\partial q_i}\\ L=T-V\end{array}$$(1)

Where, generalized coordinates, generalized forces; Total kinetic energy; V: Total potential energy.

As shown in (1)(1) $$\{\begin{array}{l}{\Gamma }_i=\frac{d}{dt}\left(\frac{\partial L}{\partial q_i}\right)-\frac{\partial L}{\partial q_i}\\ L=T-V\end{array}$$(1) , the motion equation is given by Lagrangian $L=T-V.$ Thus, the motion equations of the quadrotor can be expressed as: (2) $$\{\begin{array}{l}{\Gamma }_{\varphi }=\frac{d}{dt}\left(\frac{\partial L}{\partial \varphi }\right)-\frac{\partial L}{\partial \varphi }={\tau }_X\\ {\Gamma }_{\theta }=\frac{d}{dt}\left(\frac{\partial L}{\partial \theta }\right)-\frac{\partial L}{\partial \theta }={\tau }_Y\\ {\Gamma }_{\psi }=\frac{d}{dt}\left(\frac{\partial L}{\partial \psi }\right)-\frac{\partial L}{\partial \psi }={\tau }_Z\end{array}$$(2) where $q_i$ represent [$\varphi ,$ $\theta ,$ $\psi$] angles.

The kinetic energy is expressed as $T=\frac{1}{2}m\cdot v^2;$

where m is the mass of the quadrotor and v is the speed of the quadrotor which is given by derivate the position vector $r_{W1,W2,W3}(x,y,z).$

we note that: (3) $$r_{W1,W2,W3}\left(x,y,z\right)=R\left(\varphi ,\theta ,\psi \right)\cdot \left(\begin{array}{c}x\\ y\\ z\end{array}\right)$$(3) $R(\varphi ,\theta ,\psi )$ represents the rotational matrix and it is defined as follows: (4) $$R\left(\varphi ,\theta ,\psi \right)=\left(\begin{array}{ccc}{C}_{\psi }{C}_{\theta }& C_{\psi }{S}_{\theta }{S}_{\varphi }-S_{\psi }{C}_{\varphi }& C_{\psi }{S}_{\theta }{C}_{\varphi }+S_{\psi }{S}_{\varphi }\\ S_{\psi }{C}_{\theta }& S_{\psi }{S}_{\theta }{S}_{\varphi }+C_{\psi }{C}_{\varphi }& S_{\psi }{S}_{\theta }{C}_{\varphi }-S_{\varphi }{C}_{\psi }\\ -S_{\theta }& C_{\theta }{S}_{\varphi }& C_{\theta }{C}_{\varphi }\end{array}\right)$$(4) where $S_i:\hspace{0.17em}\sin \hspace{0.17em}(i);\begin{array}{cc}& \end{array}{C}_i:\hspace{0.17em}\cos \hspace{0.17em}(i)$

After simplification, we obtain the following expression: (5) $$\begin{array}{c}T=\frac{1}{2}{I}_x{(\varphi -\psi S_{\theta })}^2+\frac{1}{2}{I}_y{(\theta C_{\varphi }+\psi S_{\varphi }{C}_{\theta })}^2\\ +\frac{1}{2}{I}_Z{(\theta S_{\varphi }-\psi C_{\varphi }{C}_{\theta })}^2\end{array}$$(5) where I_x, I_y, and I_z are the moments of inertia in direction of x, y, and z-axis.

The potential energy V is expressed by: (6) $$V=\int \mathit{xdm}(-g\cdot S_{\theta })+\int \mathit{ydm}(g\cdot S_{\varphi }{C}_{\theta })+\int \mathit{zdm}(g\cdot C_{\varphi }{C}_{\theta })$$(6)

Using Eqs. (5)(5) $$\begin{array}{c}T=\frac{1}{2}{I}_x{(\varphi -\psi S_{\theta })}^2+\frac{1}{2}{I}_y{(\theta C_{\varphi }+\psi S_{\varphi }{C}_{\theta })}^2\\ +\frac{1}{2}{I}_Z{(\theta S_{\varphi }-\psi C_{\varphi }{C}_{\theta })}^2\end{array}$$(5) and (6)(6) $$V=\int \mathit{xdm}(-g\cdot S_{\theta })+\int \mathit{ydm}(g\cdot S_{\varphi }{C}_{\theta })+\int \mathit{zdm}(g\cdot C_{\varphi }{C}_{\theta })$$(6) , the motion equations of the quadrotor given in (2)(2) $$\{\begin{array}{l}{\Gamma }_{\varphi }=\frac{d}{dt}\left(\frac{\partial L}{\partial \varphi }\right)-\frac{\partial L}{\partial \varphi }={\tau }_X\\ {\Gamma }_{\theta }=\frac{d}{dt}\left(\frac{\partial L}{\partial \theta }\right)-\frac{\partial L}{\partial \theta }={\tau }_Y\\ {\Gamma }_{\psi }=\frac{d}{dt}\left(\frac{\partial L}{\partial \psi }\right)-\frac{\partial L}{\partial \psi }={\tau }_Z\end{array}$$(2) can be re-expressed as: (7) $$\{\begin{array}{l}{\Gamma }_{\varphi }=I_x{\stackrel{•}{\omega }}_x+(I_z-I_y){\omega }_y{\omega }_z={\tau }_X\\ {\Gamma }_{\theta }=I_y{\stackrel{•}{\omega }}_y+(I_x-I_z){\omega }_x{\omega }_z={\tau }_Y\\ {\Gamma }_{\psi }=I_z{\stackrel{•}{\omega }}_z+(I_y-I_x){\omega }_x{\omega }_y={\tau }_Z\end{array}$$(7)

All the effects of the forces applied to the quadrotor such as Thrust, Drag, and Gyroscopic effects are expressed in the system below: (8) $$\{\begin{array}{l}{\tau }_x=b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)+\\ I_r{\omega }_y({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)\\ {\tau }_y=b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)+\\ I_r{\omega }_x(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)\\ {\tau }_z=d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)\end{array}$$(8) where b: Thrust coefficient, d: Drag coefficient, I_r: Rotor Inertia, Angular velocity of motor I, l: Distance between the axis of rotation of the rotor and the center of the mass, Speed in the fixed frame.

The speed $\omega$ in the fixed frame is given by: (9) $$\left(\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right)=\left(\begin{array}{c}\stackrel{•}{\varphi }-\stackrel{•}{\psi }\hspace{0.17em}\sin \hspace{0.17em}(\theta )\\ \stackrel{•}{\theta }\hspace{0.17em}\cos \hspace{0.17em}(\varphi )+\stackrel{•}{\psi }\hspace{0.17em}\sin \hspace{0.17em}(\varphi )\hspace{0.17em}\cos \hspace{0.17em}(\theta )\\ \stackrel{•}{\theta }\hspace{0.17em}\sin \hspace{0.17em}(\varphi )-\stackrel{•}{\psi }\hspace{0.17em}\cos \hspace{0.17em}(\varphi )\hspace{0.17em}\cos \hspace{0.17em}(\theta )\end{array}\right)$$(9)

By substituting (8)(8) $$\{\begin{array}{l}{\tau }_x=b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)+\\ I_r{\omega }_y({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)\\ {\tau }_y=b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)+\\ I_r{\omega }_x(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)\\ {\tau }_z=d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)\end{array}$$(8) and (9)(9) $$\left(\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right)=\left(\begin{array}{c}\stackrel{•}{\varphi }-\stackrel{•}{\psi }\hspace{0.17em}\sin \hspace{0.17em}(\theta )\\ \stackrel{•}{\theta }\hspace{0.17em}\cos \hspace{0.17em}(\varphi )+\stackrel{•}{\psi }\hspace{0.17em}\sin \hspace{0.17em}(\varphi )\hspace{0.17em}\cos \hspace{0.17em}(\theta )\\ \stackrel{•}{\theta }\hspace{0.17em}\sin \hspace{0.17em}(\varphi )-\stackrel{•}{\psi }\hspace{0.17em}\cos \hspace{0.17em}(\varphi )\hspace{0.17em}\cos \hspace{0.17em}(\theta )\end{array}\right)$$(9) in the system (7)(7) $$\{\begin{array}{l}{\Gamma }_{\varphi }=I_x{\stackrel{•}{\omega }}_x+(I_z-I_y){\omega }_y{\omega }_z={\tau }_X\\ {\Gamma }_{\theta }=I_y{\stackrel{•}{\omega }}_y+(I_x-I_z){\omega }_x{\omega }_z={\tau }_Y\\ {\Gamma }_{\psi }=I_z{\stackrel{•}{\omega }}_z+(I_y-I_x){\omega }_x{\omega }_y={\tau }_Z\end{array}$$(7) , and assuming small angles deviations, the final attitude dynamic model is derived as follows: (10) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{I_r\stackrel{•}{\theta }({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)}{I_x}+\frac{(I_y-I_z)}{I_x}\stackrel{•}{\theta }\stackrel{•}{\psi }\\ +\frac{b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_x}\\ \stackrel{••}{\theta }=\frac{I_r\stackrel{•}{\varphi }(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)}{I_y}+\frac{(I_z-I_x)}{I_y}\stackrel{•}{\varphi }\stackrel{•}{\psi }+\\ \frac{b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)}{I_y}\\ \stackrel{••}{\psi }=\frac{(I_x-I_y)}{I_z}\stackrel{•}{\theta }\stackrel{•}{\varphi }+\frac{d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_z}\end{array}$$(10)

As mentioned above, this model is given in the direction of the axis shown in Figure 1.

The dynamics of the translation model are developed based Newton-Euler equation (Benbouali et al., 2016). (11) $$m\overrightarrow{a}=\overrightarrow{P}+{\sum }_i\hspace{0.17em}\overrightarrow{F_i}$$(11) with $\overrightarrow{P}$ the weight, ${\sum }_i\hspace{0.17em}\overrightarrow{F_i}$ represented the sum of the i^th thrust.

Assume the following assumption for the thrust: (12) $${\sum }_i\hspace{0.17em}\overrightarrow{F_i}=b({\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)\overrightarrow{z}={\sum }_i(F_i\overrightarrow{z})$$(12)

Hence, the commonly used translation dynamic model (Benbouali et al., 2016; Bouabdallah et al., 2004) is derived as: (13) $$m\left(\begin{array}{c}\stackrel{\cdot \cdot }{X}\\ \stackrel{\cdot \cdot }{Y}\\ \stackrel{\cdot \cdot }{Z}\end{array}\right)=mg\left(\begin{array}{c}0\\ 0\\ -1\end{array}\right)+{\sum }_i\hspace{0.17em}\overrightarrow{F_i}\left(\begin{array}{c}c(\psi )s(\theta )c(\varphi )+s(\psi )s(\varphi )\\ s(\psi )s(\theta )c(\varphi )-c(\psi )s(\varphi )\\ c(\theta )c(\varphi )\end{array}\right)$$(13)

Finally, the complete quadrotor dynamic model derived from Euler-Lagrange and Newton-Euler with constraints on angles variations is given in Eq. (14)(14) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{I_r\stackrel{•}{\theta }({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)}{I_x}+\frac{(I_y-I_z)}{I_x}\stackrel{•}{\theta }\stackrel{•}{\psi }\\ +\frac{b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_x}\\ \stackrel{••}{\theta }=\frac{I_r\stackrel{•}{\varphi }(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)}{I_y}+\frac{(I_z-I_x)}{I_y}\stackrel{•}{\varphi }\stackrel{•}{\psi }+\\ \frac{b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)}{I_y}\\ \stackrel{••}{\psi }=\frac{(I_x-I_y)}{I_z}\stackrel{•}{\theta }\stackrel{•}{\varphi }+\frac{d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_z}\\ \stackrel{••}{X}=\frac{(\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi +\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Y}=\frac{(\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi -\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Z}=-g+\frac{(\hspace{0.17em}\cos \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\end{array}$$(14) below. This model will be simulated in MATLAB to test and evaluate the performance of the FOPID controller proposed. (14) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{I_r\stackrel{•}{\theta }({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)}{I_x}+\frac{(I_y-I_z)}{I_x}\stackrel{•}{\theta }\stackrel{•}{\psi }\\ +\frac{b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_x}\\ \stackrel{••}{\theta }=\frac{I_r\stackrel{•}{\varphi }(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)}{I_y}+\frac{(I_z-I_x)}{I_y}\stackrel{•}{\varphi }\stackrel{•}{\psi }+\\ \frac{b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)}{I_y}\\ \stackrel{••}{\psi }=\frac{(I_x-I_y)}{I_z}\stackrel{•}{\theta }\stackrel{•}{\varphi }+\frac{d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_z}\\ \stackrel{••}{X}=\frac{(\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi +\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Y}=\frac{(\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi -\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Z}=-g+\frac{(\hspace{0.17em}\cos \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\end{array}$$(14)

The model parameters and their values used in the simulation are indicated in Table 1.

Table 1. Model parameters.

Symbol	Unit	Values
b	Kg.m/rad^2	54.2 × 10^-6
d	Kg.m^2/rad^2	1.1 × 10^-6
g	m/s^2	9.81
l	m	0.24
m	kg	1
Ix	Kg.m^2	8.1 × 10^-3
Iy	Kg.m^2	8.1 × 10^-3
Iz	Kg.m^2	14.2 × 10^-3
Jrotor	Kg.m^2	104 × 10^-6
km	N.m/A	6.3 × 10^-3
R	Ω	0.6
r	-	5.6
η	-	0.9
ΩH	Rad/s	115

3. Fractional order PID controller

Fractional order proportional-integral-derivative (PI^λD^μ) controller control is a generalization of the classical PID controller. It is introduced in control applications to improve robustness and control performance. The distinguishing feature of PI^λD^μ is its potential to realize better control performance over the classical PID controller. It has extra two parameters that is, integral order λ and differential order μ which are adjustable and hence the controller flexibility is increased.

The function of transfer of fractional-order PI^λD^μ can be expressed as follows (Liu et al., 2022; Valluru et al., 2019; Yumuk, Güzelkaya, & Eksin, 2022; Zhuo-Yun et al., 2020): (15) $$C\left(p\right)=K_p+\frac{K_i}{p^{\lambda }}+K_d{p}^{\mu }$$(15) where the five tuning parameters are: K_p proportional gain, K_i integral gain, and K_d derivative gain; λ and μ are as defined above.

Figure 2 shows the internal structure of the fractional-order PI^λD^μ controller. It consists of proportional, fractional order integral, and fractional-order derivative control actions placed in parallel (Liu et al., 2022; Valluru et al., 2019; Yumuk et al., 2022; Zhuo-Yun et al., 2020):

Figure 2. Structure of the fractional-order PI^λD^μ controller.

Currently, many researchers are interested in the tuning of the fractional-order PI^λD^μ controllers. A lot of work is achieved in developing a suitable and effective tuning method for fractional order PI^λD^μ controller (Barbosa, Tenreiro Machado, & Ferreira, 2004; Leena et al., 2018; Valério & Sà da Costa, 2006; Zhuo-Yun et al., 2020). The methods developed can be classified as time-domain methods, frequency-domain methods, and comprehensive methods. In (Halcı et al., 2019; Li et al., 2021) a tuning method inspired by the tuning method developed in (Zheng et al., 2021; Zhuo-Yun et al., 2020) for conventional PID controller is proposed. The parameters of the fractional-order PI^λD^μ controller are adjusted so that the closed-loop control system behaves like a desired reference model which meets desired closed-loop specifications.

As the starting point for the design, consider the transfer function of the closed-loop system of the form: (16) $$G_A\left(p\right)=\frac{C\left(p\right)G_p\left(p\right)}{1+C\left(p\right)G_p\left(p\right)}$$(16) where C(p) is the transfer function of the controller PI^λD^μ defined by (15)(15) $$C\left(p\right)=K_p+\frac{K_i}{p^{\lambda }}+K_d{p}^{\mu }$$(15) and G_p(p) is the open-loop function of the system to be controlled.

The setting of the five parameters of the controller C(p) must ensure that the time response of the closed-loop system with PI^λD^μ is the same as that of the reference model. The specification of this latter is the key design decision (Zheng et al., 2021; Zhuo-Yun et al., 2020).

In this work, Bode’s ideal transfer function in a closed-loop is chosen as a reference model. This choice is justified due to the following features: it is fractional order, makes the closed-loop system robust to gain variations and the step responses exhibit an over-damping property (i.e. overshoot of closed-loop system remains unchanged for different values of controller gains).

The ideal Bode’s transfer function in closed-loop is defined by: (17) $$G_d\left(p\right)=\frac{1}{1+{\left(\frac{p}{{\omega }_u}\right)}^m}$$(17) with the tuning parameters (m, ω_c) are the order and the crossover frequency, respectively. The order m may assume a real no integer value between 1 and 2. The parameters m and ${\omega }_u$ are chosen so that the reference system meets the desired specifications.

Hence, the purpose of this design approach is simply to find the five parameters of the fractional PI^λD^μ that satisfy the following relation: $G_A(p)\hspace{0.17em}≌\hspace{0.17em}{G}_d\left(p\right).$

The transfer functions $G_A(p)$ and $G_d(p)$ can be developed in the Taylor-Maclaurin series around the unity gain cross-over frequency ω_u as follows: (18) $$G_A\left(p\right)=G_A\left({\omega }_u\right)+\left(p-{\omega }_u\right){G_A}^{\left(1\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^2}{2!}{G_A}^{\left(2\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^3}{3!}{G_A}^{\left(3\right)}\left({\omega }_u\right)+\dots +\frac{{\left(p-{\omega }_u\right)}^n}{n!}{G_A}^{\left(n\right)}\left({\omega }_u\right)+\dots$$(18) (19) $$G_d\left(p\right)=G_d\left({\omega }_u\right)+\left(p-{\omega }_u\right){G_d}^{\left(1\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^2}{2!}{G_d}^{\left(2\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^3}{3!}{G_d}^{\left(3\right)}\left({\omega }_u\right)+\dots +\frac{{\left(p-{\omega }_u\right)}^n}{n!}{G_d}^{\left(n\right)}\left({\omega }_u\right)+\dots$$(19) where ${G_A}^{\left(n\right)}\left({\omega }_u\right)$ and ${G_d}^{\left(n\right)}\left({\omega }_u\right)$ are respectively the n^ths derived of the functions $G_A(p)$ and $G_d(p)$ with respect to the variable p.

Since the design objective find five parameters of FOPID that satisfy the condition $G_A(p)\hspace{0.17em}≌\hspace{0.17em}{G}_d\left(p\right).$ Comparing their Taylor-Maclaurin series (18)(18) $$G_A\left(p\right)=G_A\left({\omega }_u\right)+\left(p-{\omega }_u\right){G_A}^{\left(1\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^2}{2!}{G_A}^{\left(2\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^3}{3!}{G_A}^{\left(3\right)}\left({\omega }_u\right)+\dots +\frac{{\left(p-{\omega }_u\right)}^n}{n!}{G_A}^{\left(n\right)}\left({\omega }_u\right)+\dots$$(18) and (19)(19) $$G_d\left(p\right)=G_d\left({\omega }_u\right)+\left(p-{\omega }_u\right){G_d}^{\left(1\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^2}{2!}{G_d}^{\left(2\right)}\left({\omega }_u\right)+\frac{{(p-{\omega }_u)}^3}{3!}{G_d}^{\left(3\right)}\left({\omega }_u\right)+\dots +\frac{{\left(p-{\omega }_u\right)}^n}{n!}{G_d}^{\left(n\right)}\left({\omega }_u\right)+\dots$$(19) equations, the design is changed to the following five equations to satisfy: (20) $$\{\begin{array}{l}{\mathrm{ }G}_A\left({\omega }_u\right)=G_d\left({\omega }_u\right)\\ \begin{array}{c}{G_A}^{\left(1\right)}\left({\omega }_u\right)={G_d}^{\left(1\right)}\left({\omega }_u\right)\\ {G_A}^{\left(2\right)}\left({\omega }_u\right)={G_d}^{\left(2\right)}\left({\omega }_u\right)\end{array}\\ \begin{array}{c}{G_A}^{\left(3\right)}\left({\omega }_u\right)={G_d}^{\left(3\right)}\left({\omega }_u\right)\\ {G_A}^{\left(4\right)}\left({\omega }_u\right)={G_d}^{\left(4\right)}\left({\omega }_u\right)\end{array}\end{array}$$(20)

Hence, the five parameters of the fractional-order PI^λD^μ controller can be obtained by satisfying and solving the above five equations. To determine the parameters of the fractional order PD controller, Eq. (15)(15) $$C\left(p\right)=K_p+\frac{K_i}{p^{\lambda }}+K_d{p}^{\mu }$$(15) is replaced by the following expression: $C\left(p\right)=K_p+K_d{p}^{\mu }.$ This results in reducing equations to solve in (20)(20) $$\{\begin{array}{l}{\mathrm{ }G}_A\left({\omega }_u\right)=G_d\left({\omega }_u\right)\\ \begin{array}{c}{G_A}^{\left(1\right)}\left({\omega }_u\right)={G_d}^{\left(1\right)}\left({\omega }_u\right)\\ {G_A}^{\left(2\right)}\left({\omega }_u\right)={G_d}^{\left(2\right)}\left({\omega }_u\right)\end{array}\\ \begin{array}{c}{G_A}^{\left(3\right)}\left({\omega }_u\right)={G_d}^{\left(3\right)}\left({\omega }_u\right)\\ {G_A}^{\left(4\right)}\left({\omega }_u\right)={G_d}^{\left(4\right)}\left({\omega }_u\right)\end{array}\end{array}$$(20) to three instead of five. For PD controller the parameters are obtained using pole assignment technique. For comparative purposes, the PD controller parameters (K_p and K_d) are set to meet the same closed loop time response desired specifications (overshoot and settling time) as FOPD tuned using the proposed approach.

4. Quadrotor control strategy

As mentioned above, the tuning of the control algorithm proposed is based on the linear system model transfer function and a desired closed-loop reference model. However, the quadrotor dynamic model is nonlinear and quite complex and hence it cannot be used for the design of the controller. To simplify the control design procedure, simple transfer functions models to represent the dynamic of the quadrotor are needed. Since the model is valid only for small angles variations then the quadrotor can be assumed close to the stationary state. Hence, taking into account these latter assumptions we have:

1. $\ddot{\mathrm{z}}\approx 0$

2. $\cos \hspace{0.17em}\mathrm{\varphi }=\cos \hspace{0.17em}\mathrm{\theta }\approx 1$

3. $\sin \hspace{0.17em}\mathrm{\varphi }\approx \mathrm{\varphi }$ and $\sin \hspace{0.17em}\mathrm{\theta }\approx \mathrm{\theta }$

4. $\mathrm{\psi }$ = 0

5. ${\mathrm{\Omega }}_1-{\mathrm{\Omega }}_2+{\mathrm{\Omega }}_3-{\mathrm{\Omega }}_4$ are neglected because they are

smaller than the main ones. Using the above assumptions, the equations system (14)(14) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{I_r\stackrel{•}{\theta }({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)}{I_x}+\frac{(I_y-I_z)}{I_x}\stackrel{•}{\theta }\stackrel{•}{\psi }\\ +\frac{b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_x}\\ \stackrel{••}{\theta }=\frac{I_r\stackrel{•}{\varphi }(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)}{I_y}+\frac{(I_z-I_x)}{I_y}\stackrel{•}{\varphi }\stackrel{•}{\psi }+\\ \frac{b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)}{I_y}\\ \stackrel{••}{\psi }=\frac{(I_x-I_y)}{I_z}\stackrel{•}{\theta }\stackrel{•}{\varphi }+\frac{d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_z}\\ \stackrel{••}{X}=\frac{(\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi +\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Y}=\frac{(\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi -\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Z}=-g+\frac{(\hspace{0.17em}\cos \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\end{array}$$(14) can be rearranged and simplified as follows: (21) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{1}{I_x}{u}_3\\ \stackrel{••}{\theta }=\frac{1}{I_y}{u}_2\\ \stackrel{••}{\psi }=\frac{1}{I_z}{u}_4\\ \stackrel{••}{X}=\frac{1}{m}{u}_x{u}_1\\ \stackrel{••}{Y}=\frac{1}{m}{u}_y{u}_1\\ \stackrel{••}{Z}=\frac{1}{m}{u}_1-g\end{array}$$(21) with: $$\{\begin{array}{c}{u}_x=\theta \\ u_y=-\phi \end{array}$$ and (22) $$\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]=\left[\begin{array}{cccc}b& b& b& b\\ -bl& -bl& bl& bl\\ -bl& bl& bl& -bl\\ d& -d& d& -d\end{array}\right]\cdot \left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]$$(22)

Figure 3 shows the control scheme consisting of two loops: the inner loop for altitude and attitude control and the outer loop for controlling horizontal positions (X, Y). For a control system to function properly, the inner control loop must respond faster than the outer loop.

Figure 3. Four-rotor control scheme.

The reverse model block shown in Figure 3 is designed to adapt the control signals U_i to the voltage V_i of the robot’s motors (Figure 1). It contains the inverse motion matrix which calculates the square of propeller speeds Ω_i (Eq. (23)(23) $$\left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]={\left[\begin{array}{cccc}b& b& b& b\\ -bl& -bl& bl& bl\\ -bl& bl& bl& -bl\\ d& -d& d& -d\end{array}\right]}^{-1}\cdot \left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]$$(23) ) and the inverse linear models of the rotors (Eq. (24)(24) $$V_i=\left(\left(-A*{\Omega }_i\right)-C\right)/B$$(24) ). (23) $$\left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]={\left[\begin{array}{cccc}b& b& b& b\\ -bl& -bl& bl& bl\\ -bl& bl& bl& -bl\\ d& -d& d& -d\end{array}\right]}^{-1}\cdot \left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]$$(23)

Note that the motion matrix (Eq. (22)(22) $$\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]=\left[\begin{array}{cccc}b& b& b& b\\ -bl& -bl& bl& bl\\ -bl& bl& bl& -bl\\ d& -d& d& -d\end{array}\right]\cdot \left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]$$(22) ) has a nonzero determinant.

The inverse linear model of the rotor “i” is expressed by: (24) $$V_i=\left(\left(-A*{\Omega }_i\right)-C\right)/B$$(24) with: $\{\begin{array}{l}A=-\frac{k_m^2}{J{}_{\mathit{rotor}}R}\eta r^2-\frac{2d}{J_{\mathit{rotor}}}{\Omega }_H\\ B=\frac{k_m}{J_{\mathit{rotor}}R}\eta r\\ C=\frac{d}{J_{\mathit{rotor}}}{\Omega }_H^2\end{array}$

Note that the linearization of the rotor model is according to the hover speed Ω_H.

The block diagram of the closed loop system shown in Figure 3 is simulated and used to test the performance of the controllers. Note that the quadrotor dynamic model is simulated using Eq. (14)(14) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{I_r\stackrel{•}{\theta }({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)}{I_x}+\frac{(I_y-I_z)}{I_x}\stackrel{•}{\theta }\stackrel{•}{\psi }\\ +\frac{b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_x}\\ \stackrel{••}{\theta }=\frac{I_r\stackrel{•}{\varphi }(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)}{I_y}+\frac{(I_z-I_x)}{I_y}\stackrel{•}{\varphi }\stackrel{•}{\psi }+\\ \frac{b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)}{I_y}\\ \stackrel{••}{\psi }=\frac{(I_x-I_y)}{I_z}\stackrel{•}{\theta }\stackrel{•}{\varphi }+\frac{d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_z}\\ \stackrel{••}{X}=\frac{(\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi +\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Y}=\frac{(\hspace{0.17em}\sin \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi -\hspace{0.17em}\cos \hspace{0.17em}\psi \hspace{0.17em}\sin \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\\ \stackrel{••}{Z}=-g+\frac{(\hspace{0.17em}\cos \hspace{0.17em}\theta \hspace{0.17em}\cos \hspace{0.17em}\varphi )}{m}{\sum }_i\stackrel{}{F_i}\end{array}$$(14) and parameters given in Table 1. The reverse models are given in (23)(23) $$\left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]={\left[\begin{array}{cccc}b& b& b& b\\ -bl& -bl& bl& bl\\ -bl& bl& bl& -bl\\ d& -d& d& -d\end{array}\right]}^{-1}\cdot \left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]$$(23) and (24)(24) $$V_i=\left(\left(-A*{\Omega }_i\right)-C\right)/B$$(24) .

Recently, many tools are available for robot simulation including quadrotors. For example, we can find GAZEBO robot simulator and ARGoS large-scale robot simulator discussed respectively in (Howard & Koenig, 2022) and (Pinciroli, 2022). MATLAB adds also many dedicated libraries toolboxes for simulating robots and autonomous systems (Mathworks: MATLAB & Simulink for Robotics & Autonomous Systems, 2018). But in this study, we use the MATLAB/Simulink tools to simulate our model of the quadrotor represented by the set of differential equations (10)(10) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{I_r\stackrel{•}{\theta }({\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4)}{I_x}+\frac{(I_y-I_z)}{I_x}\stackrel{•}{\theta }\stackrel{•}{\psi }\\ +\frac{b\cdot l(-{\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_x}\\ \stackrel{••}{\theta }=\frac{I_r\stackrel{•}{\varphi }(-{\Omega }_1+{\Omega }_2-{\Omega }_3+{\Omega }_4)}{I_y}+\frac{(I_z-I_x)}{I_y}\stackrel{•}{\varphi }\stackrel{•}{\psi }+\\ \frac{b\cdot l(-{\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2)}{I_y}\\ \stackrel{••}{\psi }=\frac{(I_x-I_y)}{I_z}\stackrel{•}{\theta }\stackrel{•}{\varphi }+\frac{d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2)}{I_z}\end{array}$$(10) instead of the quadrotors’ black boxes provided by the tools mentioned previously.

4.1. Attitude tracking control

From the above differential equation system (21)(21) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{1}{I_x}{u}_3\\ \stackrel{••}{\theta }=\frac{1}{I_y}{u}_2\\ \stackrel{••}{\psi }=\frac{1}{I_z}{u}_4\\ \stackrel{••}{X}=\frac{1}{m}{u}_x{u}_1\\ \stackrel{••}{Y}=\frac{1}{m}{u}_y{u}_1\\ \stackrel{••}{Z}=\frac{1}{m}{u}_1-g\end{array}$$(21) , the system of transfer functions that describe the dynamics of each angle are obtained as follows: (25) $$\{\begin{array}{c}\frac{\varphi }{u_3}=\frac{1}{I_x{p}^2}\\ \frac{\theta }{u_2}=\frac{1}{I_y{p}^2}\\ \frac{\varphi }{u_4}=\frac{1}{I_z{p}^2}\end{array}$$(25)

As can be noticed, all the above transfer functions are of a double integrator type, owing to that only a fractional-order PD^μ controller will be used for attitude tracking control.

The transfer function of the fractional PD^μ controller is expressed as $C\left(p\right)=K_p+K_d{p}^{\mu }.$

To determine the parameters of the PD, we require a set of specifications that must depend on the dynamics of the system, especially the response time of the internal loops which must be fast to that of the external loops. For our case we take the following time response performance requirements:

1. An overshoot of 30% corresponding to a damping ratio of ξ = 0.35.

2. Response time at 63% of the steady-state value of 0.666 seconds for angles φ and θ and 1 second for angle ψ.

References (Fergani & Charef, 2016) provide expressions for setting the model reference parameters that achieve the desired time response specifications.

To achieve the above time response requirements, the parameters m and ω_c of the model reference must be set as:

1. The pulsation of unity gain in open-loop ω_u = 1.5 rad/sec for φ and θ and ω_u=1rad/sec for ψ

2. The order of the system m = 1.5.

Hence, the reference fractional-order models that meet the above specifications are: $$G_d\left(p\right)=\frac{1}{1+{\left(\frac{p}{1.5}\right)}^{1.5}} \text{for φ and θ angles}.$$ $$G_d\left(p\right)=\frac{1}{1+p^{1.5}} \text{for ψ angle}.$$

Figure 4 shows the pitch, roll, and yaw orientation set points tracking closed-loop responses of the quadrotor controlled under fractional PD^μ controllers. The initial setpoints for [φ θ ψ] are respectively set to Idrissi, Salami, et al. (2021), Idrissi and Annaz (2020), and Roy et al. (2021), then changed at t = 70s to Behrooz et al. (2018), Falanga et al. (2019), and Idrissi, Annaz, et al. (2021) and returned to the [0, 0, 0] at t = 150s. It is worth mentioning that the angles are changed within the constraints of $\pm \frac{\pi }{18}$rad.

Figure 4. Setpoint tracking responses of angles [φ θ ψ] using FOPD.

As shown, the angles follow exactly the set point without any offset at steady state. The transient responses have an overshoot of 30% and settling times of 0.6701s for both φ and θ and 1.013 s for ψ which are consistent with design specifications. These results demonstrate that the FOPD controller is able to give satisfactory control and meet desired time attitude tracking control requirements.

Figure 5 compares the closed-loop of angles under FOPD and the PD for set-point tracking and disturbance rejection. The FOPD and PD are tuned to achieve the same desired specifications. The step setpoints change for [φ θ ψ] are [5 8 10], respectively. The external step disturbances changes of 0.08Nm, −0.08Nm, and 0.08 Nm are simultaneously added at t = 25 sec in each control loop of θ, φ, and ψ respectively.

Figure 5. Comparison of FOPD and conventional PD responses for angles set-point tracking and disturbance rejection.

The simulation results in Figure 5 indicate that FOPD is able to satisfy the desired design requirements while the conventional PD controller is unable to satisfy them.

For disturbance changes, FOPDs controllers can bring the angles to their set points within 3 s while under PD it takes a longer time of about 4 sec. Hence, FOPD provides faster disturbance rejection than conventional PD controllers. The maximum errors resulting from these disturbances are 0.24° for a fractional PD^μ and 0.28° for a standard PD.

4.2. Position control

From the above differential equation system (21)(21) $$\{\begin{array}{l}\stackrel{••}{\varphi }=\frac{1}{I_x}{u}_3\\ \stackrel{••}{\theta }=\frac{1}{I_y}{u}_2\\ \stackrel{••}{\psi }=\frac{1}{I_z}{u}_4\\ \stackrel{••}{X}=\frac{1}{m}{u}_x{u}_1\\ \stackrel{••}{Y}=\frac{1}{m}{u}_y{u}_1\\ \stackrel{••}{Z}=\frac{1}{m}{u}_1-g\end{array}$$(21) , the system of transfer functions that describe the dynamics of horizontal positions are obtained as follows: (26) $$\{\begin{array}{l}\frac{Y}{\varphi }=\frac{-g}{p^2}\\ \frac{X}{\theta }=\frac{g}{p^2}\end{array}$$(26)

Each of the transfer functions can be rewritten in the following form: $G_2(p)=\frac{E_2}{u_2}=\frac{K}{p^2}.$

We obtain $E_2^{*}=\{X^{*},Y^{*}\},$ $e_2=\{e_X,e_Y\}$ $u_2=\{{\theta }^{*},{\varphi }^{*}\}$ $K\{g,-g\}$ and $E_2=\{X,Y\}.$

The closed-loop system can be represented in a block diagram as depicted in Figure 6. As mentioned above, the dynamic model used is valid only for small angles variations. To respect these constraints, a saturation block should be placed in a cascade with the angular position controller as shown in Figure 6. The outputs of the controller are angles of pitch and roll.

Figure 6. General representation of closed-loop position control.

Since the model is valid only for small-angle variations, the saturation block limits u₂ to be between $\pm \frac{\pi }{18}$ rad that corresponds to ±10°.

The desired performance requirements for this application are as follows:

1. An overshoot of 40% corresponding to a damping ratio of ξ = 0.274.

2. Time response at 63% of the steady-state value of 3.2 seconds. It is worth mentioning that this requirement is more important to achieve compared to that of the angles.

To achieve the above time response requirements, the parameters m and ω_c of the model reference must be set as (1)(1) $$\{\begin{array}{l}{\Gamma }_i=\frac{d}{dt}\left(\frac{\partial L}{\partial q_i}\right)-\frac{\partial L}{\partial q_i}\\ L=T-V\end{array}$$(1)  The gain crossover frequency in open-loop ω_u = 0.333 rad/sec.

1. The order of the system m = 1.5.

Hence, the reference fractional-order model that meets the above specifications is: (27) $$G_d\left(p\right)=\frac{1}{1+{\left(\frac{p}{0.333}\right)}^{1.6}}\text{ for both X and Y positions}$$(27)

Figure 7 compares the set-point responses of horizontal positions control of a quadrotor under FOPD along with that of the reference model. The setpoints used in this simulation for [X, Y] are initially (5, 10) then changed to (2, 3) at t = 70s and returned to the origin [0, 0] at t = 150s.

Figure 7. Setpoint responses of positions X and Y using Fractional-order PDµ.

As indicated in Figure 7, the horizontal position setpoint responses of the quadrotor under the FOPD controller are almost similar to that of the desired reference model. This demonstrates that FOPD is able to achieve satisfactorily the assigned design specifications (i.e. 40.2% overshoot and response time at 63% of the steady value, of 3.29 s).

Figure 8 shows the comparison of fractional PD and the conventional PD responses for setpoint change and disturbance rejections. Recall that both PD and FOPD are designed using the same desired requirements.

Figure 8. Comparison between FOPD and conventional PD for the control of positions X and Y of quadrotor for set-point tracking and disturbance rejection responses.

The simulation results illustrated in Figure 8, demonstrate that the fractional PD^μ controller is able to meet accurately the desired control requirements while the conventional PD controller is unable to achieve them accurately.

Figure 8 illustrates also the responses of FOPD and PD to external disturbances. The external disturbances applied at t = 50s are step changes of −0.2 Nm for the Y-axis control loop and 0.2Nm for the X-axis control loop.

Figure 8 shows that both controllers are able to reject disturbance but the FOPD rejected it faster in about 10 s while conventional PD in 13 s. However, it can be noticed that the positions responses under both controllers result in offsets after disturbance changes. The maximum offsets resulting from these disturbances are 0.225 m for a FOPD and 0.242 m for a PD.

Similarly, for the altitude control Z, the derived transfer function that describes the altitude dynamic is as follows:$\frac{Z}{r}=\frac{1}{p^2}$ with $r=\frac{1}{m}{u}_1-g.$

The desired performance requirements are as follows:

1. An overshoot of 40% corresponding to damping ratio of ξ = 0.274.

2. Time response at 63% of the steady-state value of 2.5 seconds. It is worth mentioning that this requirement is more important to achieve compared to that of the angles.

To achieve the above time response requirements, the parameters m and wc of the model reference must be set as:

1. The gain crossover frequency in open-loop ω_u = 0.42 rad/sec

2. the order of the system m = 1.6.

Hence, the reference fractional-order model that meets the above specifications is as follows: (28) $$\mathrm{ }{G}_d\left(p\right)=\frac{1}{1+{\left(\frac{p}{0.42}\right)}^{1.6}}$$(28)

Figure 9 compares the setpoint responses for altitude control of the quadrotor using the FOPD controller and the reference model. The setpoint is set initially to 12 m, then changed to 5 m at t = 70s and returned to the origin at t = 150s.

Figure 9. Altitude set-point responses of the fractional PDμ controller and reference model.

The simulation results show that the altitude setpoint response of the quadrotor using FOPD follows perfectly the response of the reference model. This indicates that the FOPD controller is able to achieve the desired requirements assigned. Figure 10 shows altitude closed-loop time response under PD controller for setpoint changes similar to that of Figure 9. As can be observed in Figure 10, the output exhibits an overshoot of 44.3% and a response time of 63% of the output steady state of 3.72 s. These response performances indicate that conventional PD is unable to achieve the desired design requirement. Comparing the responses obtained in Figures 9 and 10, it can be noticed that using FOPD the response presents less overshoot, less oscillation, and less settling time than under PD. This further demonstrates that FOPD performance is superior to that of a conventional PD controller.

Figure 10. Altitude set-point responses using conventional PD controller.

5. Conclusions

In this paper, the performance of FOPD type control to control a four-rotor type flying robot in six degrees of freedom is investigated. The FOPD parameters are adjusted using the reference model method represented by the ideal of Bode transfer function in a closed loop.

Extensive simulations are performed in order to test and evaluate the performance of the proposed controller for set point change and external disturbances rejections. The results obtained are very satisfactory and meet the assigned design requirements and specifications. In addition, a comparative study of this corrector with the classic PD was also made. The results obtained have shown that a fractional PD^μ corrector is more efficient than the standard PD corrector since it responds exactly to the requested specifications and that it ensures strong robustness for external disturbances in comparison with the conventional PD.

Disclosure statement

No potential conflict of interest was reported by the authors.