Nomenclature

$AFCS$

autonomous flight control system

${C_M}$

pitching moment coefficient

${C_x}$

axial force-body axis

${C_z}$

normal force-body axis

FOPID

fractional-order PID controller

GS

gain-scheduling

$h$

altitude

${J_j}\left( t \right)$

cost function for the switching law

${J_{LQR}}$

performance index for LQR control

LPV

linear parameter varying

LQR

linear quadratic regulator

$M$

moment along aircraft Y-axis

$m$

UAV mass

PID

proportional-integral-derivative controller

$p$

angular velocity along aircraft $X$ -axis

$q$

angular velocity along aircraft $Y$ -axis

$S$

UAV mass

UAV

unmanned aerial vehicle

$r$

angular velocity along aircraft Z-axis

$S$

wing area

${S_{prop}}$

area swept by the propeller

$SPSA$

simultaneous perturbation stochastic approximation

$u$

a scalar and represents linear velocity along ${X_b}$

${u_c}$

vector of control signal

${V_{air}}$

airspeed

Greek Symbol

$\alpha $

angle-of-attack

${\delta _e}$

elevator deflection angle

${\delta _T}$

throttle opening

${\rm{\Gamma }}\left( . \right)$

Gamma function

$\phi $

roll angle

$\psi $

yaw angle

$\omega $

linear velocity along aircraft $Z$ -axis

$\rho $

air density

$\theta $

pitch angle

2.0 Uav modelisation

UAVs are attracting great interest for their guidance and control, particularly in view of the control techniques [Reference Kissoum, Ladaci and Charef7]. So their configuration and modeling is important for a more efficient design. In this work we consider flying wings, which are more profitable in terms of aerodynamic coefficient.

A linearised equation of motion can be obtained based on the small-disturbance theory [Reference Zhiping and Fang30, Reference Ren and Beard31].

The main idea is to linearise the model around a steady-state condition (trimmed flight). Assuming that deviations from this equilibrium state are small, a linearised model is expected to provide useful and fairly accurate representation of the nonlinear system.

In this theory, each variable in the model is assumed to have a nominal value (at trimmed flight, indexed $0$ ), plus a disturbance value; for instance, ${u_0}$ is the nominal value for linear velocity along body ${X_b}$ axis, and ${\rm{\Delta }}u$ is the small perturbation.

The equation of motion of the plane is derived from the Newton’s Second Law of motion. The relationships between forces in the body frame $\left( {{F_X},{F_Y},{F_Z}} \right)$ , moments $\left( {L,M,N} \right)$ and aircraft linear $\left( {u,v,w} \right)$ and angular (p; q; r) velocities.

(1) \begin{align}{F_X} &= m\left( {\dot u + qw - rv} \right) \nonumber\\ {F_Y} &= m\left( {\dot v + ru - pw} \right) \nonumber \\ {F_Z} &= m\left( {\dot w + pv - qw} \right) \nonumber \\[-35pt] \end{align}

(2) \begin{align} L &= {I_x}\dot p - {I_{xz}}\dot r + qr\left( {{I_z} - {I_y}} \right) - {I_{xz}}pq \nonumber \\ M &= {I_y}\dot q + rp\left( {{I_x} - {I_z}} \right) + {I_{xz}}\left( {{p^2} - {r^2}} \right) \nonumber \\ N &= - {I_{xz}}\dot p + {I_z}\dot r + pq\left( {{I_y} - {I_x}} \right) + {I_{xz}}qr \end{align}

where $m$ is the UAV mass, $g$ is the gravitational acceleration and is $ \approx 9.8/{s^2}$ .

Force and moment components can be decomposed in three sub-components of thrust (created by aircraft engine), gravitational (due to the earth gravity) and aerodynamic (produced due to the governing rules of aerodynamics) [Reference Ergöçmen34].

When expressed in the body frame, the gravitational force is a function of aircraft orientation in space and depends on pitch ( $\theta $ ) and roll ( $\phi $ ) angles. The gravity creates no moment since it acts through the UAV center of gravity. Force components due to gravity expressed in body frame, can be computed as:

(3) \begin{align} {({F_X})_{gravity}} &= - mg{\rm{\;\;}}sin\left( \theta \right) \nonumber \\ {({F_Y})_{gravity}} &= mg{\rm{\;\;}}cos\left( \theta \right){\rm{\;}}sin\left( \phi \right) \nonumber \\ {({F_Z})_{gravity}} &= mg{\rm{\;\;}}cos\left( \theta \right){\rm{\;}}cos\left( \phi \right) \end{align}

where $\theta $ and $\phi $ are aircraft pitch and roll angles, respectively. Assume that the sum of aerodynamic and thrust forces are $X$ , $Y$ and $Z$ . Force components due to gravity expressed in body frame, can be computed as:

(4) \begin{align} X - mgsin\left( \theta \right) &= m\left( {\dot u + qw - rv} \right) \nonumber \\ Y + mgcos\left( \theta \right)sin\left( \phi \right) &= m\left( {\dot v + ru - pw} \right) \nonumber \\ Z + mgcos\left( \theta \right)cos\left( \phi \right) &= m\left( {\dot w + pv - qw} \right)\end{align}

