Abstract
This article researches the improvement of dynamics stability of the ducted fan unmanned aerial vehicles by optimizing its mechanical–structure parameters. The instability phenomenon of ducted fan unmanned aerial vehicles takes place frequently due to the complicated airflow in near-earth space, which easily leads to the stability problems, such as out of control, shaking, and loss accuracy of command tracking. The dynamics equations mirror its dynamics characteristics, which are primarily influenced by the mechanical–structure parameters of the whole system. Based on this, the optimization of mechanical–structure parameters has a significant to improve the dynamics stability of the whole system. Therefore, this article uses the concept of Lyapunov exponents to build the quantification relationship between system’s mechanical–structure parameters and its motion stability to enhance its stability from viewpoint of mechanical–structural parameter design. The takeoff, landing, and hovering stage are respectively studied and the conclusions suggest that the optimization of mechanical–structure parameters can be used to promote dynamics stability.
Keywords
Introduction
In recent years, due to their small sizes, compact layout, and strong maneuverability, the ducted fan unmanned aerial vehicles (UAVs) are regarded as great candidates to perform atmospheric sounding, goods delivery, and other flight tasks. 1,2 However, owing to the atmospheric disturbances, gyroscopic effects, and other unpredictable factors, the ducted fan UAVs are easily subject to shaking and other dynamics stability problem. 3 –5 Therefore, it is particularly important to study its variable–structure parameters to enhance the dynamics stability of the ducted fan UAVs.
Currently, there are two major methods in the control theory. The first is directly to solve dynamics equations. The second is represented by Lyapunov (1892) who published the famous PhD dissertation on general movement stability problems, which marked a new epoch in motion stability studies. The key to applying the Lyapunov stability analysis is to find a Lyapunov function. However, there are no general and constructive methods for deriving Lyapunov functions for nonlinear systems. Therefore, the stability of many nonlinear systems cannot be analyzed. On the other hand, an alternative tool is used to analyze the stability of the nonlinear dynamical system, which is the concept of Lyapunov exponents (LEs). 6 LEs can be described as the average exponential rates of divergence or convergence of nearby orbits in the state space, which can imply system’s stability. When the numerical value of the LE is smaller than 0, the phase orbit of the system will be attracted to a stable fixed point which suggests the whole systems should be stable. The essential characteristics of negative LE indicate that the whole system is dissipative system or nonconservative. With the increasing of negative value, the convergence rate of phase trajectory increases. The system is super stable as the negative value tends to infinity. On the contrary, when the numerical value of the LE is more than 0, the system is chaotic. In addition, the numerical value is close to 0 and the system moves periodically. 7,8
Methods for calculating LEs based on a mathematical model have been well developed and widely used for diagnosing chaotic system. 9 Besides, many complex nonlinear systems have been used for stability using the LEs. For instance, Sadri, Yang and Wu 10 –12 conducted the stability analysis of biped standing using the concept of LEs. Omar et al. 13 improved the robustness using Lyapunov’s stability theory, Zevin and Pinski 14 developed absolute stability criterion for a nonautonomous linear system controlled by a nonlinear feedback control with a time-varying delay based on LEs. The optimization of mechanical–structure parameters is used to investigate a series of accomplishments to improve the stability of whole system. 15,16 Therefore, it is reasonable to use LEs to quantitative the dynamics stability of the ducted fan UAVs to increase dynamics stability via the optimization of the mechanical–structure parameters.
Methodology
Dynamic analysis and modeling
Before constructing a dynamics equation of the system, two reference frames will be defined, the inertial reference frame and the body fixed reference frame. A ground frame could be closed to the inertial reference frame for local flight to derive the equations of the ducted fan UAVs motion. The ground reference frame has its origin fixed to the ground, generally at the starting point, and the three axes constitute north-east-down coordinate. 17 The body fixed reference frame has its origin fixed in the center of gravity of the ducted fan UAVs. The body fixed reference frame is shown in Figure 1. The ϕ, θ, and ψ are the Euler angles; x, y, and z are linear displacement; p, q, and r are angular rates; and u, v, and w are linear rates around x-, y-, and z-axes, respectively. However, some assumptions are regarded as follows: (1) the body is considered as rigid body and (2) ignoring air friction and other tiny resistance.

