2.1. Body dynamics

Two coordinate systems are introduced in Fig. 1 for dynamical analysis of the quadrotor’s body: the first is the earth-fixed inertial frame defined by $\left\{ E \right\}\,:\,\{ {\vec O,( {\hat X,\hat Y,\hat Z} )} \}$ , and the second is the body-fixed frame defined by $\left\{ B \right\}\,:\,\{ {\vec P,( {\hat x,\hat y,\hat z} )} \}$ . Vector $\vec O = {\left[ {0,0,0} \right]^T}$ indicates the origin of the earth-fixed inertial frame $\left\{ E \right\}$ . Also, vector $\vec P = {\left[ {X,Y,Z} \right]^T}$ indicates the origin of the frame $\left\{ B \right\}$ in the earth-fixed inertial frame $\left\{ E \right\}$ . Orientation of the body-fixed frame $\left\{ B \right\}$ with respect to the frame is defined by roll–pitch–yaw angles which are denoted by $\phi ,\theta $ , and $\psi $ , respectively. The quadrotor translational dynamics is obtained in the inertial frame $\left\{ E \right\}$ , while its rotational dynamics is obtained in the body frame $\left\{ B \right\}$ .

The thrust force ( ${T_i}$ ) and the aerodynamic torque ( ${Q_i}$ ) are produced by the four propeller spinning rotors which are given by Eqs. (1) and (2), respectively [Reference Kayacan and Maslim29]:

(1) \begin{align}{T_i} = b{\omega _i}^2 \qquad \quad i = 1, \cdots ,4\end{align}

(2) \begin{align}{Q_i} = d{\omega _i}^2\qquad \quad i = 1, \cdots ,4\end{align}

where ${\omega _i}$ , $b$ , and $d$ are angular speed of the i-th rotor, thrust factor, and drag factor, respectively. According to Eqs. (3)–(6), combination of the thrust forces and the aerodynamic torques produce four generalized forces in the body frame $\left\{ B \right\}$ [Reference Kayacan and Maslim29]:

(3) \begin{align}{U_1} = \mathop \sum \limits_{i = 1}^4 {T_i} = b\left( {\omega _1^2 + \omega _2^2 + \omega _3^2 + \omega _4^2} \right)\end{align}

(4) \begin{align}{U_2} = {T_4}L - {T_2}L = Lb\left( {\omega _4^2 - \omega _2^2} \right)\end{align}

(5) \begin{align}{U_3} = {T_3}L - {T_1}L = Lb\left( {\omega _3^2 - \omega _1^2} \right)\end{align}

(6) \begin{align}{U_4} = \mathop \sum \limits_{i = 1}^4 {\left( { - 1} \right)^{i + 1}}{Q_i} = d\left( {\omega _1^2 - \omega _2^2 + \omega _3^2 - \omega _4^2} \right),\end{align}

where ${U_1}$ is the actuation force along $\hat z$ axis, ${U_2}$ is the actuation torque about $\hat x$ axis, ${U_3}$ is the actuation torque about $\hat y$ axis, and ${U_4}$ is the actuation torque about $\hat z$ axis. Quadrotor’s rotational dynamics is represented in Eq. (7) which is obtained using Euler’s rotational motion equations [Reference Kayacan and Maslim29]:

(7) \begin{align}\left[\begin{array}{c} {\ddot \phi } \\[6pt] {\ddot \theta } \\[4pt] {\ddot \psi } \end{array}\right] = \left[\begin{array}{c} \dfrac{{U_2} - {J_r}\dot \theta {\omega _r} + \left( {{I_{yy}} - {I_{zz}}} \right)\dot \psi \dot \theta }{{I_{xx}}}\\[13pt] \dfrac{{U_3} + {J_r}\dot \phi {\omega _r} + \left( {{I_{zz}} - {I_{xx}}} \right)\dot \psi \dot \phi }{{I_{yy}}}\\[13pt] \dfrac{{U_4} + \left( {{I_{xx}} - {I_{yy}}} \right)\dot \theta \dot \phi }{{I_{zz}}}\end{array}\right],\end{align}