Equation (4) represent the dynamics and kinematics of the UAV.

If we consider $X$ (the force along the body ${X_b}$ axis) as described in Equation (4). By introducing the small-disturbance notation into this equation and simplifying the result, following equation will be obtained as,

(5) \begin{align}{\rm{\Delta }}X - mg{\rm{\;\Delta }}\theta {\rm{\;\;}}cos{\rm{\;}}{\theta _0} = m{\rm{\;\Delta }}\dot u + m{\rm{\;}}{w_0}{\rm{\;\;\Delta }}q\end{align}

Where ${\rm{\Delta }}X$ is the change in the force along the ${X_b}$ direction (contribution of aerodynamic force and thrust). It is possible to express ${\rm{\Delta }}X$ in terms of perturbation variables using Taylor series [Reference Nelson32].

Considering the longitudinal axis, the $X$ force equation can be rewritten as

(6) \begin{align}\dot u = rv - qw - gsin\theta + \frac{X}{m}\end{align}

where $X$ represents the combination of aerodynamic and propulsive forces (i.e. $X = {F_{{X_{aero}}}} + {F_{{X_{prop}}}}$ ). Assuming that Ref. (Reference Beard and McLain33) $\phi = p = r = \beta = v = 0$ , ${F_{{X_{aero}}}}$ and ${F_{{X_{prop}}}}$ can be computed as:

(7) \begin{align}{F_{{X_{aero}}}} &= {{\rho \left( {{u^2} + {w^2}} \right)S} \over m}\left[ {{C_{{x_0}}} + {C_{{x_\alpha }}}\alpha + {C_{{x_{{\delta _e}}}}}{\delta _e}} ) \right] + {{\rho \sqrt {{u^2} + {w^2}} S} \over m}{C_{{X_q}}} \hat{c}q \nonumber \\ {F_{{X_{prop}}}} &= \frac{1}{2}\rho {S_{prop}}(({K_{motor}}{\delta _T}{)^2} - V_{air}^2)\end{align}

where $\rho $ is the air density, $u$ and $w$ are velocities along body $x$ -axis and $z$ -axis, respectively, ${C_{{X_0}}}$ is the axial force coefficient, ${C_{{X_\alpha }}}$ is the derivative of the axial force with respect to $\alpha $ , ${C_{{x_{{\delta _e}}}}}$ is the derivative of the axial force with respect to ${\delta _e}$ and $\hat c$ is the wing chord.

Substituting back Equation (7) into (6) with the assumption that $\phi = p = r = \beta = v = 0$ , we will get:

(8) \begin{align} \dot u = - qw - gsin\theta &+ {{\rho \left( {{u^2} + {w^2}} \right)S} \over m}\left[ {{C_{{x_0}}} + {C_{{x_\alpha }}}\alpha + {C_{{x_{{\delta _e}}}}}{\delta _e}} )\right] + {{\rho \sqrt {{u^2} + {w^2}} S} \over m}{C_{{X_q}}}\hat{c}q + \nonumber \\ & \quad + \frac{1}{2}\rho {S_{prop}}(({K_{motor}}{\delta _T}{)^2} - V_{air}^2)\end{align}

This procedure should be repeated for all non-linear forces and moments. The results (longitudinal terms) would be as follows (see Fig. 1) [Reference Ergöçmen34]:

(9) \begin{align}\left\{ {\begin{array}{l}{\dot u = - qw - g{\rm{sin}}\theta + \dfrac{{\rho \left( {{u^2} + {w^2}} \right)s}}{{2m}}\left[ {{C_{{x_0}}} + {C_{{x_\alpha }}}\alpha + {C_{{x_{{\delta _e}}}}}{\delta _e}} \right] + \dfrac{{\rho \sqrt {{u^2} + {w^2}s} }}{{2m}}{C_{{x_q}}}\hat {c}q } \\[7pt] { \qquad + \dfrac{1}{{2m}}\rho {s_{prop}}(({k_{motor}}{\delta _T}){^2}) - V_\alpha ^2)} \\[7pt]{\dot w = - qu - g{\rm{cos}}\theta {\rm{cos}}\phi + \dfrac{{\rho \left( {{u^2} + {w^2}} \right)s}}{{2m}}\left[ {{C_{{z_0}}} + {C_{{z_\alpha }}}\alpha + {C_{{z_{{\delta _e}}}}}{\delta _e}} \right] + \dfrac{{\rho \sqrt {{u^2} + {w^2}s} }}{{2m}}{C_{{z_q}}}\bar cq} \\[7pt]{\dot q = \dfrac{1}{{2{I_y}}}\rho \left( {{u^2} + {w^2}} \right)\bar cs\left[ {{C_{{m_0}}} + {C_{{m_\alpha }}}\alpha + {C_{{m_{{\delta _e}}}}}{\delta _e}} \right] + \dfrac{1}{{2{I_y}}}\rho \sqrt {{u^2} + {w^2}} \bar cs{C_{{m_q}}}\bar cq} \\[7pt]{\dot \theta = q{\rm{cos}}\phi - r{\rm{sin}}\phi = q}\end{array}} \right. \nonumber \\[-15pt]\end{align}

Figure 1. UAV forces in stability and body axes.