Simplified vehicle model of the 6-DOF. DOF: degree of freedom.
The six degrees of freedom the ducted fan UAV model has 12 states when Euler angles are used to parameterize attitude. The model can be generated in the form of Poincare’s equation in this case 18
where
where m is the mass of the ducted fan UAVs; g is gravity acceleration; and Ix , Iy , and Iz are the mass moment of inertia. ∑Mx , ∑My , ∑Mz , and ∑Fx , ∑Fy , and ∑Fz are total rolling, pitching, yawing moments, and total forces on the aircraft about body axes in x-, y-, and z-axes, respectively. Moreover, C θ = cos (θ), C −1 θ = sec (θ), Tθ = tan (θ), Cψ = cos (ψ), S ϕ = sin (ϕ), S ψ = sin (ψ), C ϕ = cos (ϕ), and S θ = sin (θ).
Meanwhile, the proportional-derivative (PD) controller will be used to control the ducted fan UAVs when the hovering state. For the PD controller, the torque at each joint is
where θdi is the desired position, kpi and kdi are the proportional and derivative gains.
There are several forces and moments: the gravity, lift, control surface, and fuselage acting on the vehicle, respectively. 19 –21
Gravity
Because the roll attitude ϕ and the pitch attitude θ are not zero, the gravity can be decomposed to three forces along the x-, y-, and z-axes
Airfoil
The lift of airfoil can be generated as
where vi represents the induced velocity via the rotor; ao stand for rotor lift curve coefficient; b denotes the number of the rotor blades; and cr , ρ, Nb , r, ωr , k twist, and CD are the rotor blade chord, air density, number of airfoil, radius of the rotor, angular velocity, twist of the blades, and drag coefficients of the vehicle, respectively.
Control surface
All of the forces and moments acting on the vehicle by the control surface can be expressed as follows
where δ 1, δ 2, and δ 3 are control vane bias, l 1, l 2, and l 3 are the arms of forces, respectively.
Fuselage
As in the hovering, u = 0, v = 0, w = 0, owing to the fuselage aerodynamic drag, the forces and the moments can be generated as
where Cx , Cy , and Cz are the coefficient of lift; Sx , Sy , and Sz are vehicle platform area.
The equations (4) to (11) can be expressed as
Combining the kinematics equation (1), dynamics equation (2) with equations (12) and (13) to obtain the state equations
Numerical calculation of LEs
When equation (14) is combined with equation (15), the nonlinear equations of the ducted fan UAVs can be written as
where X = [q, p] T = (ϕ, θ, ψ, x, y, z, p, q, r, u, v, w) T is state variable.
The calculation of LEs based on dynamics equation can be written as
To calculate LEs spectrum λ 1, λ 2,…, λ 6. Letting time T is equivalent to 0.01 s and iterations k is equivalent to 100. With regard to k th iteration (k = 1, 2,…, k), the initial conditions of the variation equation (16) is set to {v 1 (k−1), v 2 (k−1),…, v 6 (k−1)}. The {v 1 (k−1), v 2 (k−1),…, v 6 (k−1)} is translated into {w 1 (k−1), w 2 (k−1),…, w 6 (k−1)} after the Ts integration. Then, the {w 1 (k−1), w 2 (k−1),…, w 6 (k−1)} is converted into {u 1 (k−1), u 2 (k−1),…, u 6 (k−1)} through the Gram–Schm-orthogonal, the ultimate vector value are the {u 1 k , u 2 k ,…, u 6 k }. Repeat this process and the k chosen large enough. 5
LEs theory can be applied to build the quantitative relationship between the dynamics stability of the whole system and mechanical–structure parameters shown in Figure 2.

Between Lyapunov exponents and structural parameters.
Simulation results
Parameter identification
To obtain numerical dynamics polynomial of the ducted fan UAVs, m, g, Ix , Iy , Iz , ao , b, Nb and so on need to be acquired. Here, mechanical–structure parameters of the vehicle existed in Table 1. In simulation analysis, the Mathematica software is used to develop the dynamics model for the ducted fan UAVs during takeoff, landing, and hovering stage, respectively.
The value of mechanical–structure parameters of the system.
Takeoff stage
With the propeller speed increased slowly, when the total tensile force in the z-axis is larger than its gravity, the ducted fan UAVs display takeoff state. It can been shown that the attitude of the curve y (t), z (t) display raise, and x (t) displays decline, which is shown in Figure 3.

The curve of attitude at takeoff state. (a) ϕ (t), (b) θ(t), (c) ψ(t), (d) x (t), (e) y (t), and (f) z(t).
Landing stage
When the speed of rotor of the vehicle decreases, the system emerges landing state, as shown in Figure 4. The curve x (t), y (t), and z (t) will be decreasing.

The curve of attitude at landing state. (a) ϕ (t), (b) θ (t), (c) ψ (t), (d) x (t), (e) y (t), and (f) z (t).
Seen from Figures 5 and 6, it can be drawn that the convergence rate of the LEs spectrum to zero or negative during the landing stage is much faster than that of takeoff stage, which shows that the dynamic stability of the ducted fan UAVs at landing state is better than at takeoff. It can be used to interpret as the stable aircraft can land even when the engines have no power, only control surfaces are also enough for landing.