where ${J_r}$ is the rotor moment of inertia and ${\omega _r}$ is the overall speed of rotors which is given by Eq. (8).

(8) \begin{align}{\omega _r} = - {\omega _1} + {\omega _2} - {\omega _3} + {\omega _4}.\end{align}

The translational dynamics of quadrotor is given from Newton approach [Reference Qi, Jiang, Wu, Wang, Wang and Shan25]:

(9) \begin{align}\left[\begin{array}{c} {\ddot X}\\[4pt] {\ddot Y}\\[4pt] {\ddot Z}\end{array}\right] = \left[\begin{array}{c}\dfrac{\left( {{C_\phi }{S_\theta }{C_\psi } + {S_\phi }{S_\psi }} \right){U_1}}{m}\\[13pt]\dfrac{\left( {{C_\phi }{S_\theta }{S_\psi } - {S_\phi }{C_\psi }} \right){U_1}}{m}\\[13pt]\dfrac{\left( {{C_\phi }{C_\theta }} \right){U_1}}{m} - g \end{array}\right],\end{align}

where ${S_.}$ and ${C_.}$ denote ${\rm{sin}}\left(\cdot \right)$ and ${\rm{cos}}\left( \cdot \right)$ functions, respectively. Also, the parameter $g = 9.81{\rm{m}}/{{\rm{s}}^2}$ is the gravity acceleration in the $\hat Z$ direction. Note that effect of aerodynamic torques on the translational motion and external disturbances are neglected. According to Eqs. (7) and (9), the quadrotor is an under-actuated system because it has six DOF but only four actual inputs. As a result, only four DOFs can be controlled independently (i.e., $Z,\psi $ , and one from each pairs of $\left\{ {X,\theta } \right\}$ , $\left\{ {Y,\phi } \right\}$ ). In this research, these four DOFs $\left\{ {X,Y,Z,\psi } \right\}$ are considered to be controlled independently.

2.2. Actuator dynamics

In this study, quadrotor actuators are considered to be brushed direct current (DC) electrical motors without speed reduction gears. Assuming that motor inductance has negligible effect, an approximated first-order equation can be obtained for actuator dynamics:

(10) \begin{align}{\dot \omega _m} = - \frac{d}{{{J_r}}}\omega _m^2 - \frac{{{K_{tm}}{K_{em}}}}{{{R_m}{J_r}}}{\omega _m} + \frac{{{K_{tm}}}}{{{R_m}{J_r}}}{V_a},\end{align}

where, ${\omega _m}$ is angular speed of rotor, ${V_a}$ is supply voltage of motor, ${K_{tm}}$ is motor torque constant, ${K_{em}}$ is electrical constant of motor, ${R_m}$ is motor resistance, ${J_r}$ is rotor inertia, and $d$ is drag factor. Electrical DC motors have a limited input voltage and generate a limited output torque/speed. Consequently, they generate a saturated thrust output which affects motion control of quadrotor in the presence of external disturbances. On the other hand, transient dynamic response of DC motors is very important when control system generates fast and high amplitude speed commands. The input voltage and output velocity limits are included using saturation function in the actuator dynamics.

In general, the actuator dynamics is included to simulate more realistic model of quadrotor. Indeed, it is included to make a challenge with transient response and saturated output thrust of actuators. While it is not considered in the control design, controllers are tuned with the limitations to obtain best possible fast and robust response. This is to ensure that stability and tracking performance are realized with reasonable available power of the actuators.

Figure 2. Architecture of the proposed control system.

3.1. Control system architecture