In order to obtain an operational model, we will use the so-called small-disturbance method. It means that the variables are replaced by their nominal values plus a small disturbance value, then by taking a linear approximation we obtain the following equations using the Jacobians:

(10) \begin{align}\dfrac{{df}}{{dx}} = \left[ {\begin{array}{*{20}{l}}{\dfrac{{\partial \dot u}}{{\partial u}}} \quad {}{\dfrac{{\partial \dot u}}{{\partial w}}} \quad {}{\dfrac{{\partial \dot u}}{{\partial q}}} \quad {}{\dfrac{{\partial \dot u}}{{\partial \theta }}}\\{} {}{} {}{} {}{}\\{\dfrac{{\partial \dot w}}{{\partial u}}} \quad {}{\dfrac{{\partial \dot w}}{{\partial w}}} \quad {} {\dfrac{{\partial \dot w}}{{\partial q}}} \quad {}{\dfrac{{\partial \dot w}}{{\partial \theta }}}\\{} {}{} {}{} {}{}\\{\dfrac{{\partial \dot q}}{{\partial u}}} \quad {}{\dfrac{{\partial \dot q}}{{\partial w}}} \quad {}{\dfrac{{\partial \dot q}}{{\partial q}}} \quad {}{\dfrac{{\partial \dot q}}{{\partial \theta }}}\\{} {}{} {}{} {}{}\\{\dfrac{{\partial \dot \theta }}{{\partial u}}} {} \quad {\dfrac{{\partial \dot \theta }}{{\partial w}}} \quad {}{\dfrac{{\partial \dot \theta }}{{\partial q}}} {} \quad {\dfrac{{\partial \dot \theta }}{{\partial \theta }}} \end{array}} \right]\end{align}

(11) \begin{align}\dfrac{{df}}{{d{u_c}}} = \left[ {\begin{array}{*{20}{l}}{\dfrac{{\partial \dot u}}{{\partial {\delta _e}}}} \quad {}{\dfrac{{\partial \dot u}}{{\partial {\delta _T}}}}\\{} {}{}\\{\dfrac{{\partial \dot w}}{{\partial {\delta _e}}}} {} \quad {\dfrac{{\partial \dot w}}{{\partial {\delta _T}}}}\\{} {}{}\\{\dfrac{{\partial \dot q}}{{\partial {\delta _e}}}} {} \quad {\dfrac{{\partial \dot q}}{{\partial {\delta _T}}}}\\{} {}{}\\{\dfrac{{\partial \dot \theta }}{{\partial {\delta _e}}}} \quad {}{\dfrac{{\partial \dot \theta }}{{\partial {\delta _T}}}} \end{array}} \right] \\[8pt] \nonumber \end{align}

Note that ${u_c}$ is a vector of control signal ${u_c} = {[{\delta _e}{\rm{\;\;}}{\delta _T}]^{\rm{T}}}$ , where ${\delta _e}$ is the elevator deflection angle and ${\delta _T}$ is the throttle opening, such that $u$ is a scalar and represents linear velocity along ${X_b}$ , $w$ is the linear velocity along aircraft $Z$ -axis, $\theta $ is the pitch angle and $q$ is the angular velocity along aircraft $Y$ -axis.

The UAV system outputs are the airspeed ${V_{air}}$ and the altitude $h$ .

However, controlled outputs are airspeed ${\rm{\Delta }}{V_{air}}$ and aircraft’s height ${\rm{\Delta }}h$ .

Finally, longitudinal linearised model will be given by:

(12) \begin{align}\begin{array}{*{20}{l}}{\Delta \, \dot u = {X_u} \, \Delta \, u + {X_w} \, \Delta \, w + {X_q} \, \Delta \, q - g \, {\rm{cos}} \, \theta \, \Delta \, \theta + {X_{{\delta _e}}} \, \Delta \, {\delta _e} + {X_{{\delta _T}}} \, \Delta \, {\delta _{}} + {w_0} \, \Delta \, q} {}{} {}{} {}{}\\[3pt]{\Delta \, \dot w = {Z_u} \, \Delta \, u + {Z_w} \, \Delta \, w + {Z_q} \, \Delta \, q - g \, {\rm{sin}}\theta \, \Delta \, \theta + {Z_{{\delta _e}}} \, \Delta \, {\delta _e} + {w_0} \, \Delta \, q} {}{} {}{} {}{}\\[3pt]{\Delta \, \dot q = {M_u} \, \Delta \, u + {M_w}\, \Delta \, w + {X_q} \, \Delta \, q + {M_{{\delta _e}}} \, \Delta \, {\delta _e}} {}{} {}{} {}{}\\[3pt]{\Delta \, \dot \theta = \Delta \, q} {}{} {}{} {}{}\end{array}\end{align}

2.1 State-space representation of the linearised model

Linearised equations are based on simple linear ordinary differential equations. The parameters are given in Table 1, and the numerical values for longitudinal stability derivatives are presented in Table 2 [Reference Fahlstrom and Gleason35].

With $\rho $ , ${C_M},{C_x},{C_Z},M,m,p,S + $ and ${S_{prop}}$ are air density, pitching moment coefficient, axial force-body axis, normal force-body axis, moment along aircraft Y-axis, UAV mass, angular velocity along aircraft X-axis, wing area and area swept by the propeller, respectively. We are now able to transform them into a state space model with a number of ordinary first-order differential equations [Reference Ahmad, Kumar, Bhandari and Patil9].