Lyapunov exponents spectrum during takeoff state. (a) λ (ϕ), (b) λ (θ), (c) λ (ψ), (d) λ (x), (e) λ (y), and (f) λ (z).

Lyapunov exponents spectrum during landing stage. (a) λ (ϕ), (b) λ (θ), (c) λ (ψ), (d) λ (x), (e) λ (y), and (f) λ (z).
For the parametric diversity, based on the mechanical–structural stability analysis, significant insights into the effects of the radius of rotor on the system have been studied. Other mechanical–structure parameters are not unchanged and the value of the rotor radius r is enhanced to 0.6 m (the modified value r is equivalent to 0.6 m, and the original value is equivalent to 0.4 m). Modified the curve of attitude and LEs spectrum shown in Figures 7 and 8, respectively. Looking from Figure 7, it has been revealed that the attitude of the curve x (t), y (t), and z (t) will be greatly decreased compared with Figure 5. Compared with Figures 6 and 8, it can be discover that the value of LEs spectrum at the modified landing state to negative values are larger than during the original landing state. Therefore, changing mechanical–structure parameters can be used to improve its dynamics stability.

Modified the curve of attitude at landing state. (a) ϕ (t), (b) θ (t), (c) ψ (t), (d) x (t), (e) y (t), and (f) z (t).

Modified Lyapunov exponents spectrum during landing stage. (a) λ (ϕ), (b) λ (θ), (c) λ (ψ), (d) λ (x), (e) λ (y), and (f) λ (z).
Hovering stage
In this section, by setting the θd 1 = ϕ = 0.7sin(2ϕ), θd 2 = θ = 0.7sin(2θ), and θd 3 = ψ = 0.7sin(2ψ), the controller gains were chosen as the kp 1 = 170 (N/rad), kd 1 = 40 (NS/rad); kp 2 = 90 (N/rad), kd 2 = 10 (NS/rad); and kp 3 = 110 (N/rad), kd 3 = 10 (NS/rad) in the simulations. Figure 10(a) to (c) reveals that the three LEs take about 40 s converge to negative values, indicating the system has been successfully driven to the steady state. Therefore, the controller gain are reasonable. Figure 9(a) to (c) shows that the attitude of the curve ϕ (t), θ (t), and ψ (t) maintains the level trajectory (the angle is approximately equal to 0) after around 3 s (the dash lines are the desired trajectories and the solid lines are the actual trajectories).

The curve of attitude at hovering state. (a) ϕ (t), (b) θ (t), (c) ψ (t), (d) x (t), (e) y (t), and (f) z (t).

Lyapunov exponents spectrum during hovering state. (a) λ (ϕ), (b) λ (θ), (c) λ (ψ), (d) λ (x), (e) λ (y), and (f) λ (z).
Meanwhile, letting the θd 4 = x = 1.7sin (2x), θd 5 = y = 1.7sin (2y), and θd 6 = z = 1.7sin (2z); also, the controller gains are selected as kp 4 = 1 (N/rad), kd 4 = 300(NS/rad); kp 5 = 2 (N/rad), kd 5 = 350 (NS/rad); and kp 6 = 2 (N/rad), kd 6 = 250 (NS/rad), it has been shown that all of LEs are negative values, suggesting the system is stable, as shown in Figure 10(d) to (f). Figure 9(d) to (f) also shows that the linear displacement can successfully keep the horizontal trajectory (the linear displacement is approximately equal to 0). On the above analysis in this section point that LE spectrum takes to converge to 0 or negative values, suggesting that the system has been successfully driven to the steady state.
The initial angular displacement was changed into 10−1 (rad) with the same initial angular velocities p, q, and r are equivalent to 0. The results of the attractor are not changed significantly, which demonstrates that the system is not sensitive to the initial condition. This result denotes the system is stable, as shown in Figure 11(a) to (c) (the red line is original stage, and the blue line is modified initial angular displacement).

Attractor of the state space. (a) (ϕ-p), (b) (θ-q), and (c) (ψ-r).
Experimental results
In this section, the simulation results of different stages are revealed. The ducted fan UAV’s flight experience is used to Pixhawk’s flight control development environment, as shown in Figure 12, where it can monitor the flight status of the aircraft, data storage, and so on. When the value of rotor radius r is equivalent to 0.6 m, the dynamics stability of the vehicle is much better than that of r is equivalent to 0.4 m in the case of the landing state.

