2.2.2. Aircraft-Body Coordinate Frame *OBxByBz<sup>B</sup>*

The selected mass center of the aircraft is the origin *O<sup>B</sup>* of the coordinates. The coordinate system is fixed to the aircraft body, and the axis *x<sup>B</sup>* is along the axis of the aircraft's symmetry plane, which points to the nose of the aircraft. *y<sup>B</sup>* points to the starboard side of the aircraft, while *z<sup>B</sup>* is perpendicular to *x<sup>B</sup>* in the aircraft's symmetry plane, which points to the bottom of the body.

The Earth-surface reference frame was converted to the aircraft-body coordinate frame as follows:

$$\begin{pmatrix} R\_E^B = \begin{pmatrix} \cos\theta\cos\psi & \cos\theta\sin\psi & -\sin\theta \\ \sin\phi\sin\theta\cos\psi - \cos\phi\sin\psi & \sin\phi\sin\theta\sin\psi + \cos\phi\cos\psi & \sin\phi\cos\theta \\ \cos\phi\sin\theta\cos\psi + \sin\phi\sin\psi & \cos\phi\sin\theta\sin\psi - \sin\phi\cos\psi & \cos\phi\cos\theta \end{pmatrix} \tag{24}$$

#### 2.2.3. Velocity of Aircraft Relative to the Air

The velocity of the aircraft relative to the air is defined as follows:

$$
\stackrel{\rightarrow}{V}\_{a} = \begin{pmatrix} V\_{ax} \\ V\_{ay} \\ V\_{az} \end{pmatrix} = \begin{pmatrix} V\_{G^x}^E - V\_{wind^x}^E \\ V\_{G^y}^E - V\_{wind^y}^E \\ V\_{G^z}^E - V\_{wind^z}^E \end{pmatrix} \tag{25}
$$

where → *V<sup>a</sup>* is the airspeed of the aircraft, *V E <sup>G</sup> <sup>x</sup>*, *V E <sup>G</sup> <sup>y</sup>*, *V E <sup>G</sup> <sup>z</sup>* is the relative velocities to the Earth on the axis *x*, *y*, *z*, and *V E windx*, *V E windy*, *V E windz* is the relative wind speed to the Earth on the axis *x*, *y*, *z*.

#### 2.2.4. Angle of Attach and Sideslip Angle

The aircraft should make its wings fly at a positive angle with respect to the airspeed vector in order to rise. The positive angle is the angle of attack, noted as *α*; the angle between the velocity vector and the plane *xBz<sup>B</sup>* is deemed as the sideslip angle, noted as *β*.

The velocity under the aircraft-body coordinate frame was adopted to calculate the angle of attach and the sideslip angle:

$$\begin{cases} \alpha = \arctan\left(\frac{V\_a^B}{V\_{ax}^B}\right) \\ \beta = \arcsin\left(\frac{V\_{ay}^B}{\left|\frac{\overrightarrow{\cdot}}{V\_a}\right|}\right) \end{cases} \tag{26}$$

#### 2.2.5. Force and Moment

Based on aerodynamic principles, forces and moments acting on a DEP aircraft can be summarized as follows:

$$
\begin{pmatrix} F\_x \\ F\_y \\ F\_z \end{pmatrix} = \begin{pmatrix} -mg\sin\theta \\ mg\cos\theta\sin\phi \\ mg\cos\theta\cos\phi \end{pmatrix} + \begin{pmatrix} -\frac{1}{2}\rho S V\_a^2 \mathbb{C}\_x \\ -\frac{1}{2}\rho S V\_a^2 \mathbb{C}\_y \\ -\frac{1}{2}\rho S V\_a^2 \mathbb{C}\_z \end{pmatrix} + \begin{pmatrix} F\_{P\_x} \\ F\_{P\_y} \\ F\_{P\_z} \end{pmatrix} \tag{27}
$$

$$
\begin{pmatrix} M\_x \\ M\_y \\ M\_z \end{pmatrix} = \begin{pmatrix} \frac{1}{2}\rho SbV\_a^2\mathcal{C}\_l \\ \frac{1}{2}\rho S\overline{c}V\_a^2\mathcal{C}\_m \\ \frac{1}{2}\rho SbV\_a^2\mathcal{C}\_n \end{pmatrix} + \begin{pmatrix} T\_{P\_x} \\ T\_{P\_y} \\ T\_{P\_z} \end{pmatrix} \tag{28}
$$

where *m* is the mass of aircraft, *ρ* is the air density, *S* is the wing area for reference, *V<sup>a</sup>* is the norm of the airspeed vector → *Va*, *b* is the wingspan for reference, and *c* is the mean aerodynamic wing chord. *F<sup>P</sup>* and *T<sup>P</sup>* are the thrust and torque generated by the thrusters above. *C<sup>x</sup>* represents the drag coefficient, *C<sup>y</sup>* is the lateral force coefficient, and *C<sup>z</sup>* is the lift coefficient. *C<sup>l</sup>* is the roll moment coefficient, *C<sup>m</sup>* is the pitch moment coefficient, and *C<sup>n</sup>* is the yaw moment coefficient, as calculated and shown below.