Table 1. Longitudinal stability derivatives

The description of the UAV longitudinal movement is described by a set of three equations:

(13) \begin{align}\left[ {\begin{array}{*{20}{l}}{\Delta \, \dot u}\\{}\\\Delta \\{}\\{\Delta\, \dot q}\\{}\\{\Delta \, \dot \theta }\\ \end{array}} \right] = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{X_u}} & {} {{X_w}} & {} {{X_q} + {w_0}} & {}{ - g \, {\rm{cos}} \, {\theta _0}}\\{} {}{} {}{} {}{}\\{{Z_u}} & {}{{Z_w}} & {}{{Z_q} - {w_0}} & {}{ - g \, {\rm{sin}} \, {\theta _0}}\\{} {}{} {}{} {}{}\\{{M_u}} & {}{{M_w}} & {}{{M_q} - {w_0}} & {}0\\{} {}{} {}{} {}{}\\0 {}& 0 & 1 & {} 0\\ \end{array}} \right]\left[ {\begin{array}{c}{\Delta \, u}\\{}\\{\Delta \, w}\\{}\\{\Delta \, q}\\{}\\{\Delta \, \theta }\\ \end{array}} \right] + \left[ {\begin{array}{c@{\quad}c} {{X_{{\delta _{}}}}} & {}{{X_{{\delta _e}}}}\\{} {}{} {}{} {}{}\\{{Z_{{\delta _e}}}} & {}0\\{} {}{} {}{} {}{}\\{{M_{{\delta _e}}}} & {}0\\{} {}{} {}{} {}{}\\0 & {}0 \\ \end{array}} \right]\left[ {\begin{array}{c} {{\Delta _{{\delta _e}}}}\\{} {}{}\\{{\Delta _{{\delta _T}}}}\\ \end{array}} \right]\end{align}

Table 2. Numerical values for longitudinal stability derivatives

Table 3 tabulates trim condition values. It is possible to compute small variations in each force and moment element $\left( {\Delta X,\Delta Y,\Delta Z,\Delta L,\Delta M,\Delta N} \right)$ as an effect of variation in a particular variable $\left( {u,w,\delta , \ldots } \right)$ [Reference Botez36]:

Table 3. Trim condition values

In this equation, the zero index denotes the trimmed light condition ( ${\theta _0}$ is the pitch angle for trim condition, ${u_0}$ is the linear velocity along body ${X_b}$ axis etc.).

2.2 LPV model based on the angle-of-attack ${\bf{\alpha }}$

In order to obtain the LPV model set, we consider the scheduling variable here as the angle-of-attack $\alpha $ as mentioned in Table 1. By taking five different values of the parameter $\alpha $ (which can be considered as the centres of five intervals of variation for the angle-of-attack), we obtain five state-space models reported in Table 4.

Table 4. State-space models vs angle-of-attack ${\rm{\alpha }}$

3.0 Fractional-order pid controller

Since the early days of fractional calculus theory, many have seen the value of this new tool in improving and positively affecting all aspects of modern engineering. The mathematical definitions of fractional order differentiation are numerous and diverse, the best known and used of which are those of Grünwald-Letnikov, Riemann-Liouville and Caputo [Reference Podlubny37].

3.1 Definition of fractional integration

Riemann-Liouville definition for non-integer integral is the following:

(14) \begin{align}D_a^{ - \nu }f\left( t \right) = {1 \over {{\rm{\Gamma }}\left( \nu \right)}}\mathop \smallint \nolimits_a^t {(t - \zeta )^{\nu - 1}}f\left( \zeta \right)d\left( \zeta \right)\end{align}

with the Gamma function presented as

(15) \begin{align}{\rm{\Gamma }}\left( x \right) = \mathop \smallint \nolimits_0^\infty {y^{x - 1}}{e^{ - y}}dy,\end{align}

such that $\left( {a,t} \right) \in {\Re ^2}$ with $a \lt t$ and ( $0 \lt \nu \lt 1$ , $\nu \in \Re $ ).

3.2 Definition of fractional derivation

The fractional-order derivative of order $\mu \in R$ in the sense of Riemann–Liouville is given as:

(16) \begin{align}{{\rm{\;}}_{RL}}D_{{t_0}}^\mu f\left( t \right) = {1 \over {{\rm{\Gamma }}\left( {n - \mu } \right)}}{{{d^n}} \over {d{t^n}}}\mathop \smallint \nolimits_{{t_0}}^t {(t - \tau )^{n - \mu - 1}}f\left( \tau \right)d\tau \end{align}

where ${\rm{\Gamma }}\left( . \right)$ is Euler’s gamma function and $n - 1 \lt \mu \lt n$ , $n \in {R^ + }$ .

It is well known that industrial controlled processes are habitually discretised, that justifies the approximation of Equation (16). The literature proposes many numerical tools for fractional-order continuous and discrete differential equations [Reference Hajiyev, Soken and Vural38]. In our work, we are concerning only with integrator and derivative approximation.

3.3 Singularity functions