Online interface of the Pixhawk’s flight control development environment.
Figures 13 to 15 are derived from the online interface of the Pixhawk’s flight control and Figures 16 to 18 are linear displacement parameters in x-, y-, and z-axes, respectively It can be easily observed that the attitude angle of the ducted fan UAVs shows much better convergence to 0 in the value of rotor radius r which is equivalent to 0.6 m than that of r which is equivalent to 0.4 m according to Figures 16 to 18 in landing stage. This conclusion can be used to verify the compared results of Figure 6 with Figure 8. The ducted fan UAV’s stability in landing is better than in takeoff according to the comparison of Figures 4 and 6. Figures 19 to 21 are also obtained from the online interface of the Pixhawk’s flight control development environment of the ducted fan UAV, which show the three Eluer angular of the system in x-, y-, and z-axes, respectively. The experimental results show that the ducted fan UAV has much better flight performance in the landing stage than in the takeoff stage, which can be proved in the results of Figures 4 and 6. Figures 19 to 21(a) show that the attitude of the roll, pitch, and yaw maintain the level trajectory (the angle is approximately equal to 0) over time than Figures 19 to 21(b). In summary, the optimization of mechanical–structure parameters can be used to promote dynamics stability.

Experimental attitude of the curve location parameter of in the x-axis during landing state. (a) r = 0.4 m and (b) r = 0.6 m.

Experimental attitude of the curve location parameter of in the y-axis during landing state. (a) r = 0.4 m and (b) r = 0.6 m.

Experimental attitude of the curve location parameter of in the z-axis during landing state. (a) r = 0.4 m and (b) r = 0.6 m.

Experimental attitude of the curve angular in the x-axis during landing state. (a) r = 0.4 m and (b) r = 0.6 m.

Experimental attitude of the curve angular in the y-axis during landing state. (a) r = 0.4 m and (b) r = 0.6 m.

Experimental attitude of the curve angular in the z-axis during landing state. (a) r = 0.4 m and (b) r = 0.6 m.

The curve angular in the x-axis. (a) Roll angular during takeoff stage and (b) roll angular during landing stage.

The curve angular in the y-axis. (a) Pitch angular during takeoff stage and (b) pitch angular during landing stage.

The curve angular in the y-axis. (a) Yaw angular during takeoff stage and (b) yaw angular during landing stage.
Conclusions
The stability analysis is very important for the development of the ducted fan UAVs. This article has proposed to use the concept of LEs to establish the quantification relationship between the ducted fan UAV’s structural parameters and its motion stability to optimize structural parameters for increased dynamics stability. This concept is one of the most powerful tools for analyzing complex nonlinear dynamical system. Furthermore, LEs are “invariant” measures of dynamics characteristics, that is, independent of the initial condition. However, existing methods are either not constructed, for example, Lyapunov’s second methods, or sensitive to the individual testing or initial conditions, such as simulation-based methods. The proposed method via the concept of LEs is constructive and theoretically sound in this article, owing to the exponents are “invariant” measures of system dynamics.
Three different stage of stability, such as takeoff, landing, and hovering stage, respectively, have been rigorously analyzed in this article using the concept of LEs for the ducted fan UAVs. In the case of landing stage stability analysis, the radius range of rotor where the vehicle stability is guaranteed has been changed. It has shown significant insights into the effects of mechanical–structure parameters on dynamics stability. Moreover, the state space is another important property of a dynamics system. To investigate the effect of the initial condition on state space. The PD controller has been presented for the hovering state. The disturbance of system stability comes from initial conditions. It can be derived that the vehicle is stable and will return to its stable fixed point, which verifies the stable under the suitable PD controller. The conclusion is the optimization of mechanical–structure parameters has exerted significant effects on its stability.
The main contribution of this article lies in the fact that it starts from symbol dynamic modeling to stability analysis using the concept of LEs. The optimization of structure parameters is of great significance to improve the movement stability of the ducted fan UAVs. Therefore, this article builds the quantitative relationship between the structural parameter of the ducted fan aerial vehicles and their system stability based on LEs theory from viewpoint of structural parameters design. Moreover, the reliability and stability of ducted fan UAVs under service could be promoted by optimizing the mechanical–structure parameters.
As discussed earlier, the presented method for calculating LEs is applicable to general dynamics systems provided that the mathematical models are available. Therefore, the proposed method can be extended to more complex models, especially to multi-body dynamics system, which is the subject of future work.
Footnotes
Acknowledgements
We thank the anonymous reviewers for helpful and insightful remarks. In addition, I want to thank, in particular, the patience, care, and support from XX Chen over the passed years.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research is supported by the National Key Research and Development Plan, the On-line monitoring technology and equipment of contact acoustic flow (Grant No. 2017YFC0405702), and central public welfare basic scientific research institutes special funds (Grant No. Y919008).