The aerodynamic coefficient is as follows:

$$\begin{cases} \mathsf{C}\_{x} = \mathsf{C}\_{x0} + \mathsf{C}\_{xa}\mathfrak{a} + \mathsf{C}\_{xq}\frac{\mathsf{c}}{2\mathsf{V}\_{a}}q + \mathsf{C}\_{x\delta\_{E}}\Big|\delta\_{E}\Big|\\\mathsf{C}\_{y} = \mathsf{C}\_{y\beta}\mathfrak{f} + \mathsf{C}\_{y\eta}\frac{\mathsf{c}}{2\mathsf{V}\_{a}}p + \mathsf{C}\_{yr}\frac{\mathsf{c}}{2\mathsf{V}\_{a}}r + \mathsf{C}\_{y\delta\_{A}}\delta\_{A} + \mathsf{C}\_{y\delta\_{R}}\delta\_{R}\\\ \mathsf{C}\_{z} = \mathsf{C}\_{z0} + \mathsf{C}\_{za}\mathfrak{a} + \mathsf{C}\_{zq}\frac{\mathsf{c}}{2\mathsf{V}\_{a}}q + \mathsf{C}\_{z\delta\_{E}}\delta\_{E} \end{cases} \tag{29}$$

The pneumatic moment coefficient is as follows:

$$\begin{cases} \mathsf{C}\_{l} = \mathsf{C}\_{l\not\beta}\mathfrak{f} + \mathsf{C}\_{l\not p}\frac{\mathsf{c}}{\overline{2V\_{a}}}p + \mathsf{C}\_{l r}\frac{\mathsf{c}}{\overline{2V\_{a}}}r + \mathsf{C}\_{l\delta\_{A}}\delta\_{A} + \mathsf{C}\_{l\delta\_{R}}\delta\_{R} \\ \mathsf{C}\_{m} = \mathsf{C}\_{m0} + \mathsf{C}\_{m a}a + \mathsf{C}\_{m q}\frac{\mathsf{c}}{\overline{2V\_{a}}}q + \mathsf{C}\_{m\delta\_{E}}\delta\_{E} \\ \mathsf{C}\_{n} = \mathsf{C}\_{n\beta}\mathfrak{f} + \mathsf{C}\_{n p}\frac{\mathsf{c}}{\overline{2V\_{a}}}p + \mathsf{C}\_{n r}\frac{\mathsf{c}}{\overline{2V\_{a}}}r + \mathsf{C}\_{n\delta\_{A}}\delta\_{A} + \mathsf{C}\_{n\delta\_{R}}\delta\_{R} \end{cases} \tag{30}$$

In the abovementioned equation, coefficients such as *Cx*0, *Cxα*, *Cxq*, *Cxδ<sup>E</sup>* are derived from the partial derivatives in a Taylor series approximation process and are dimensionless values, which are determined by the aircraft's parameters.

#### 2.2.6. Flight Dynamics Equations

From the momentum theorem, the following can be obtained:

$$
\begin{pmatrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{pmatrix} = \frac{1}{m} \begin{pmatrix} F\_x \\ F\_y \\ F\_z \end{pmatrix} + R\_E^B \begin{pmatrix} 0 \\ 0 \\ g \end{pmatrix} - \begin{pmatrix} 0 & -r & q \\ r & 0 & -p \\ -q & p & 0 \end{pmatrix} \begin{pmatrix} u \\ v \\ w \end{pmatrix} \tag{31}
$$

From the moment of momentum theorem, the following can be obtained:

$$
\begin{pmatrix}
\dot{p} \\
\dot{q} \\
\dot{r}
\end{pmatrix} = I^{-1} \begin{pmatrix}
M\_{\mathcal{X}} \\
M\_{\mathcal{Y}} \\
M\_{z}
\end{pmatrix} - I^{-1} \begin{pmatrix}
0 & -r & q \\
r & 0 & -p \\
\end{pmatrix} I
\tag{32}
$$

In the equations, *I* is the moment of inertia of the aircraft:

$$I = \begin{pmatrix} I\_{xx} & -I\_{xy} & -I\_{xz} \\ -I\_{yx} & I\_{yy} & -I\_{yz} \\ -I\_{zx} & -I\_{zy} & I\_{zz} \end{pmatrix} \tag{33}$$

The set of supplementary kinematic equations is presented as follows:

$$
\begin{pmatrix}
\dot{\phi} \\
\dot{\theta} \\
\dot{\psi}
\end{pmatrix} = \begin{pmatrix}
1 & \sin\phi\tan\theta & \cos\phi\tan\theta \\
0 & \cos\phi & -\sin\phi \\
0 & \sin\phi\sec\theta & \cos\phi\sec\theta
\end{pmatrix} \begin{pmatrix}
p \\ q \\ r
\end{pmatrix} \tag{34}
$$

Equations (31), (32), and (34) are the six-degree-of-freedom flight dynamics equations of DEP aircrafts, and the flight state of an aircraft can be obtained by solving the above equations. *Drones* **2022**, *6*, x FOR PEER REVIEW 10 of 25

Simulation was carried out in order to verify the correctness of the mathematical model of DEP aircrafts. The inputs of the model were the target roll angle and pitch angle, and the aircraft was controlled to fly in a steady state with zero sideslip angle. The attitude response inputs of the aircraft are shown in Figures 3 and 4. *Drones* **2022**, *6*, x FOR PEER REVIEW 10 of 25

**Figure 3.** The pitch angle response of the DEP aircraft. **Figure 3.** The pitch angle response of the DEP aircraft. **Figure 3.** The pitch angle response of the DEP aircraft.

**Figure 4.** The roll angle response of the DEP aircraft. **Figure 4.** The roll angle response of the DEP aircraft.

**Figure 4.** The roll angle response of the DEP aircraft.

ing it had a good control effect.

ing it had a good control effect.

angle of the aircraft changed accordingly. As the roll channel is coupled with the yaw channel, in order to ensure zero sideslip angle flight, the yaw channel must respond to meet the control requirements. Figure 5 shows the response curve of the aircraft's sideslip angle. It can be seen that the sideslip angle caused by the roll channel only changed slightly. The aircraft returned to the steady flight with zero sideslip angle quickly, mean-

channel, in order to ensure zero sideslip angle flight, the yaw channel must respond to meet the control requirements. Figure 5 shows the response curve of the aircraft's sideslip angle. It can be seen that the sideslip angle caused by the roll channel only changed slightly. The aircraft returned to the steady flight with zero sideslip angle quickly, mean-

At the fifth second, the input command of the roll angle changed, and the sideslip angle of the aircraft changed accordingly. As the roll channel is coupled with the yaw channel, in order to ensure zero sideslip angle flight, the yaw channel must respond to meet the control requirements. Figure 5 shows the response curve of the aircraft's sideslip angle. It can be seen that the sideslip angle caused by the roll channel only changed slightly. The aircraft returned to the steady flight with zero sideslip angle quickly, meaning it had a good control effect. *Drones* **2022**, *6*, x FOR PEER REVIEW 11 of 25

**Figure 5.** The sideslip angle response of the DEP aircraft.

#### **3. Coordinated Thrust Control and Fault-Tolerant Control of the DEP Aircraft**

#### **3. Coordinated Thrust Control and Fault-Tolerant Control of the DEP Aircraft** *3.1. Coordinated Thrust Control of the DEP Aircraft*

**Figure 5.** The sideslip angle response of the DEP aircraft.

#### *3.1. Coordinated Thrust Control of the DEP Aircraft* 3.1.1. Longitudinal Control Loop

3.1.1. Longitudinal Control Loop In this study, a total energy control system (TECS) was designed for the longitudinal control of DEP aircrafts, which controls the entire flight of climb, cruise and descent with the best goal of minimizing the aircraft's energy consumption [24]. The system solves the coupling problems concerning the power lever angle and the elevator. Based on the total aircraft energy, the throttle directly corresponds to the increase or decrease in the overall aircraft energy, the rise and fall directly correspond to the distribution of the aircraft's kinetic and potential energy, and the altitude and airspeed are the results produced by In this study, a total energy control system (TECS) was designed for the longitudinal control of DEP aircrafts, which controls the entire flight of climb, cruise and descent with the best goal of minimizing the aircraft's energy consumption [24]. The system solves the coupling problems concerning the power lever angle and the elevator. Based on the total aircraft energy, the throttle directly corresponds to the increase or decrease in the overall aircraft energy, the rise and fall directly correspond to the distribution of the aircraft's kinetic and potential energy, and the altitude and airspeed are the results produced by the joint action of the power lever angle and the elevator. TECS was derived as follows:

$$E\_{\text{tot}} = E\_{\text{kin}} + E\_{\text{pot}} = \frac{1}{2} m V\_a^2 + mgH \tag{35}$$

$$\frac{\dot{E}\_{\text{tot}}}{\text{mg}} = \frac{m V\_a \dot{V}\_a}{\text{mg}} + \frac{\dot{H} \text{mg}}{\text{mg}} = \frac{V\_a \dot{V}\_a}{\text{g}} + \dot{H} \tag{36}$$

$$\dot{E}\_{\text{spec}} = \frac{\dot{E}\_{\text{tot}}}{\text{mg}V\_a} = \frac{\dot{V}\_a}{\text{g}} + \frac{\dot{H}}{V\_a} = \frac{\dot{V}\_a}{\text{g}} + \sin \gamma \approx \frac{\dot{V}\_a}{\text{g}} + \gamma \tag{37}$$

$$
\dot{E}\_{\text{dist}} = \gamma - \frac{\dot{V}\_a}{\mathcal{S}} \tag{38}
$$

(39)

*g* In the equations, *Etot* is the total energy of the aircraft, *Ekin* is the kinetic energy of the aircraft, *Epot* is the potential energy of the aircraft, *m* is the aircraft mass, *H* is the altitude, *g* is the gravitational acceleration, *V<sup>a</sup>* is the airspeed, is the track angle, and *Edist* is the specific energy distribution rate. In addition, the incremental thrust *T<sup>c</sup>* In the equations, *Etot* is the total energy of the aircraft, *Ekin* is the kinetic energy of the aircraft, *Epot* is the potential energy of the aircraft, *m* is the aircraft mass, *H* is the altitude, *g* is the gravitational acceleration, *V<sup>a</sup>* is the airspeed, *γ*is the track angle, and . *Edist* is the specific energy distribution rate. In addition, the incremental thrust ∆*T<sup>c</sup>* is associated with the specific energy gradient . *Espec*.

The TECS approach connects the change in commanded thrust to the change in spe-

.

*TI c TP spec K T K E s*

is associated with the specific energy gradient *Espec*

cific energy rate as follows:

The TECS approach connects the change in commanded thrust to the change in specific energy rate as follows: *Drones* **2022**, *6*, x FOR PEER REVIEW 12 of 25

$$
\Delta T\_c = \left(K\_{TP} + \frac{K\_{TI}}{s}\right) \dot{E}\_{spec} \tag{39}
$$

In the above equation in *s* domain, *KTP* is the proportional gain of the thrust control loop, while *KT I* is the integral gain of the thrust control loop that drives the steady-state error to zero. It was assumed that the elevator control is under energy conservation and the elevator can convert kinetic energy to potential energy, so the specific energy distribution rate is presented as follows: . In the above equation in *s* domain, *KTP* is the proportional gain of the thrust control loop, while *KTI* is the integral gain of the thrust control loop that drives the steadystate error to zero. It was assumed that the elevator control is under energy conservation and the elevator can convert kinetic energy to potential energy, so the specific energy distribution rate is presented as follows:

$$
\dot{E}\_{\text{dist}} = \gamma - \frac{V\_a}{8} \tag{40}
$$

Based on that, changes in pitch angle command ∆*θ<sup>c</sup>* are related to changes in . *Edist*: Based on that, changes in pitch angle command are related to changes in *Edist* :

$$
\Delta\theta\_{\rm c} = \left(K\_{EP} + \frac{K\_{EI}}{s}\right) \dot{E}\_{dist} \tag{41}
$$

*c*

where *KEP* is the proportional gain of the pitch angle control loop, and *KEI* is the integral gain of the pitch angle control loop. where *KEP* is the proportional gain of the pitch angle control loop, and *KEI* is the integral gain of the pitch angle control loop.

The aircraft's thrust is associated with the thrust command, and the change of elevator deflection angle ∆*δ<sup>e</sup>* is related to the pitch command: The aircraft's thrust is associated with the thrust command, and the change of elevator deflection angle *<sup>e</sup>* is related to the pitch command:

$$
\Delta T = \mathcal{G}\_{\text{eff}}(\mathbf{s}) \Delta T\_{\text{c}} \,\,\Delta \delta\_{\text{c}} = \mathcal{G}\_{\text{elev}}(\mathbf{s}) \Delta \theta\_{\text{c}} \tag{42}
$$

where *Gthr*(*s*) denotes the combined thrust control function, and *Gelev*(*s*) is the combined pitch control and elevator actuator dynamics function. Based on the above derivation, the functional block diagram of TECS can be represented as in Figure 6. where ( ) *G s thr* denotes the combined thrust control function, and ( ) *G s elev* is the combined pitch control and elevator actuator dynamics function. Based on the above derivation, the functional block diagram of TECS can be represented as in Figure 6.

**Figure 6.** Functional block diagram of TECS. **Figure 6.** Functional block diagram of TECS.

#### 3.1.2. Lateral Control Loop 3.1.2. Lateral Control Loop

The total heading control system (THCS) is leveraged for the lateral control of DEP aircrafts [22]. The error signals between the commanded and actual rate of change of heading (*c* and ) and between the commanded and actual rate of change of sideslip ( *c* and ) are computed as follows: The total heading control system (THCS) is leveraged for the lateral control of DEP aircrafts [24]. The error signals between the commanded and actual rate of change of heading ( . *ψc* and . *ψ*) and between the commanded and actual rate of change of sideslip ( . *βc* and . *β*) are computed as follows:

$$
\Delta \dot{\psi} = \dot{\psi}\_c - \dot{\psi} \tag{43}
$$

$$
\Delta \dot{\beta} = \dot{\beta}\_c - \dot{\beta} \tag{44}
$$

The commanded roll angle changes *c* and the yaw rate changes *<sup>c</sup> r* based on The commanded roll angle changes ∆*φ<sup>c</sup>* and the yaw rate changes ∆*r<sup>c</sup>* based on these errors are calculated as follows:

*c*

$$
\Delta\phi\_c = \frac{V\_a}{\mathcal{g}} \left( K\_{RP} + \frac{K\_{RI}}{s} \right) \left( \Delta\dot{\psi} + \Delta\dot{\beta} \right) \tag{45}
$$

$$
\Delta r\_c = \frac{V\_a}{\mathcal{S}} \left( K\_{YP} + \frac{K\_{YI}}{s} \right) \left( \Delta \dot{\psi} - \Delta \dot{\beta} \right) \tag{46}
$$

where *KRP*, *KRI* is the proportional gain and integral gain of the roll angle control loop, and *KYP*, *KY I* is the proportional gain and integral gain of the yaw rate control loop.

In terms of aircraft-related components, deflection changes in ailerons ∆*δ<sup>a</sup>* and the deflection changes in rudder ∆*δ<sup>r</sup>* are calculated in response to roll angle commands and roll angle speed variations, respectively.

$$
\Delta \delta\_{\mathfrak{a}} = \mathcal{G}\_{\text{ail}}(\mathbf{s}) \Delta \phi\_{\mathfrak{c}} \,\,\Delta \delta\_{\mathfrak{r}} = \mathcal{G}\_{\text{rad}}(\mathbf{s}) \Delta r\_{\mathfrak{c}} \tag{47}
$$

where *Gail*(*s*) and *Grud*(*s*) are ailerons and the rudder controller and the actuator dynamics function, respectively.

#### *3.2. Fault Response Strategy and Fault-Tolerant Control of DEP Aircrafts*

This study focused on stuck and failed thrusters. Causes of thruster fault include decreased gain of a brushless motor due to aging of the motor stator coils, excessive friction of the motor rotor's shaft, and degradation of the motor's magnet performance, leading to the output deviating from the normal one. Macroscopically, when the output of a brushless motor is weak during the actual flight, changes in attitude angle of the motor is reduced with the same control amount, and the entire aircraft becomes "sluggish". In this study, 16 electric thrusters were adopted for the model object, with a symmetric distribution of eight thrusters on the left and eight on the right. The thruster near the center was numbered 1, and the outermost thruster was numbered 8. The state matrices of the thrusters on the left and the right were expressed by *X<sup>L</sup>* and *XR*. In the preliminary design, the total thrust of the system is given by Equation (50), assuming that the thrust of all thrusters on the same side is equal [25].

$$T\_{\mathbb{R}\_{-}\dot{\iota}} = T\_{\mathbb{R}\_{-}\dot{\jmath}\_{\prime}} \forall i \in [1, 8], j \in [1, 8] \tag{48}$$

$$T\_{L\_{\text{-}}j} = T\_{L\_{\text{-}}j\prime} \forall i \in [1, 8] , j \in [1, 8] \tag{49}$$

$$T\_{total} = T\_L \sum\_{i=1}^{8} X\_L(i) + T\_R \sum\_{i=1}^{8} X\_R(i) \tag{50}$$

Through this thruster counting method, the total yaw moment provided by this propulsion system is given by Equation (51), where the diameter of each thruster is *D*:

$$M\_{\rm tot} = T\_L \frac{D}{2} \sum\_{i=1}^{8} (2i - 1) X\_L(i) - T\_R \frac{D}{2} \sum\_{i=1}^{8} (2i - 1) X\_R(i) \tag{51}$$

Solving *T<sup>R</sup>* in (50) and (51), *TLe f t* and *TRight* can be obtained as shown in Equations (53) and (54) below:

$$M\_{tot} = T\_L \frac{D}{2} \sum\_{i=1}^{8} (2i - 1) X\_L(i) - (\frac{T\_{tot}}{\frac{8}{2}} - \frac{\sum\_{i=1}^{8} X\_L(i)}{\frac{8}{2}}) \frac{D}{2} \sum\_{i=1}^{8} (2i - 1) X\_R(i) \tag{52}$$

$$T\_{Leff} = \frac{\frac{2}{\Sigma}M\_{tot} + \frac{T\_{lat}}{\frac{\Sigma}{\Sigma}X\_R(i)}\sum\_{i=1}^{8}(2i-1)X\_R(i)}{\frac{8}{\Sigma}(2i-1)[X\_L(i)] + \frac{\frac{8}{\Sigma}X\_L(i)}{\frac{8}{\Sigma}X\_R(i)}\sum\_{i=1}^{8}(2i-1)[X\_R(i)]}}\tag{53}$$