In this work, and in order to implement fractional-order models, we will use the so-called singularity function method developed by Charef et al. [Reference Charef, Sun, Tsao and Onaral39].

3.3.1 Fractional-order integrator approximation

The transmittance representing the fractional-order integrator is approximated in the frequency domain by an irrational transfer function as:

(17) \begin{align}{G_I}\left( s \right) = \frac{1}{{{s^\mu }}}\end{align}

such that $0 \lt \mu \lt 1$ .

In an operational frequency band $\left[ {{\omega _b},{\omega _h}} \right]$ , it can be modelled by a fractional power pole (FPP) whose transmittance is:

(18) \begin{align}{G_I}\left( s \right) = \frac{{{K_I}}}{{{{\left( {1 + \frac{s}{{{\omega _c}}}} \right)}^\mu }}}\end{align}

Suppose that for $\omega \in \left[ {{\omega _b},{\omega _h}} \right]$ ,with $\omega \gt \gt {\omega _c}$ therefore

(19) \begin{align}{G_I}\left( s \right) = \frac{{{K_I}}}{{{{\left( {\frac{s}{{{\omega _c}}}} \right)}^\mu }}} = \frac{{{K_I}\omega _c^\mu }}{{{s^\mu }}} = \frac{1}{{{s^\mu }}}\end{align}

were ${K_I} = \frac{1}{{\omega _c^\mu }}$ and ${\omega _c}$ is the $ - 3mdB$ frequency corner of the FPP, which is computed based on the low frequency ${\omega _l}$ , as ${\omega _c} = \sqrt {{{10}^{\frac{\varepsilon }{{10\alpha }}}} - 1} $ were $\varepsilon $ is the highest allowed gap between the slopes of the fractional-order integrator of Equation (17) and the FPP of Equation (18) in the considered frequency operational band $\left( {{\omega _l},{\omega _h}} \right)$ .

The FPP of Equation (18) and equally the integrator have to be modeled by means of a linear time-invariant system and for that aim we have to approximate its transmittance by a rational one. This method implies the approximation the $ - 20\mu {\rm{\;\;}}dB/dec$ slope on the Bode plot of the FPP by a set of alternative slopes of $ - 20{\rm{\;\;}}dB/dec$ and $0{\rm{\;\;}}dB/dec$ equivalent to alternative poles and zeros on the negative part of the real axis of the complex plane such that ${p_0} \lt {z_0} \lt {p_1} \lt {z_1} \lt \ldots \lt {z_{N - 1}} \lt {p_N}$ ; Hence, the approximation is given by

(20) \begin{align}{G_I}\left( s \right) = {{{K_I}} \over {{{\left( {1 + {s \over {{\omega _c}}}} \right)}^\mu }}} = {K_I}{{\mathop \prod \nolimits_{i = 0}^{N - 1} \left( {1 + {s \over {{z_i}}}} \right)} \over {\mathop \prod \nolimits_{i = 0}^N \left( {1 + {s \over {{p_i}}}} \right)}}\end{align}

Now by making use of a simple graphical technique that begins with the error specification $y$ in decibels and the frequency band ${\omega _{max}}$ , which can be $100{\omega _h}$ , the parameters $a,b,{p_0},{z_0}$ can be computed as [Reference Ladaci and Bensafia40]

\begin{align*}a = {10^{\frac{y}{{10\left( {1 - \mu } \right)}}}},b = {10^{\frac{y}{{10\mu }}}}\end{align*}

with:

\begin{align*}{p_0} = {\omega _c}\sqrt b ,{\rm{\;\;}}{z_0} = a{p_0}\end{align*}

The poles ${p_i}$ and the zeros ${z_i}$ of Equation (20) are in a form of geometric progression. We can then obtain them based on the above parameters as

(21) \begin{align}{p_i} &= {(ab)^i}{p_0},{\rm{\;\;}}for{\rm{\;\;}}i = 0,1, \ldots ,N \nonumber \\[4pt]{z_i} &= {(ab)^i}a{p_0},{\rm{\;\;}}for{\rm{\;\;}}i = 0,1, \ldots ,N - 1\end{align}

with

(22) \begin{align}N = {\rm{Integer}}\left[ {\frac{{log\left( {\frac{{{\omega _{max}}}}{{{p_0}}}} \right)}}{{log\left( {ab} \right)}} + 1} \right] + 1\end{align}

It results that the fractional-order integrator is equivalent to a rational transfer function expressed as

(23) \begin{align}{G_I}\left( s \right) = \frac{1}{{{s^\mu }}} = \frac{{{K_I}}}{{{{\left( {1 + \left( {\frac{s}{{{\omega _c}}}} \right)} \right)}^\mu }}}\end{align}

So ${G_I}\left( s \right)$ can be approximated as follow

(24) \begin{align}{G_I}\left( s \right) = {K_I}{{\mathop \prod \nolimits_{i = 0}^{N - 1} \left( {1 + {s \over {{{(ab)}^i}a{p_0}}}} \right)} \over {\mathop \prod \nolimits_{i = 0}^N \left( {1 + {s \over {{{(ab)}^i}{p_0}}}} \right)}} = \mathop \sum \limits_{i = 0}^N {{{h_i}} \over {\left( {1 + {s \over {{{(ab)}^i}{p_0}}}} \right)}}\end{align}