Architecture of the proposed vision-based closed-loop control system for autonomous landing maneuver of the quadrotor is depicted in Fig. 2. In the presented closed-loop control system, a vision positioning system computes the relative position of the compound AprilTag using the acquired image from virtual reality scene. Relative position is used as the only feedback source for the position controller subsystem and also as an input to the autonomous landing path generator subsystem.

The autonomous landing path generator is a supervisory control algorithm that is realized by a programmed state machine and manages the process of autonomous landing maneuver by simply generating conditional set points.

The quadrotor control subsystem gets relative position feedback and desired set points and generates motors speed commands using the proposed control law. Equations (7) and (9) can be expressed in the closed form of Eq. (12):

(12) \begin{align}\dot X = f\!\left( X \right) + b\!\left( X \right)U + D,\end{align}

with the state vector $X = \left[ {x,y,z,\phi ,\theta ,\psi } \right]$ and the disturbance vector $D$ , the control law is presented in Eq. (13):

(13) \begin{align}U_i^d = \frac{1}{{{b_i}\left( x \right)}}\left( {{\nu _\cdot } - {f_i}\left( x \right)} \right),\end{align}

where a simple feedback linearization term is considered to compensate model nonlinearities and an error dependent term ( ${\nu.}$ ) is included to achieve reliable tracking accuracy. It should be noted that the error dependent term ${\nu _\cdot }$ is the main part of the proposed control law and then is designed as the proposed FOFPID controller. Subsequently, the position and attitude control can be obtained by using the following control laws:

(14) \begin{align}U_y^d & = F{L_y} + m{\nu _y};\qquad\qquad {F{l_y} = 0}\nonumber\\[5pt] U_x^d & = F{L_x} + m{\nu _x}; \qquad\qquad {F{l_x} = 0}\\[5pt] U_1^d & = F{L_z} + {m \over {{C_\phi }{C_\theta }}}{\nu _z};\qquad {F{l_z} = {m \over {{C_\phi }{C_\theta }}}g}\nonumber \end{align}

(15) \begin{align}U_2^d & = F{L_\phi } + {I_{xx}}{\nu _\phi };\qquad {F{l_\phi } = \left( {{I_{zz}} - {I_{yy}}} \right)\dot \psi \dot \theta + {J_r}\dot \theta {\omega _r}}\nonumber\\[6pt] U_3^d & = F{L_\theta } + {I_{yy}}{\nu _\theta };\qquad {F{L_\theta } = \left( {{I_{xx}} - {I_{zz}}} \right)\dot \psi \dot \theta - {J_r}\dot \phi {\omega _r},} \\[6pt] U_4^d & = F{L_\psi } + {I_{zz}}{\nu _\psi };\qquad {F{l_\psi } = ({I_{yy}} - {I_{xx}})\dot \theta \dot \phi }\nonumber \end{align}

where $U_.^d$ is desired generalized force acting on the axis (.), $F{L_.}$ consist feedback linearization part of the control law, and term ${\nu.}$ denotes error dependent part of the control law. Finally, from desired generalized forces, velocity commands of the rotors are given by

(16) \begin{align}\left[ \begin{array}{c}{\omega _1^2}\\[5pt]{\omega _2^2}\\[5pt]{\omega _3^2}\\[5pt]{\omega _4^2}\end{array}\right] = \left[ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c}b & b & b& b\\[4pt]0 & { - Lb} & 0& {Lb}\\[4pt]{ - Lb}& 0 & {Lb}& 0\\[4pt]d& { - d} & d& { - d}\\[4pt]\end{array} \right]^{ - 1}\left[ \begin{array}{c}{U_1^d}\\[5pt]{U_2^d}\\[5pt]{U_3^d}\\[5pt]{U_4^d}\end{array}\right].\end{align}

3.2. FOFPID controller

Structure of the proposed FOFPID controller is shown in Fig. 3. The proposed FOFPID controller is formed as fuzzy fractional-order PI + fuzzy fractional-order PD subsystem with a common core of fuzzy logic controller.