$$T\_{Rlight} = \frac{T\_{tot}}{\sum\_{i=1}^{8} X\_R(i)} - \frac{\frac{2}{\Sigma} M\_{tot} + \frac{T\_{tot}}{\sum\_{i=1}^{8} X\_R(i)} \sum\_{i=1}^{8} (2i - 1) X\_R(i)}{\frac{8}{\frac{8}{\pi}} X\_R(i)}\tag{54}$$

*TLe f t* and *TRight* of the above equations were input as thrust commands to the electric thrusters on both sides, where the state matrices of the thrusters were considered for monitoring the minimum thrust demand, turbine engine state, generator state, power bus state, and electric thruster state. If any of the components fail, the corresponding variable in the state matrices degrades to 0. In addition to enabling coordinated control of the electric thrusters on both sides, the thrust can be redistributed to maintain stability and maneuverability of the aircraft in case of a component fault. The fault-tolerant controller of the DEP system designed in this study features a fault injection module. The function developed so far allows the remaining thrusters to make corresponding changes to recover the aircraft's thrust to the prefault level when a single thruster on the left/right fails and the torque and the rotational speed fail to reach the normal operational level. thrusters on both sides, where the state matrices of the thrusters were considered for monitoring the minimum thrust demand, turbine engine state, generator state, power bus state, and electric thruster state. If any of the components fail, the corresponding variable in the state matrices degrades to 0. In addition to enabling coordinated control of the electric thrusters on both sides, the thrust can be redistributed to maintain stability and maneuverability of the aircraft in case of a component fault. The fault-tolerant controller of the DEP system designed in this study features a fault injection module. The function developed so far allows the remaining thrusters to make corresponding changes to recover the aircraft's thrust to the prefault level when a single thruster on the left/right fails and the torque and the rotational speed fail to reach the normal operational level. When the *i* th thruster fails, the mathematical form of the rotational speed of the