3.3.2 Fractional-order derivative’s approximation

The fractional-order derivative equivalent Laplace-domain transmittance is expressed by the following irrational transmittance:

(25) \begin{align}{G_D}\left( s \right) = {s^\mu }\end{align}

where $s = j\omega $ is the complex frequency and $\mu $ is a positive real number such that $0 \lt \mu \lt 1$ .

In the considered operational frequency band of $\left[ {{\omega _b},{\omega _h}} \right]$ , this fractional-order operator can be represented by a fractional power zero (FPZ) whose transmittance is the following,

(26) \begin{align}G\left( s \right) = {K_D}{\left( {1 + \frac{s}{{{\omega _c}}}} \right)^\mu }\end{align}

Suppose that for $\omega \in \left[ {{\omega _b},{\omega _h}} \right]$ , with $\omega \gt \gt {\omega _c}$ therefore

(27) \begin{align}{G_D}\left( s \right) = {K_D}{\left( {1 + \frac{s}{{{\omega _c}}}} \right)^\mu }\end{align}

Were ${K_D} = \omega _c^\mu $ .

Following the preceding approach, the approximation of the FPP [Reference Charef, Sun, Tsao and Onaral39], the fractional-order derivative can be approximated by a rational transmittance represented in the operational frequency band by

(28) \begin{align}{G_D}\left( s \right) = {s^\mu } = {K_D}{\left( {1 + {s \over {{\omega _c}}}} \right)^\mu } \approx {K_D}{{\mathop \prod \nolimits_{i = 0}^N \left( {1 + {s \over {{{(ab)}^i}{z_0}}}} \right)} \over {\mathop \prod \nolimits_{i = 0}^N \left( {1 + {s \over {{{(ab)}^i}a{z_0}}}} \right)}}\end{align}

with

\begin{align*}a = {10^{\frac{y}{{10\left( {1 - \mu } \right)}}}},b = {10^{\frac{y}{{10\mu }}}}\end{align*}

and:

\begin{align*}{z_0} = {\omega _c}{10^{\left( {\frac{y}{{20\mu }}} \right)}},{\rm{\;\;}}{p_0} = a{z_0}\end{align*}

The zeros ${p_i}$ and the poles ${z_i}$ of the rational function approximation are given as:

(29) \begin{align} {z_i} &= {(ab)^i}{z_0},{\rm{\;\;}}pouri = 0,1, \ldots ,N \nonumber \\[3pt] {p_i} &= {(ab)^i}a{z_0},{\rm{\;\;}}pouri = 0,1, \ldots ,N\end{align}

with

(30) \begin{align}N = {\rm{Integer}}\left[ {\frac{{log\left( {\frac{{{\omega _{max}}}}{{{p_0}}}} \right)}}{{log\left( {ab} \right)}} + 1} \right] + 1\end{align}

3.4 Fractional-order PID controller (FOPID)

Since the pioneering work of Podlubny [Reference Podlubny37] in 1999, many researchers and engineers have been attracted to the fractional order PID controller for its augmented flexibility and performance enhancement [Reference Edet and Katebi41–Reference Krishna and Rao43]. The FOPID controller used in this work is given by the following transfer function:

(31) \begin{align}{C_{PID}}\left( s \right) = {K_p}\left( {1 + {{\left( {\frac{1}{{{T_i}s}}} \right)}^\lambda } + {{\left( {{T_d}s} \right)}^\mu }} \right)\end{align}

where $\lambda $ , $\mu $ are positive real numbers.

4.0 Lqr controller

The well-known LQR is able to stabilise the feedback control closed-loop with optimal gains and offer high-performance synthesis of control systems [Reference Chacko, Neeraj and Abraham44]. In our context, it is used with the aim of controlling the UAV altitude. MATLAB/Simulink is the simulation platform allowing the design and tuning of the LQR regulator and the simulation to UAV identified model. LQR is based on the theory of optimal control whose objective is to minimise a cost function by optimising a dynamic system, precisely, minimising the deviation of the UAV altitude.

In fact, the main control objective is to achieve certain closed-loop poles that define some desired performance measures (response time, overshoot level, etc.) by placing the closed-loop system eigenvalues in a desired location using a feedback gain. Unfortunately, this statement is not so useful for practical application. The reason is that such a level of performance might require a control signal that is too large for a physical actuator to generate. Such a feedback gain would not be realisable from implementation point of view. In order to find a trade-off between performance and control signal magnitude, LQR control design is a good solution to overcome such problems [Reference Willis, Johnson and Beard45].

One of its main advantages is its easy design and increased state variables accuracy. This property guarantees the stability of the control system [Reference Wang, Ni, Zhou, Pardalos, Fang and Fei46, Reference Zhao, Yang, Li, Xu, She and Yan47].

Let us consider an LTI system, expressed as

(32) \begin{align}\dot x = Ax + Bu + G\delta \end{align}

where $x$ is the state vector, $u$ the control signal and $\delta $ a disturbance.

Let a state variable feedback regulator be,

(33) \begin{align}u = - Kx\end{align}

where $K$ is the matrix of the state feedback gain.