The fuzzy logic controller is a Mamdani-type fuzzy system with two scaled inputs (i.e., scaled tracking error and its scaled fractional-order derivative ( $E,\dot E$ )). The Caputo fractional-order derivative is presented in Eq. (17):

(17) \begin{align}{}_0^CD_t^\mu U\left( t \right) = \frac{1}{{{\rm{\Gamma }}\left( {a - \mu } \right)}}\mathop \smallint \limits_0^t {\left( {t - \tau } \right)^{a - \mu - 1}}{U^{\left( a \right)}}\left( \tau \right)d\tau ,\end{align}

where, $\mu $ is the order of derivative, $a$ is integer supremum or least upper bound of $\mu $ ( $a - 1 \lt \mu \lt a$ ), and $U( t )$ is a continuous time-variant function. The Gamma function is denoted by ${{\Gamma }}({.})$ . Seven fuzzy sets with Gaussian membership function of Eq. (18) are considered for the fuzzy inputs and the fuzzy output in the range between $\left[ { - 1\ {\rm{and}}\ 1} \right]$ .

(18) \begin{align}\mu\! \left( x \right) = {e^{\frac{{ - {{\left( {x - c} \right)}^2}}}{{2{\sigma ^2}}}}}.\end{align}

The membership functions are depicted in Fig. 4. The proposed fuzzy rule base consists of 49 rules, and it is demonstrated in Table I. The fuzzy rule base is generally obtained by try and error in our previous work [Reference Parivash and Ghasemi30]. In the manual tuning procedure, fuzzy rule base is somehow rearranged to obtain smooth and homolographic surface for the fuzzy system output while counteracting the error growth. Abbreviations NL, NM, NS, Z, PS, PM, and PL mean negative large, negative medium, negative small, zero, positive small, positive medium, and positive large, respectively.

Table I. Fuzzy rule base.

Figure 3. Structure of the proposed fuzzy FOPID controller.

Figure 4. Membership function of fuzzy input–output.

Every rule is expressed in the form of Eq. (19), and the output rule firing strength is obtained by using Eq. (20). Kindly, see ref. [Reference Gottwald31] for more detailed descriptions on the fuzzy sets and fuzzy functions.

(19) \begin{align}{\rm{If}}\,E\,{\rm{is\ NL}}\,{\rm{And}}\,\dot E\,{\rm{is}}\,{\rm{NL}}\,{\rm{Then}}\,{U_{FLC}}\,{\rm{is}}\,{\rm{NL}}\end{align}

(20) \begin{align}{\mu _{N{L_E} \cap N{L_{\dot E}}}} = {\rm{min}}\left( {{\mu _{NL}}\left( E \right),{\mu _{NL}}\left( {\dot E} \right)} \right).\end{align}

In the proposed fuzzy logic controller, Mamdani min–max method is considered for inference engine. Consequently, a rule output and interfered output of the whole rules can be obtained by Eqs. (21) and (22), respectively. Kindly, see ref. [Reference Niemiec32] for more detailed descriptions on the Mamdani min–max inference method.