*Drones* **2022**, *6*, x FOR PEER REVIEW 14 of 25

1 1

 

*R*

*tot*

ight 8 8

*X i*

*T*

*i*

*R*

*T*

where 0 1 

*i*

1

( ) ( )

 *R*

 *X i*

*X i*

*L*

2

*D*

*tot*

*M*

( )

*i*

 *i*

When the *i*th thruster fails, the mathematical form of the rotational speed of the thruster can be expressed as follows: thruster can be expressed as follows: = *f* 

$$
\omega\_{\rm i}^f = \sigma\_{\rm l} \omega\_{\rm l} \tag{55}
$$

*i i i* (55)

*i*

(2 1)[ ( )] (2 1)[ ( )]

 *i*

 *i*

(2 1) ( )

 *i*

 *X i*

 *R*

> *R*

(54)

 *X i*

8

( )

8 1

*tot*

*T*

*X i*

*L*

*i*

8 1 1

*i*

*T<sup>L</sup>*eft and *T<sup>R</sup>*ight of the above equations were input as thrust commands to the electric

1

 *X i*

8 8

*R*

where 0 ≤ *σ<sup>i</sup>* < 1 denotes the fault rate of the *i*th thruster under a fault. When the motor is completely jammed, then *σ<sup>i</sup>* = 0. is completely jammed, then 0 *<sup>i</sup>* . When a simulation test of thruster fault-tolerant control is conducted in this simula-