In order to optimise it, one has to determine the control input $u$ , which minimises the performance index ${J_{LQR}}$ . This later represents the necessary performance characteristic as well as the controller input constraint. The optimal controller has to minimise the following performance index [Reference Rui, Wang, Song, Luo, Huang, textscYu, Liu, Wu and Xiao48],

(34) \begin{align}{J_{LQR}} = \mathop \smallint \nolimits_0^\infty \left( {x{\rm{'}}Qx + u{\rm{'}}Ru} \right)dt\end{align}

where $K$ is given by;

(35) \begin{align}K = {R^{ - 1}}B{\rm{'}}P\end{align}

and $P$ is a matrix that satisfies the reduced-matrix Riccati equation,

(36) \begin{align}A{\rm{'}}P + PA - PB{R^{ - 1}}B{\rm{'}}P + Q = 0\end{align}

Then the feedback regulator,

(37) \begin{align}u = - \left( {{R^{ - 1}}B{\rm{'}}P} \right)x\end{align}

Altitude and airspeed outputs are subtracted from the desired reference altitude and airspeed to produce an error signal. Then, the control signals generated by the LQR and FOPID controllers are combined at the input of the UAV system.

These final control signals are further passed through saturation blocks, which represent limitations in aileron deflection angle ( $ - 12{\rm{\;\;}}deg \lt {\delta _e} \lt 12\;\;deg$ ) and throttle opening ( $0 \lt {\delta _T} \lt 1$ ).

5.0 Robust gain-scheduled control system design

Gain scheduling is a particularly appreciated control technique in industrial applications. This is due to the possible application of intuitive approaches and its simple implementation. However, its principal drawback is the lack of theoretical stability tools [Reference Kissoum and Ladaci49, Reference Hencey and Alleyne50].

Automation engineers often deal with systems whose model is unique and fixed or slowly varying over time. But when tackling more complex systems, models can vary in a variety of other ways: configuration changes, system faults, measurement and actuator noise, external disturbances and parameter variations.

In general, complex systems operate in multiple environments that may abruptly change from a context to another as is often the case with aircraft and UAV. The system model depends on a variable parameter, which allows the use of LPV systems theory leading to a set of linear models that become active for certain values of this parameter [Reference Leulmi, Ladaci and Schulte51].

If models are available for the different discretised parameters values, the corresponding regulators may be designed apriori. While the system is operating, we have to identify the actual parameter value in order to determine the adapted regulator.

Based on these two ideas, the proposed control strategy is to determine the best model for the current measured or estimated parameter value at every moment and use the corresponding regulator.

For the switching, the problem is to determine when the current parameters values are not satisfactory (that is, when to switch) and with which one we have to replace it (that is, switch to what). For the adjustment, the problem is to determine the law with which the parameter value has to be adjusted at every time.

In this work, linear control solutions are used to supervise the UAV’s longitudinal movements with LPV representation for a set of its possible models. LQR method with the PID controller is used to guarantee the stability of the system and a good tracking performance. This technique has been proposed in several studies for processes with invariant models [Reference He, Wang and Lee52, Reference Kim, Kim, Song and Paek53].

The main objective of this mixed control is to take benefit from the LQR optimal setting that stabilises the first control loop. Then, the fractional-order PID controller becomes easier to adjust with better performance in tracking the reference input signals. Figure 2 represents the block diagram of the proposed controller for longitudinal controller.

Figure 2. Block diagram of the proposed controller.

Figure 3 represents the proposed control method. The plant to be controlled has an input $u$ and an output $y$ . The aim of this work is to force the tracking error $e = {y_r} - y$ to converge to zero, where ${y_r}$ is the reference signal.

Figure 3. Gain-scheduled control algorithm.

5.2 Moving time switching law

In control and measurement, a switching law typically dictates which number of control, from a given list of controllers, is chosen to be used at a particular time instant, depending on particular given conditions.

In the present work, the decision when and to which controller should be the switching is based on a cost function minimisation $mi{n_j}\left( {{J_j}\left( t \right)} \right)$ for each regulator ${C_j}$ at every instant. It implies the evaluation of each model for each sample time. We propose the following cost function represented by a moving time window of length $T$ formulated as:

(38) \begin{align}{J_j}\left( t \right) = \mathop \smallint \nolimits_t^{t - T} e_j^2\left( t \right)\end{align}

At every instant, the performance factors $\left\{ {{J_j}\left( t \right)} \right\}$ are compared allowing the switching decision, as explained in the algorithm of Fig. 3.

Following the switch, the system is blocked for a ${T_{min}} \gt 0$ length waiting period, and after that, the regulator that indicates the minimum cost function is pointed out (switching target) to be applied to the plant. This waiting period is motivated by preventing the system from too fast arbitrary switching [Reference Kissoum, Ladaci and Charef7, Reference Ladaci and Khettab54].

6.0 Application for uav airspeed and altitude control

In this section we will apply the proposed gain scheduled strategy to the control of the UAV Altitude and velocity. ${\rm{Five}}$ different models are considered for modelling the LPV system. Both integer order and fractional order controllers are considered for either the altitude or velocity of the plane, and the results are given in the ideal case and in presence of random disturbances.

Consider a real LPV system given by its state space linearised model presented in Section 2.1 and given in the form of Equation (31).

We suppose that the plant parameters will change periodically between five constant values of the angle-of-attack $\alpha $ , such that: $\alpha = \left[ {6,{9^ \circ }\left( {plant} \right){\rm{\;\;}}7,{5^ \circ }{\rm{\;\;}}7,{9^ \circ }{\rm{\;\;}}13,{2^ \circ }{\rm{\;\;}}{{14}^ \circ }{\rm{\;\;}}{{20}^ \circ }}\right]$ .

Then we have five linear models representing the LPV plant whose different linear transmittances are given in Table 5.

Table 5. Different models of the LPV models of altitude and velocity of the UAV

First of all, we apply LQR feedback controller for the ${\rm{five}}$ system’s models in (using state space formula), we get stability of the system, like mentioned in Table 6.

Table 6. FOPID controllers set for UAV altitude

Then, for each LQR closed-loop transfer function, we compute a fractional-order PID controller as in Equation (30).

Table 7 gives the regulators of the UAV altitude. Wheras, Table 8 gives the regulators of the UAV velocity.

Table 7. FOPID controllers set for UAV altitude

Table 8. FOPID controllers set for UAV velocity

6.1 Control with ideal conditions (without noises)

We begin the UAV control experiment by considering that the plant is in ideal conditions. This means that we minimise external factors like wind, obstacles, additive actuators’ and sensors’ noises and signal interference, which enables more precise navigation and manoeuvers.

6.1.1. Simulations of UAV altitude

In nominal flight conditions, the designed controllers must achieve stabilisation of the UAV altitude around $h \approx 150m$ . This is the normal flying altitude.

In our experimental test we will apply gradual reference steps for the desired altitude ${h_{des}}$ to test the reference tracking and closed-loop system performance.

The simulation results of the UAV altitude $h$ control in ideal conditions (without additive noises) using the proposed robust gain-scheduling controller are given in Fig. 4.

Figure 4. UAV altitude response in ideal conditions with integer and fractional PID IDLQR controller: (a) Output signal ${\rm{h}}$ (b) Control signal (c) Switching index.

It can be observed from Fig. 4(a) that the altitude response tracks the reference signal with rapidity and precision. The fractional order case does not present any overshoot, while in the integer case, there are important peaks at the steps instants. The same remark is observable on the control signal in Fig. 4(b) which is smoother with the fractional order controller which implies a reduction of the energy costs. Besides, one can see that the switching system is more active in the fractional case in Fig. 4(c), which improves the precision of the control comparatively to the integer case where the controller stabilises with the model 1 after a short time.

6.1.2. Simulations of UAV velocity

In nominal flight conditions, the designed controllers must achieve stabilisation of airspeed (velocity) around ${V_{air}} \approx 10m/s$ . This is the normal cruise speed.

In our experimental test we will apply gradual reference steps for ${V_{ai{r_{des}}}}$ to test the reference tracking and closed-loop system performance.

The simulation results of the UAV velocity (airspeed ${V_{air}}$ ) control in ideal conditions (without additive noises) using the proposed robust gain-scheduling controller are given in Fig. 5.

Figure 5. UAV airspeed response in ideal conditions with integer and fractional PID IDLQR controller: (a) Output signal ${{\rm{V}}_{{\rm{air}}}}$ (b) Control signal (c) Switching index.

It can be noticed again from Fig. 5(a) that the velocity response tracks the reference signal with good temporal performance. The overshoot in the fractional-order response is very small and negligible comparatively to the integer case. The same comment can be stated for the control signal in Fig. 5(b) which is smoother with the fractional order controller which also implies a reduction of the energy costs. Besides, the switching system is more dynamic in the fractional-order case in Fig. 5(c), which allows better adaptation and tracking performance.

6.2 Case of UAV with sensors random noises

In this second experiment setting, the system is supposed to be subject to additive output random measure noises.

6.2.1 Simulations of UAV altitude (noisy case)

The simulation results of the UAV altitude control that is affected by random noises using the proposed robust gain-scheduling controller are given in Fig. 6.

Figure 6. UAV altitude response presence of sensors random noises with integer and fractional PID IDLQR controller: (a) Output signal ${\rm{h}}$ (b) Control signal (c) Switching index.

The altitude output signal in Fig. 6(a) shows a lower level of oscillations in the fractional-order response with the absence of picks that are remarkable in the integer-order case. This fact is still more visible for the control signals in Fig. 6(b). Same for the switching system, which is more active in the fractional-order case as can be seen in Fig. 6(c),

Figure 7. UAV airspeed response in ideal conditions with integer and fractional PID IDLQR controller: (a) Output signal ${{\rm{V}}_{{\rm{air}}}}$ (b) Control signal (c) Switching index.

6.2.2 Simulations of UAV velocity (noisy case)

The simulation results of the UAV airspeed control affected with sensor’s random noises using the proposed robust gain-scheduling controller are given in Fig. 7.

The airspeed output signal in Fig. 7(a) shows a similar level of oscillations in the fractional-order response with the absence of high picks. The same remark can be given for the control signals in Fig. 7(b). However, it is obvious that there only a few switches in the beginning of simulation for the integer order controller, whereas many controllers and models are involved with the fractional order response as could be seen in Fig. 7(c),