(21) \begin{align}\mu _{N{L^{'}}}^k\left( {{U_{FLC}}} \right) = {\rm{min}}\left({\mu _{N{L_E} \cap N{L_{\dot{E}}}}},{\mu _{NL}}\left( {{U_{FLC}}} \right) \right)\end{align}

(22) \begin{align}{\mu _{FLC}}\left( {{U_{FLC}}} \right) = \max \left\{ {\mu _\cdot ^1, \cdots ,\mu _{N{L^{'}}}^k, \cdots ,\mu _{\cdot}^N} \right\}k = 1, \cdots N.\end{align}

The centroid method is used for defuzzification to obtain the fuzzy system output ( ${U_{FLC}}$ ) according to Eq. (23).

(23) \begin{align}{U_{FLC}} = \frac{{\mathop \smallint \nolimits_U^{} y \cdot {\mu _{FLC}}\left( y \right)dy}}{{\mathop \smallint \nolimits_U^{} {\mu _{FLC}}\left( y \right)dy}}.\end{align}

Finally, output of the proposed fuzzy FOPID controller is given by

(24) \begin{align}{\nu _ \cdot } = \alpha {\rm{U}} + \beta D_{0,t}^{ - \lambda }\left( U \right),\end{align}

where $D_{0,t}^{ - \lambda }\left( U \right)$ denotes a fractional-order integration which is defined in Eq. (25).

(25) \begin{align}D_{0,t}^{ - \lambda }U\!\left( t \right) = \frac{1}{{{{\Gamma }}\!\left( \lambda \right)}}\mathop \int \limits_0^t {\left( {t - \tau } \right)^{\lambda - 1}}U\!\left( \tau \right)d\tau .\end{align}

There are six variables for tuning of the proposed FOFPID controller: variables ${K_e}$ and ${K_d}$ are the scaling factors of the tracking error ( ${e_\cdot}$ ) and its fractional-order derivative ( ${}_0^CD_t^\mu e$ ), respectively. Also, variables $\alpha $ and $\beta $ are proportional and integral gains, respectively. Coefficients $\mu $ and $\lambda $ are derivative and integration fractional order, respectively. Kindly, see refs. [Reference Parivash and Ghasemi33] and [Reference Sharma, Rana and Kumar34] for general stability analysis of the proposed FOFPID formulation.

In this research, it is assumed that The UGV move independently, and there is no communication and collaboration between two vehicles. Consequently, the onboard vision positioning system is the only available feedback source for position control subsystem and supervisory control algorithm.

3.3. Vision positioning system

Vision positioning systems are utilized to obtain the precise position and orientation of a desired object. For this purpose, a certain pattern is placed on the desired object as the target of vision positioning system. To simplify the detection, recognition, and tracking of moving object, the pattern can be a simple color pattern or a complex one. Using the camera parameters and patterns geometry, the pose of the camera with respect to the center of the pattern can be estimated. Pose estimation methods are critically dependent on the accuracy of camera calibration and object geometry [Reference Corke35].

In this paper, the PBVS method is employed for visual position estimation of the landing pad. A compound fiducial marker is introduced to be used as the target of PBVS method. The proposed compound fiducial marker is made of two concentric AprilTag fiducial markers as shown in Fig. 5.

Figure 5. AprilTag fiducial markers.

As shown in Fig. 5(c), the compound AprilTag consists of a large tag and a small tag which are shown separately in Fig. 5(a) and 5(b), respectively. The large AprilTag with a size of 50 × 50 cm is used for detection in far distances, and the small AprilTag with a size of 3 × 3 cm tag is used for near distances and high-precision landing. The small AprilTag, because of its small size, does not distort the overall pattern of the large AprilTag. The proposed compound fiducial marker enables high precision relative positioning of the landing pad using the PBVS method in the range between 10 and 350 cm height. It should be noted that the upper range can be extended easily by using a bigger large AprilTag.

The lower range is very important for autonomous landing maneuver since it implies minimum detectable relative height. Autonomous landing controller has to turn off the engines for a free fall and hopes a successful landing when reaches the minimum detectable relative height. There is a high risk for long free falls during vertical landing since it implies undesirable interaction forces with landing pad or the ground for unsuccessful landings (land outside of the landing pad). Also, long free fall cause inaccurate landing since the landing pad is moving by ground vehicle.

Using the proposed compound fiducial marker with PBVS method provides suitable detection ranges that eliminates the risk of free falls and can be employed by a supervisory landing controller to perform fast, smooth, and accurate landings.

3.4. Supervisory control algorithm

The supervisory control algorithm motivates the quadrotor UAV in the autonomous landing maneuver according to a mission scenario. The proposed supervisory control algorithm is depicted in Fig. 6. In the mission scenario, the quadrotor takes off from the landing pad surface and goes to 3 (m) height and hovers for 5 s. UGV starts moving within these 5 s and goes on an elliptic trajectory defined in Eq. (26).

(26) \begin{align} X\!\left( t \right) & = \left\{\begin{array}{c@{\quad}c} 0 & {0 \le t \le 20}\\[6pt] {2\cos \left( {\dfrac{\pi}{{30}}t - \dfrac{\pi}{2}} \right)} & {20 \lt t}\end{array}\right. \nonumber\\[8pt] Y\!\left(t\right) & = \left\{\begin{array}{c@{\quad}c} 0 & {0 \le t \le 20} \\[6pt] {1 + \sin \left( {\dfrac{\pi}{{30}}t - \dfrac{\pi}{2}} \right)} & {20{\rm{ }} \lt t} \end{array}\right. . \end{align}

Then quadrotor starts following ground vehicle and tries to achieve rendezvous conditions of Eq. (27):

(27) \begin{align}{\mathbb{C}_R} = \left\{\begin{array}{lc}\mathrm{Case}1:\,\mathrm{if}\left\{\begin{array}{c} \left| {{r_x}} \right| \le 0.2\\[3pt] {|{r_y}| \le 0.2}\\[3pt] {|{{\dot{r}}_x}| \le 0.2}\\[3pt] {|{{\dot{r}}_y}| \le 0.2}\\[3pt] {{r_z} \ge 0.5} \end{array}\right. & {\mathrm{Then}\ {Z_d} = {R_z} - 0.10}\\[40pt] \mathrm{Case}2:\,\mathrm{if} \left\{\begin{array}{c}\left| {{r_x}} \right| \le 0.15\\[3pt] {|{r_y}| \le 0.15}\\[3pt] {|{{\dot{r}}_x}| \le 0.15}\\[3pt] {|{{\dot{r}}_y}| \le 0.15}\\[3pt] {0.32 \lt {r_z} \lt 0.5} \end{array}\right. & {\mathrm{Then}\ {Z_d} = {R_z} - 0.01}\\[40pt] \mathrm{Case}3:\,\mathrm{if} \left\{\begin{array}{c} \left| {{r_x}} \right| \le 0.15\\[3pt] {\left| {{r_y}} \right| \le 0.15}\\[3pt] {\left| {{{\dot{r}}_x}} \right| \le 0.15}\\[3pt] {\left| {{{\dot{r}}_y}} \right| \le 0.15}\\[3pt] {{r_z} \le 0.5}\end{array}\right. & {\mathrm{Then}\ {Z_d} = 0.25\ (m)} \end{array}\right.,\end{align}

where, ${r_x}$ , ${r_y}$ , and ${r_z}$ are relative positions obtained from the vision positioning system along X, Y, and Z axes, respectively. Rendezvous conditions are defined on the relative position and velocity states. Whenever rendezvous conditions are satisfied, quadrotor starts vertical landing while tries to preserve rendezvous conditions. If rendezvous conditions are lost, the quadrotor preserves its height and tries to achieve rendezvous condition again. Rendezvous margins and vertical decent rate are decreased going step by step for smooth and accurate landings. When the quadrotor reaches the landing conditions of Eq. (28), supervisory controller commands turn off engines for a vertical touchdown. It should be noted that the final set point for relative height should be considered according to landing gear height for a smooth touchdown.

(28) \begin{align}{\mathbb{C}_L} = {\rm{If}}\left\{\begin{array}{l}\left| {{r_x}} \right| \lt 0.15\ (m)\\[6pt]{\left| {{r_y}} \right| \lt 0.15\ (m)}\ {\rm{Then}}\,{\rm{Turn}}\,{\rm{off}}\,{\rm{engines}}.\\[6pt]{{r_z} \le 0.3\ (m)} \end{array}\right.\end{align}

Figure 6. Proposed supervisory control algorithm.