When a simulation test of thruster fault-tolerant control is conducted in this simulation platform, a random fault thruster ID number, that is *nFault*, will be randomly generated in the *τ*th second in order to simulate a thruster fault more realistically. tion platform, a random fault thruster ID number, that is *Fault n* , will be randomly generated in the th second in order to simulate a thruster fault more realistically. When a thruster numbered *Fault n* fails in the th second, the torque of the thruster

When a thruster numbered *nFault* fails in the *τ*th second, the torque of the thruster corresponding to the failed thruster should be first controlled to the torque value of the failed one, i.e., *T<sup>i</sup>* ∗ *σ<sup>i</sup>* , at which point the difference between the thrust in a steady-state flight and that of a failed aircraft is deemed as the control error *ψerr*. In this case, the thruster control torque needed to recover the prefault thrust can be calculated by PID control, which will be fed back to the aircraft control input, in order to achieve the fault-tolerant control of thrusters in the DEP system. The functional block diagram of the designed fault-tolerant control methods for DEP aircrafts in this study is shown in Figure 7. corresponding to the failed thruster should be first controlled to the torque value of the failed one, i.e., *Ti i* , at which point the difference between the thrust in a steady-state flight and that of a failed aircraft is deemed as the control error *err* . In this case, the thruster control torque needed to recover the prefault thrust can be calculated by PID control, which will be fed back to the aircraft control input, in order to achieve the fault-tolerant control of thrusters in the DEP system. The functional block diagram of the designed fault-tolerant control methods for DEP aircrafts in this study is shown in Figure 7.

**Figure 7.** Functional block diagram of the fault-tolerant control of the DEP aircraft. **Figure 7.** Functional block diagram of the fault-tolerant control of the DEP aircraft.

#### **4. Simulation Results and Discussion**

*4.1. Simulation Tests Carried out within Mission Segments*

Simulation tests of the DEP aircraft on the coordinated and comprehensive control of thrust were conducted during the entire process in the mission profile, namely takeoff, cruise, and descent. The set flight conditions are shown in Table 1.


Simulation tests of the DEP aircraft on the coordinated and comprehensive control of thrust were conducted during the entire process in the mission profile, namely takeoff,

**Table 1.** Parameter setting in mission segments of takeoff/cruise/descent.

cruise, and descent. The set flight conditions are shown in Table 1.

**Table 1.** Parameter setting in mission segments of takeoff/cruise/descent.

*Drones* **2022**, *6*, x FOR PEER REVIEW 15 of 25

*4.1. Simulation Tests Carried out within Mission Segments*

**4. Simulation Results and Discussion**

The flight simulation test results of the control system within the full mission segments are presented in Figure 8, with the response curve of the flight altitude showing good tracking effects. good tracking effects.

ments are presented in Figure 8, with the response curve of the flight altitude showing

**Figure 8.** Response curve of the flight altitude. **Figure 8.** Response curve of the flight altitude.

and overshoot.

For the altitude control within the flight mission segment, the quantitative description of the control effect is shown in Table 2, including rise time, peak time, settling time, For the altitude control within the flight mission segment, the quantitative description of the control effect is shown in Table 2, including rise time, peak time, settling time, and overshoot.



Figure 9 is the acceleration response curve of the *z-*axis. As can be seen, there is a change of acceleration when the aircraft's flight state changes. The curve then converges to zero. Figure 10 is the velocity response curve of the *z-*axis. When the aircraft enters cruise from climb, changes in acceleration results in the aircraft's velocity in the *z-*axis Figure 9 is the acceleration response curve of the *z*-axis. As can be seen, there is a change of acceleration when the aircraft's flight state changes. The curve then converges to zero. Figure 10 is the velocity response curve of the *z*-axis. When the aircraft enters cruise from climb, changes in acceleration results in the aircraft's velocity in the *z*-axis reaching almost zero in order to maintain a flight state with constant height and uniform speed.

Overshoot 3.46 percent

Variation trend of the pitch angle of the aircraft in the corresponding mission segments is shown in Figure 11. During the cruise phase, the pitch angle of the aircraft returns to zero degrees.

speed.

speed.

*Drones* **2022**, *6*, x FOR PEER REVIEW 16 of 25

reaching almost zero in order to maintain a flight state with constant height and uniform

reaching almost zero in order to maintain a flight state with constant height and uniform

**Figure 9.** Acceleration response curve of the *z*-axis. **Figure 9.** Acceleration response curve of the *z*-axis. **Figure 9.** Acceleration response curve of the *z*-axis.

**Figure 10.** Velocity response curve of the *z*-axis. **Figure 10.** Velocity response curve of the *z*-axis.

**Figure 11.** Variation trend of the pitch angle in the mission segment. **Figure 11.** Variation trend of the pitch angle in the mission segment.

Figures 12 and 13 show the variation trend of the roll angle and yaw angle of an aircraft in the corresponding mission segments. The roll angle and the yaw angle will wit-

thrust and the attitude of thrusters in the DEP system. They will then return to a flight state without roll and deviation. The test results verifies the stability of the control system

**Figure 12.** Variation trend of the roll angle in the mission segment.

designed in this study.

Figures 12 and 13 show the variation trend of the roll angle and yaw angle of an aircraft in the corresponding mission segments. The roll angle and the yaw angle will witness some small changes at the moment the flight state switches due to changes in the thrust and the attitude of thrusters in the DEP system. They will then return to a flight state without roll and deviation. The test results verifies the stability of the control system designed in this study. Figures 12 and 13 show the variation trend of the roll angle and yaw angle of an aircraft in the corresponding mission segments. The roll angle and the yaw angle will witness some small changes at the moment the flight state switches due to changes in the thrust and the attitude of thrusters in the DEP system. They will then return to a flight state without roll and deviation. The test results verifies the stability of the control system designed in this study.

**Figure 11.** Variation trend of the pitch angle in the mission segment.

*Drones* **2022**, *6*, x FOR PEER REVIEW 17 of 25

**Figure 12.** Variation trend of the roll angle in the mission segment. **Figure 12.** Variation trend of the roll angle in the mission segment.

**Figure 13.** Variation trend of the yaw angle in the mission segment.

**Figure 13.** Variation trend of the yaw angle in the mission segment.

The thrust controller solves the control input torque of a corresponding single thruster to produce thrust in different stages. The curve of the total thrust variations generated by all thrusters in different mission segments is displayed in Figure 14. When the aircraft enters cruise in 300 s, the thrust required by the aircraft decreases, and the thrust is further reduced after it descends. The total power generated by the thruster module in The thrust controller solves the control input torque of a corresponding single thruster to produce thrust in different stages. The curve of the total thrust variations generated by all thrusters in different mission segments is displayed in Figure 14. When the aircraft enters cruise in 300 s, the thrust required by the aircraft decreases, and the thrust is further reduced after it descends. The total power generated by the thruster module in the entire mission segments is demonstrated in Figure 15.

#### *4.2. DEP System Thruster Fault-Tolerant Control Simulation Test*

the entire mission segments is demonstrated in Figure 15.

Propellers of thrusters on the left and right wings of an aircraft are designed to be righthanded and left-handed, and the torque direction is also symmetrical. When a thruster fails when *nFault* = 2 and *σ<sup>i</sup>* = 0.2 is randomly generated at the 200th second, the torque of the failed thruster instantly drops to the moment value of *T<sup>i</sup>* ∗ *σ<sup>i</sup>* , as shown by the red

**Figure 14.** Curve of the total thrust variations generated by the thruster module.

curve in Figure 16. Therefore, the torque of the symmetrical thruster No. 15 should change symmetrically in order to first ensure the balance of moment, as shown by the blue curve in Figure 16. aircraft enters cruise in 300 s, the thrust required by the aircraft decreases, and the thrust is further reduced after it descends. The total power generated by the thruster module in the entire mission segments is demonstrated in Figure 15.

The thrust controller solves the control input torque of a corresponding single thruster to produce thrust in different stages. The curve of the total thrust variations generated by all thrusters in different mission segments is displayed in Figure 14. When the

**Figure 13.** Variation trend of the yaw angle in the mission segment.

*Drones* **2022**, *6*, x FOR PEER REVIEW 18 of 25

**Figure 14.** Curve of the total thrust variations generated by the thruster module. **Figure 14.** Curve of the total thrust variations generated by the thruster module.

**Figure 15.** Curve of the total power variations generated by the thruster module. **Figure 15.** Curve of the total power variations generated by the thruster module.

*4.2. DEP System Thruster Fault-Tolerant Control Simulation Test* Propellers of thrusters on the left and right wings of an aircraft are designed to be right-handed and left-handed, and the torque direction is also symmetrical. When a thruster fails when 2 *Fault n* and 0.2 *<sup>i</sup>* is randomly generated at the 200th second, the torque of the failed thruster instantly drops to the moment value of *Ti i* , as shown by In order to recover a stable flight state, all thrusters, except the thrusters symmetrical to the failed ones, should increase their thrust, so the torque input of the rest thrusters should be up. As the remaining thrusters change in the same way, the response value of the thruster torque can be observed with thruster No. 1 as an example, and the input control torque of thruster No. 1 gradually increases after the fault occurs in the 200th second. The generated thrust also grows at the 200th second, as shown in Figure 17.

the red curve in Figure 16. Therefore, the torque of the symmetrical thruster No. 15 should change symmetrically in order to first ensure the balance of moment, as shown by the blue curve in Figure 16. The variation curve of the total thrust of the DEP aircraft after the fault is shown in Figure 18. As can be seen, the total thrust decreases after the thruster fault occurs in the 200th second, and the thrust of each thruster on the left and right is then altered by coordinate control to recover the thrust to a level that can maintain a stable flight of the aircraft.

**Figure 16.** Variation curve of moments of the faulty thruster and symmetrical thruster.

In order to recover a stable flight state, all thrusters, except the thrusters symmetrical to the failed ones, should increase their thrust, so the torque input of the rest thrusters should be up. As the remaining thrusters change in the same way, the response value of the thruster torque can be observed with thruster No. 1 as an example, and the input

**Figure 16.** Variation curve of moments of the faulty thruster and symmetrical thruster. **Figure 16.** Variation curve of torque of the faulty thruster and symmetrical thruster.

**Figure 15.** Curve of the total power variations generated by the thruster module.

torque of the failed thruster instantly drops to the moment value of *Ti i*

Propellers of thrusters on the left and right wings of an aircraft are designed to be right-handed and left-handed, and the torque direction is also symmetrical. When a

the red curve in Figure 16. Therefore, the torque of the symmetrical thruster No. 15 should change symmetrically in order to first ensure the balance of moment, as shown by the blue

*<sup>i</sup>* is randomly generated at the 200th second, the

, as shown by

*4.2. DEP System Thruster Fault-Tolerant Control Simulation Test*

thruster fails when 2 *Fault n* and 0.2

curve in Figure 16.

**Figure 17.** Variation curve of the remaining normal thrusters with thruster No. 1 as an example. dinate control to recover the thrust to a level that can maintain a stable flight of the aircraft.

The controller can better control the thrust to the prefault level without overshoot

**Figure 18.** Total thrust curve of the DEP aircraft after thruster fault. **Figure 18.** Total thrust curve of the DEP aircraft after thruster fault.

distributed by coordinated control after the fault.

1

The controller can better control the thrust to the prefault level without overshoot after the total thrust changes in the 200th second, taking 0.3 s for adjustment. In other words, thrust generated by the thrusters on the left and right wing of the aircraft is evenly distributed by coordinated control after the fault.

#### **5. Conclusions**

First, a mathematical model of the DEP aircraft's propulsion system, including the engine module, the generator and energy storage system module, and the thruster module, was established. Then, a mathematical model of the six-degree-of-freedom DEP aircraft was built based on the principles of aerodynamics and flight dynamics, which laid a theoretical foundation for subsequent simulation experiment.

Research on control methods to coordinate thrust from multiple thrusters were carried out based on the mathematical model of DEP aircrafts. The lateral and longitudinal control loops of DEP aircrafts were set up based on the principles of total energy and total heading control, and a simulation experiment was carried out in the mission segment of the DEP aircraft. The effects of aircraft attitude control and altitude control verified the stability and accuracy of the mathematical model of the aircraft.

Furthermore, a fault-tolerant control method was developed for the case where a thruster of a DEP aircraft has failed. Experiments simulating flight tests and fault-tolerant control within the mission segment were conducted, and the experimental results verified the effectiveness of the designed coordinated thrust control system and the fault-tolerant control method. The controller could control the thrust to the prefault level.

The correctness and effectiveness of the designed coordinated thrust control method and fault-tolerant control method for DEP aircrafts were theoretically verified, providing a theoretical basis for future engineering application and development of the control system for DEP aircrafts.

**Author Contributions:** Conceptualization, J.L.; methodology, J.L.; validation, J.L.; formal analysis, J.L.; investigation, J.L.; data curation, J.L.; writing—original draft preparation, J.L.; writing—review and editing, J.Y.; supervision, H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **Nomenclature**





## **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Drones* Editorial Office E-mail: drones@mdpi.com www.mdpi.com/journal/drones

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-7560-5