2.1 Forward flight aerodynamic model: wing-body aerodynamics at high speeds

The USAF Digital DATCOM is a computer program that estimates aerodynamic coefficients using semi-empiric methods [16, Reference Rosema, Doyle, Auman, Underwood and Blake18]. DATCOM has been used in literature and industry for a long time to build preliminary flight dynamics models of several aircraft and missiles [Reference Roskam19–Reference Vasile, Bryson and Fresconi23]. It gives very fast and accurate enough results for the initial design phase [Reference Rosema, Doyle, Auman, Underwood and Blake18–Reference Roskam20]. Detailed verification of DATCOM methods are also studied in the literature [Reference Sooy and Schmidt21, Reference Segui, Mantilla and Botez24, Reference Doyle and Rosema25]. To conclude, Digital DATCOM program is considered as an effective way to estimate the wing-body aerodynamics with acceptable accuracy since the main focus of the study is not detailed aerodynamic modeling. At this point, it is good to emphasise that detailed aerodynamic modeling of this unique aircraft is still an open area for researchers.

The input file of DATCOM for our case mainly contains 3-D geometry and desired flight conditions. The wing-body geometry of the aircraft is constructed iteratively to satisfy the following stability requirements [Reference McLean26] and aerodynamic efficiency at cruise represented as the lift-to-drag ratio ( ${C}_{L}/{C}_{D}$ ).

1. 1. ${C}_{m\alpha }\lt 0$ $\equiv $ longitudinal static stability

2. 2. ${C}_{n\beta }\gt 0$ $\equiv $ static directional(weathercock) stability (could not be satisfied due to the tailless design)

3. 3. ${C}_{l\beta }\lt 0$ $\equiv $ lateral static stability

4. 4. ${C}_{Y\beta }\lt 0$ $\equiv $ sideslip stability

5. 5. ${C}_{L}/{C}_{D}\gt 5$ for 0 < ${\alpha }_{cruise}\lt 10deg$ , ${\alpha }_{cruise}$ is the angle-of-attack at cruise

Directional stability at forward flight could not be satisfied ( ${C}_{n\beta }$ is negative, see Fig. 3). This is an expected result since the aircraft has no vertical tail [Reference McLean26]. As mentioned previously, not having conventional control and stability surfaces comes with several advantages. On the other hand, there are some drawbacks regarding the flight control such as open loop unstability in the directional channel, limited control authority and problems introduced by the couplings due to pure thrust vector control (i.e. couplings between control axis and control effectors).

Figure 3. Non-dimensional coefficients of the forward flight aerodynamic model obtained via Digital DATCOM.

Flight conditions are chosen according to the expected operation of the air taxi (see Table 1). Maximum airspeed is considered as 0.3 Mach ( $\approx $ 100 m/s). Operations are assumed to be performed at low altitudes (0–1,000 m). It is observed that aerodynamic coefficients don’t change significantly with respect to altitude for the considered altitude range. Therefore, sea level altitude is used in the DATCOM input file for simplicity. Digital DATCOM gives results at discrete flight conditions, which are represented as functions via curve-fitting (Fig. 3). The forward-flight model is dominant after 20 m/s airspeed. Therefore, the validity ranges of the functions are $-15\mathrm{}\mathrm{}\mathrm{d}\mathrm{e}\mathrm{g}\lt \alpha \lt +15deg,-15deg\lt \beta \lt +15,0.05\lt \mathrm{M}\lt 0.5$ . Details of the DATCOM modeling, input file and non-dimensional coefficients (output of DATCOM) can be found in [Reference Suicmez27].

Table 1. General parameters of the air taxi

Overall static and dynamic contributions of the coefficients are given in Equation (1).

(1) \begin{align}\begin{array}{r}\left.\begin{array}{r}{C}_{{D}_{,sta}}={C}_{D}\\[5pt] {C}_{{y}_{,sta}}={C}_{y\beta }\cdot \beta \cdot 180/\pi \\[5pt] {C}_{{L}_{,sta}}={C}_{L}\\[5pt] {C}_{{l}_{,sta}}={C}_{l\beta }\cdot \beta \cdot 180/\pi \\[5pt] {C}_{{m}_{,sta}}={C}_{m}\\[5pt] {C}_{{n}_{,sta}}={C}_{n\beta }\cdot \beta \cdot 180/\pi \end{array}\right\}\end{array}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{l}\left.\begin{array}{r}{C}_{{D}_{,dyn}}=0\\[5pt] {C}_{{y}_{,dyn}}={C}_{yp}\cdot p\cdot 180/\pi \cdot ({b}_{ref}/(2{V}_{\mathrm{\infty }}\left)\right)\\[5pt] {C}_{{L}_{,dyn}}=({C}_{Lq}\cdot q+{C}_{L\dot{\alpha }}\cdot \dot{\alpha })\cdot 180/\pi \cdot (\stackrel{-}{c}/(2{V}_{\mathrm{\infty }}\left)\right)\\[5pt] {C}_{{l}_{,dyn}}=({C}_{lp}\cdot p+{C}_{lr}\cdot r)\cdot 180/\pi \cdot ({b}_{ref}/(2{V}_{\mathrm{\infty }}\left)\right)\\[5pt] {C}_{{m}_{,dy}}=({C}_{mq}\cdot q+{C}_{m\dot{\alpha }}\cdot \dot{\alpha })\cdot 180/\pi \cdot (\stackrel{-}{c}/(2{V}_{\mathrm{\infty }}\left)\right)\\[5pt] {C}_{{n}_{,dyn}}=({C}_{np}\cdot p+{C}_{nr}\cdot r)\cdot 180/\pi \cdot ({b}_{ref}/(2{V}_{\mathrm{\infty }}\left)\right)\end{array}\right\}\end{array}\end{align}

Subscripts $sta$ and $dyn$ represent the static and dynamic contributions. ${C}_{D},{C}_{y},{C}_{L}$ are the non-dimensional drag, side-force and lift coefficients; ${C}_{l},{C}_{m},{C}_{n}$ are the non-dimensional moment coefficients in 3-axis; ${V}_{\mathrm{\infty }}$ is the freestream velocity; $\alpha $ and $\beta $ represent the angle-of-attack and sideslip angle; $p,q,r$ are the body rotational rates; and ${b}_{ref}$ and $\stackrel{-}{c}$ are the reference lengths defined in Table 1.

Equation (2) gives the overall forward-flight aerodynamic force and moment expressed in the body frame. It is noted that the DATCOM gives the coefficients in wind axes; however, the transformation from wind axes to body axes is considered insignificant due to the small angle-of-attack and sideslip angles.

(2) \begin{align} {F}_{aero,\;forward}^{b}=\left[\begin{array}{c}-\stackrel{-}{\mathrm{q}}\mathrm{S}{\mathrm{C}}_{{\mathrm{D}}_{,\mathrm{s}\mathrm{t}\mathrm{a}}}\\[5pt] \stackrel{-}{\mathrm{q}}\mathrm{S}({\mathrm{C}}_{{\mathrm{y}}_{,\mathrm{s}\mathrm{t}\mathrm{a}}}+{\mathrm{C}}_{{\mathrm{y}}_{,\mathrm{d}\mathrm{y}\mathrm{n}}})\\[5pt] -\stackrel{-}{\mathrm{q}}\mathrm{S}({\mathrm{C}}_{{\mathrm{L}}_{,\mathrm{s}\mathrm{t}\mathrm{a}}}+{\mathrm{C}}_{{\mathrm{L}}_{,\mathrm{d}\mathrm{y}\mathrm{n}}})\\[5pt] \end{array}\right]\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad {M}_{aero,forward}^{b}=\left[\begin{array}{c}\stackrel{-}{\mathrm{q}}\mathrm{S}({\mathrm{C}}_{{\mathrm{l}}_{,\mathrm{s}\mathrm{t}\mathrm{a}}}+{\mathrm{C}}_{{\mathrm{l}}_{,\mathrm{d}\mathrm{y}\mathrm{n}}}){\mathrm{b}}_{\mathrm{r}\mathrm{e}\mathrm{f}}\\[5pt] \stackrel{-}{\mathrm{q}}\mathrm{S}({\mathrm{C}}_{{\mathrm{m}}_{,\mathrm{s}\mathrm{t}\mathrm{a}}}+{\mathrm{C}}_{{\mathrm{m}}_{,\mathrm{d}\mathrm{y}\mathrm{n}}})\stackrel{-}{\mathrm{c}}\\[5pt] \stackrel{-}{\mathrm{q}}\mathrm{S}({\mathrm{C}}_{{\mathrm{n}}_{,\mathrm{s}\mathrm{t}\mathrm{a}}}+{\mathrm{C}}_{{\mathrm{n}}_{,\mathrm{d}\mathrm{y}\mathrm{n}}}){\mathrm{b}}_{\mathrm{r}\mathrm{e}\mathrm{f}}\\[5pt] \end{array}\right]\end{align}

$\stackrel{-}{q}=1/2\rho {V}_{\mathrm{\infty }}^{2}$ is the dynamic pressure with $\rho $ being the air-density. $S$ is the reference area defined in Table 1.

2.2 Hover flight aerodynamic model: resistance force at low speeds around hover

During the low-speed flight around hover, the most dominant aerodynamic effect is the resistance drag force since aerodynamic lift and moments are negligible [Reference Lombaerts, Kaneshige, Schuet, Hardy, Aponso and Shish28]. Ground effect is also relevant but it is not modeled for the scope of this study. The low-speed flight around hover is restricted to 10 m/s for horizontal flight and 5 m/s for vertical flight.

In a study about modeling and flight control of a fixed wing VTOL aircraft, the preliminary hover model is based on flat plate aerodynamic modeling [Reference Lombaerts, Kaneshige, Schuet, Hardy, Aponso and Shish28]. A similar approach is used to generate a preliminary aerodynamic model at hover. Drag resistance is roughly estimated for each direction based on a rectangular flat-plate model for the y and z directions, and a circular cylinder model for the x-direction. The frontal areas used in the estimations are calculated as follows. ${l}_{fus}$ , ${h}_{fus}$ , and $S$ are defined in Table 1.

(3) \begin{align}{S}_{x}=\pi \mathrm{}\mathrm{}{h}_{fus}^{2}/4\approx 3\mathrm{}\mathrm{}{m}^{2},\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{S}_{y}={l}_{fus}\mathrm{}\mathrm{}{h}_{fus}=8\mathrm{}\mathrm{}{m}^{2},\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{S}_{z}={l}_{fus}\mathrm{}\mathrm{}{h}_{fus}+S\approx 10\mathrm{}\mathrm{}{m}^{2}\end{align}

Based on the frontal areas, aerodynamic drag resistance force around hover is expressed in body coordinates as in Equation (4). Aircraft can move forward/backward and sideward around the hover so that the sign function is added to include the direction of motion into the equations. Body velocities ( $u,v,w$ ) are used in the equations assuming that roll and pitch angles are small for the low-speed flight around hover. The drag coefficients for the y and z directions ( ${C}_{d,y}$ and ${C}_{d,z}$ ) are equal to 1.2 for a rectangular plate, and the drag coefficient for the x-direction ( ${C}_{d,x}$ ) is equal to 0.74 for a circular cylinder [Reference Prasuhn29].

(4) \begin{align}{F}_{aero,hover}^{b}=\left[\begin{array}{l}-sign\left(u\right)\mathrm{}\mathrm{}0.5\mathrm{}\mathrm{}\rho \mathrm{}\mathrm{}{u}^{2}{S}_{x}{C}_{d,x}\\[5pt] -sign\left(v\right)\mathrm{}\mathrm{}0.5\mathrm{}\mathrm{}\rho \mathrm{}\mathrm{}{v}^{2}{S}_{y}{C}_{d,y}\\[5pt] -sign\left(w\right)\mathrm{}\mathrm{}0.5\mathrm{}\mathrm{}\rho \mathrm{}\mathrm{}{w}^{2}{S}_{z}{C}_{d,z}\\[5pt] \end{array}\right]\end{align}

2.3 Complete aerodynamic model

For high- and low-speed regions, aerodynamic models are generated separately in previous sections. These models are merged to express the complete aerodynamic model considering the transition. The merging/decision parameter is chosen as body x velocity ( $u$ ). The airspeed ( ${V}_{\mathrm{i}\mathrm{n}\mathrm{f}}$ ) is not used as a decision parameter since it increases during the sideward motion, but the wing does not generate lift.

Based on the DATCOM results, the wing-body aerodynamics seems to be dominant after 20 m/s (0.05 Mach) so that the forward flight model is fully active after $u=20\mathrm{}\mathrm{m}/\mathrm{s}$ . Transition from hover to forward aerodynamic model starts at $u=10\mathrm{}\mathrm{m}/\mathrm{s}$ .

The complete aerodynamic model is given as a combination of the hover and forward aerodynamic models in Equation (5) via introducing the merging coefficient ${k}_{aero}$ .

(5) \begin{align}\begin{array}{r}\left.\begin{array}{r}{F}_{aero}^{b}=(1-{k}_{aero})\mathrm{}\mathrm{}{F}_{aero,forward}^{b}+{k}_{aero}\mathrm{}\mathrm{}{F}_{aero,hover}^{b}\\[5pt] {M}_{aero}^{b}=(1-{k}_{aero})\mathrm{}\mathrm{}{M}_{aero,forward}^{b}\end{array}\right\}\end{array}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad {k}_{aero}=\left\{\begin{array}{l@{\quad}l}1 & \mathrm{u}\lt 10\\[5pt] (20-u)/10 & 10\le \mathrm{u}\le 20\\[5pt] 0 & \mathrm{u}\gt 20\end{array}\right.\end{align}

2.4 Propulsion system modelling

The propulsion system consists of several EDFs distributed over the wing and front sections (Fig. 2). The datasheet of a commercially available EDF named Schubeler DS-215-DIA HST is used [17] to estimate the thrust and torque generated by each EDF. A similar EDF is also used in a NASA founded DEP aircraft project [Reference Freeman30]. The basic parameters of the EDF are given in Table 2.

Table 2. Parameters of the Schubeler DS-215-DIA HST EDF

Based on the datasheet, thrust ( $T$ ) and torque ( $Q$ ) of each EDF are modeled as follows. Values of ${C}_{T}$ and ${C}_{Q}$ are defined in Table 1.

(6) \begin{align}T={C}_{T}\mathrm{}\mathrm{}(rpm\mathrm{}\mathrm{}2\pi /60{)}^{2}\mathrm{}\mathrm{}\mathrm{}\mathrm{},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad Q={C}_{Q}\mathrm{}\mathrm{}T\end{align}

Derivation of the total force and moment generated by the propulsion systems depends on the distribution of EDFs over the wing and front sections (Fig. 2) and the thrust vector control concept (Fig. 4). Tilt angle ( $\delta $ ) and thrust ( $T$ ) of EDFs are adjusted to control the air taxi. Since rpm is related to thrust with a constant coefficient (i.e. ${C}_{T}$ ), thrust is used as a control parameter instead of rpm for simplicity. There are 26 EDFs with 18 distributed over the wings and 8 on the front. Distribution of 26 EDFs between the front and wing sections is based on the pitch moment balance at hover (see EDF lever arms info in Table 3, and top view in Fig. 2). There are six EDF sets (each has three EDF) on the wings and four EDF sets (each has two EDF) on the front sections (Fig. 2). EDF sets are combined for left/right front and wing sections to simplify the controller design. By this way number of control inputs is reduced to eight (i.e. four tilt angle ( ${\delta }_{wl},{\delta }_{wr},{\delta }_{fl},{\delta }_{fr}$ ) and four thrust ( ${T}_{wl},{T}_{wr},{T}_{fl},{T}_{fr}$ ) for each EDF section). It is also noted that the average lever arms defined in Table 3 are used in the controller. Lever arm information is a direct input to the controller so that it is possible to manipulate this input properly in case of failures or any other required corrections. Table 3 gives the detailed EDF parameters for the left/right wing and front sections.

Figure 4. Thrust vector control concept, side view.

Based on the thrust vector concept illustrated in Fig. 4, total force and moment generated at each section (front left, front right, wing left and, wing right) in the body coordinate system are represented in a general equation form as follows. The parameters used in the following equations are given in Table 3.

(7) \begin{align}{F}_{i}^{b}=\left[\begin{array}{c@{\quad}c@{\quad}c}\mathrm{c}\mathrm{o}\mathrm{s}\left({\delta }_{\mathrm{i}}\right) &0& \mathrm{s}\mathrm{i}\mathrm{n}\left({\delta }_{\mathrm{i}}\right)\\[5pt] 0 & 1& 0\\[5pt] -\mathrm{s}\mathrm{i}\mathrm{n}\left({\delta }_{\mathrm{i}}\right) &0& \mathrm{c}\mathrm{o}\mathrm{s}\left({\delta }_{\mathrm{i}}\right)\\[5pt] \end{array}\right]\left[\begin{array}{c}{\mathrm{T}}_{\mathrm{i}}\\[5pt] 0\\[5pt] 0\\[5pt] \end{array}\right]\mathrm{}\mathrm{}{n}_{i}=\left[\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\left({\delta }_{\mathrm{i}}\right)\mathrm{}\mathrm{}{\mathrm{T}}_{\mathrm{i}}\\[5pt] 0\\[5pt] -\mathrm{s}\mathrm{i}\mathrm{n}\left({\delta }_{\mathrm{i}}\right)\mathrm{}\mathrm{}{\mathrm{T}}_{\mathrm{i}}\\[5pt] \end{array}\right]\mathrm{}\mathrm{}{n}_{i}\mathrm{}\mathrm{},\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}with\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad i=fl,fr,wl,wr\end{align}

(8) \begin{align}\begin{array}{c}{M}_{i}^{b}=\left[\begin{array}{c@{\quad}c@{\quad}c}\mathrm{c}\mathrm{o}\mathrm{s}\left({\delta }_{\mathrm{i}}\right) &0& \mathrm{s}\mathrm{i}\mathrm{n}\left({\delta }_{\mathrm{i}}\right)\\[5pt] 0 &1& 0\\[5pt] -\mathrm{s}\mathrm{i}\mathrm{n}\left({\delta }_{\mathrm{i}}\right) &0& \mathrm{c}\mathrm{o}\mathrm{s}\left({\delta }_{\mathrm{i}}\right)\\[5pt] \end{array}\right]\left[\begin{array}{c}{\mathrm{Q}}_{\mathrm{i}}\\[5pt] 0\\[5pt] 0\\[5pt] \end{array}\right]\mathrm{}\mathrm{}t{d}_{i}\mathrm{}\mathrm{}{n}_{i}+\left[\begin{array}{c}\mathrm{\Delta }{\mathrm{x}}_{\mathrm{i}}\\[5pt] \mathrm{\Delta }{\mathrm{y}}_{\mathrm{i}}\\[5pt] \mathrm{\Delta }{\mathrm{z}}_{\mathrm{i}}\\[5pt] \end{array}\right]\times \left[\begin{array}{c}\mathrm{c}\mathrm{o}\mathrm{s}\left({\delta }_{\mathrm{i}}\right)\mathrm{}\mathrm{}{\mathrm{T}}_{\mathrm{i}}\\[5pt] 0\\[5pt] -\mathrm{s}\mathrm{i}\mathrm{n}\left({\delta }_{\mathrm{i}}\right)\mathrm{}\mathrm{}{\mathrm{T}}_{\mathrm{i}}\\[5pt] \end{array}\right]\mathrm{}\mathrm{}{n}_{i}\mathrm{}\mathrm{},\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}with\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad i=fl,fr,wl,wr\end{array}\end{align}

Then, the overall propulsion force and moment in the body coordinate system ( ${F}_{prop}^{b}$ and ${M}_{prop}^{b}$ ) are the sum of each section.

(9) \begin{align}\begin{array}{r}\left.\begin{array}{r}{F}_{prop}^{b}=[F{x}_{prop}\quad 0\quad F{z}_{prop}{]}^{T}={F}_{fl}^{b}+{F}_{fr}^{b}+{F}_{wl}^{b}+{F}_{wr}^{b}\\[5pt] {M}_{prop}^{b}=[{L}_{prop}\quad {M}_{prop}\quad {N}_{prop}{]}^{T}={M}_{fl}^{b}+{M}_{fr}^{b}+{M}_{wl}^{b}+{M}_{wr}^{b}\end{array}\right\}\end{array}\end{align}

2.5 Actuator dynamics

There are two control inputs for each EDF, which are the thrust ( $T$ ) and the tilt angle ( $\delta $ ). Thrust dynamics of a similar Schubeler EDF considered in this study is modeled in a work based on the wind-tunnel data as a second order system with a natural frequency of 18.85 rad/s and damping ratio of 1 [Reference Freeman and Klunk31]. Moreover, the maximum thrust of a Schubeler EDF is given as 250 N in the datasheet (Table 2). Based on these reference values, a little improved EDF system with a maximum thrust of 300 N ( $\approx $ 15,000 $\mathrm{r}\mathrm{p}\mathrm{m}$ ) and a natural frequency of 25 rad/s is used in this study. The overall parameters of the thrust dynamics are given in Table 4, considering the number of EDFs in the front and wing sections (Table 3).

Table 3. Parameters of the EDF sections

Table 4. Parameters of actuator dynamics

There is no data to be used for the tilting dynamics of EDFs. Considering the small size/mass of Schuler EDF (Table 2), the required hinge-moments to rotate the EDF sets are relatively small. For this reason, it is reasonable to assume quite fast tilting dynamics. Regarding the tilt angle limits, the wing EDF sets are constrained to 0° tilt angle, which is a limitation applied due to the wing structure. For EDFs at the front section there is no fixed structure, so that the minimum tilt angle is constrained to −30°. For both the front and wing EDFs, the maximum tilt angle is chosen as 120°. The overall parameters of the tilt angle actuator dynamics are also given in Table 4.

2.6 Six DOF equations of motion

Six degree of freedom (6-DOF) equations of motion are obtained by using the overall aerodynamic and propulsion forces/moments defined in Equations (5) and (9), respectively.

(10) \begin{align}\begin{array}{c}\left[\begin{array}{c}\dot{\mathrm{u}}\\[5pt] \dot{\mathrm{v}}\\[5pt] \dot{\mathrm{w}}\\[5pt] \end{array}\right]=({F}_{aero}^{b}+{F}_{prop}^{b}+{F}_{grav}^{b}+{F}_{dist}^{b})/m-\left[\begin{array}{c}\mathrm{p}\\[5pt] \mathrm{q}\\[5pt] \mathrm{r}\\[5pt] \end{array}\right]\times \left[\begin{array}{c}\mathrm{u}\\[5pt] \mathrm{v}\\[5pt] \mathrm{w}\\[5pt] \end{array}\right],\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad with\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}{F}_{grav}^{b}=\left[\begin{array}{c}-\mathrm{g}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\\[5pt] \mathrm{g}\mathrm{s}\mathrm{i}\mathrm{n}\left(\phi \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\\[5pt] \mathrm{g}\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi \right)\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\\[5pt] \end{array}\right]\end{array}\end{align}

(11) \begin{align}\begin{array}{l}\left[\begin{array}{l}\dot{\mathrm{p}}\\[5pt] \dot{\mathrm{q}}\\[5pt] \dot{\mathrm{r}}\\[5pt] \end{array}\right]={J}^{-1}\left({M}_{aero}^{b}+{M}_{prop}^{b}+{M}_{dist}^{b}-\left[\begin{array}{l}\mathrm{p}\\[5pt] \mathrm{q}\\[5pt] \mathrm{r}\\[5pt] \end{array}\right]\times J\left[\begin{array}{l}\mathrm{p}\\[5pt] \mathrm{q}\\[5pt] \mathrm{r}\\[5pt] \end{array}\right]\right)\end{array}\end{align}

${F}_{grav}^{b}$ represents the gravitational force in the body coordinates, with $\phi $ and $\theta $ being the roll and pitch angles. $m$ and $J$ are the mass and inertia matrix defined in Table 1. ${F}_{dist}^{b}$ and ${M}_{dist}^{b}$ are the disturbance force and moment introduced in the simulations to test the disturbance rejection characteristics of the designed controller.

3.1 Control input definitions

Before diving into the formulation of the INDI control law, control inputs are defined first. As shown in Fig. 4, the thrust vector is adjusted to control both translational and rotational dynamics. Tilt angle and rpm are the physical actuator states, but for simplification, thrust is used instead of rpm based on the relation given in Equation (6). Using this simplification, the control input vector on the physical actuator level includes the tilt angles( $\delta $ ) and thrust magnitudes( $T$ ). ${U}_{act}$ and $U$ represent the physical level control inputs with and without considering the actuator dynamics, respectively. ${U}_{act}$ and $U$ are separately defined to check the CA performance (see Section 4 for details). On the aircraft level, our aim is to generate desired forces and moments via controlling the thrust vector. Therefore, a virtual control input vector ( ${v}_{INDI}$ ) that includes the forces and moments generated by the propulsion system on body axes, defined in Equation (9), is introduced. Moreover, to get rid of the nonlinearities of the thrust vector concept, the INDI control input vector ( ${U}_{INDI}$ ) is defined considering the thrust contribution on the body x and z directions (Fig. 4). In this way, it is possible to relate the virtual and INDI control inputs using only the lever arm information, data input to the controller, and calculated using simple geometry. Equation (12) gives the definitions and relations of the different control inputs. The naming convention of the control inputs and other data such as the mean lever arms are given in Table 3.

(12) \begin{align}\begin{array}{l} U={\left[\begin{array}{l}{T}_{fl}\quad {T}_{fr}\quad {T}_{wl}\quad {T}_{wr}\quad {\delta }_{fl}\quad {\delta }_{fr}\quad {\delta }_{wl}\quad {\delta }_{wr}\end{array}\right]}^{T}\nonumber \\[5pt] {U}_{act}={\left[\begin{array}{l}{T}_{fl,act}\quad {T}_{fr,act}\quad {T}_{wl,act}\quad {T}_{wr,act}\quad {\delta }_{fl,act}\quad {\delta }_{fr,act}\quad {\delta }_{wl,act}\quad {\delta }_{wr,act}\end{array}\right]}^{T}\nonumber\\[5pt] {T}_{x}=sin(\pi /2-\delta )T,\quad {T}_{z}=cos(\pi /2-\delta )T,\quad T=\sqrt{{T}_{x}^{2}+{T}_{z}^{2}}\nonumber\\[5pt] {v}_{INDI}={\left[\begin{array}{l}{L}_{prop}\quad {M}_{prop}\quad {N}_{prop}\quad F{z}_{prop}\quad F{x}_{prop}\end{array}\right]}^{T}\nonumber\\[5pt] {U}_{INDI}={\left[\begin{array}{l}{T}_{x,fl}\quad {T}_{x,fr}\quad {T}_{x,wl}\quad {T}_{x,wr}\quad {T}_{z,fl}\quad {T}_{z,fr}\quad {T}_{z,wl}\quad {T}_{z,wr}\end{array}\right]}^{T}\nonumber\\[5pt] {v}_{INDI}={T}_{v,INDI}{U}_{INDI},\quad {U}_{INDI}={T}_{v,INDI}^{-1}{v}_{INDI},\quad {T}_{v,INDI}^{-1}={T}_{v,INDI}^{T}({T}_{v,INDI}{T}_{v,INDI}^{T}{)}^{-1}\nonumber\\[5pt]\end{array}\\ \left[\begin{array}{c}{L}_{prop}\\[5pt] {M}_{prop}\\[5pt] {N}_{prop}\\[5pt] F{z}_{prop}\\[5pt] F{x}_{prop}\\[5pt] \end{array}\right]=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0 & 0& 0& 0 & \Delta {y}_{fl} & -\Delta {y}_{fr} & \Delta {y}_{wl} & -\Delta {y}_{wr} \\[5pt] 0 & 0 & 0 & 0 & \Delta {x}_{fl} & \Delta {x}_{fr} & -\Delta {x}_{wl} & -\Delta {x}_{wr} \\[5pt] \Delta {y}_{fl} & -\Delta {y}_{fr} & \Delta {y}_{wl} & -\Delta {y}_{wr} & 0 & 0 & 0 & 0 \\[5pt] 0 &0 & 0 & 0 & -1 & -1 & -1 & -1 \\[5pt] 1 &1 & 1& 1& 0& 0& 0& 0 \end{array}\right]\left[\begin{array}{l}{T}_{x,fl}\\[5pt] {T}_{x,fr}\\[5pt] {T}_{x,wl}\\[5pt] {T}_{x,wr}\\[5pt] {T}_{z,fl}\\[5pt] {T}_{z,fr}\\[5pt] {T}_{z,wl}\\[5pt] {T}_{z,wr} \end{array}\right] \end{align}

As formulated in Equation (12), system is over-actuated since there are eight control effectors to control five axis/channels ( $L,M,N,Fz,Fx$ and noted that side-force $Fy$ is indirectly controlled via roll angle ( $\phi $ )-roll moment( $L$ ), see Section 3.2.1 for details). Pseudo-inverse approach is used to find solution for the over-actuated system. However, pseudo-inverse approach is not sufficient to find a solution that allocates the limited control authority properly to guarantee stable flight in case of actuator saturation. Therefore, additional consideration is needed regarding the CA design to take into account actuator saturation. Details of the CA design and its criticality are discussed in Section 3.3.

3.2 Formulation of the INDI control law

State space representation of the EOM defined in Equations (10) and (11) can be rewritten in a control-affine form.

(13) \begin{align}\begin{array}{l}\dot{x}=f\left(x\right)+g\quad {v}_{INDI},\quad x={\left[\begin{array}{l}\mathrm{p}\quad \mathrm{q}\quad\mathrm{r}\quad\mathrm{w}\quad \mathrm{u}\end{array}\right]}^{T},\quad g=diag\!\left(1/{I}_{x},\quad 1/{I}_{y},\quad 1/{I}_{z},\quad 1/m,\quad 1/m\right)\end{array}\end{align}

In Equation (13), $f\left(x\right)$ includes all the components other than the propulsion forces and moments that are defined as the virtual control input ( ${v}_{INDI}$ ) in Equation (12). State vector ( $x$ ) does not contain the body y-velocity ( $v$ ) since the propulsion system can not directly generate forces on the body y-direction. As explained later, the roll and yaw angle outer loop controllers are used to indirectly control the motion on body y-direction for low and high speed flight.

Once the control-affine state space representation is defined, it is possible to write the first-order Taylor series approximation evaluated at the previous time step (subscript “0") of the states and the virtual control input.

(14) \begin{align}\begin{array}{l}\dot{x}\approx {\dot{x}}_{0}+{\left.\frac{\partial \left[f\right(x)+g\mathrm{}\mathrm{}{v}_{INDI}]}{\partial x}\right|}_{\substack{x={x}_{0} \\ {v}_{INDI}={v}_{INDI,0}}}(x-{x}_{0})\mathrm{}\mathrm{}+{\left.\frac{\partial \left[f\right(x)+g\mathrm{}\mathrm{}{v}_{INDI}]}{\partial {v}_{INDI}}\right|}_{\substack{x={x}_{0}\\ {v}_{INDI}={v}_{INDI,0}}}({v}_{INDI}-{v}_{INDI,0})\end{array}\end{align}

The primary assumption of the INDI approach is the time-scale separation principle. Suppose the actuator dynamics are relatively fast, and the time step is small enough. In that case, it is possible to assume that the changes in states are relatively small compared to the changes in control inputs during a single time step [Reference Smeur, Chu and de Croon35, Reference Pollack, Looye and Van der Linden39]. Considering the fast thrust/rpm dynamics of EDFs and 100 Hz sample time, time scale separation fits well with the problem. Regarding the tilting dynamics of EDFs (see Table 4), it is much slower than the thrust/rpm dynamics, so that time scale separation principle may not be satisfied. This was one of the open questions at the beginning of the study. However, simulation results show that the proposed tilt angle dynamics is also fast enough to satisfy the time scale separation principle. Then, state term is assumed negligible and removed from the Taylor series expansion based on the time-scale separation principle.

(15) \begin{align}\begin{array}{l}\dot{x}\approx {\dot{x}}_{0}+g\mathrm{}\mathrm{}({v}_{INDI}-{v}_{INDI,0})={\dot{x}}_{0}+g\mathrm{}\mathrm{}\mathrm{\Delta }{v}_{INDI}\end{array}\end{align}

Inverting Equation (15) and also using the relations given in Equation (12), the incremental form of the virtual and INDI control inputs ( $\mathrm{\Delta }{v}_{INDI}$ and $\mathrm{\Delta }{U}_{INDI}$ ) are found.

(16) \begin{align} \begin{array}{l}\mathrm{\Delta }{v}_{INDI}={g}^{-1}(\dot{x}-{\dot{x}}_{0})\mathrm{}\mathrm{}\mathrm{}\mathrm{},\mathrm{}\mathrm{}\;\;where\;\;\mathrm{}\mathrm{}\mathrm{}\mathrm{}\dot{x}={\dot{x}}_{req}\\[5pt] \mathrm{\Delta }{U}_{INDI}={T}_{v,INDI}^{-1}\mathrm{\Delta }{v}_{INDI}\end{array}\end{align}

The incremental INDI control input depends on the following terms:

• ${\dot{x}}_{0}$ : The state derivative vector evaluated using the current and previous values of sensor outputs. A simple 2nd order filter is used to filter the noise due to the derivative. See Section 3.2.2.

• ${\dot{x}}_{req}$ : The required state derivative vector to track the attitude and velocity commands. A simple linear controller based on kinematic relations is defined in Section 3.2.1 to generate ${\dot{x}}_{req}$ .

• ${g}^{-1}$ : The inverse of $g$ matrix, which depends only on the inertia and mass. See Equation (13).

• ${T}_{v,INDI}^{-1}$ : Pseudo inverse of ${T}_{v,INDI}$ matrix, which depends only on the mean lever arms. See Equation (12).

The mean lever arms are calculated based on the geometry (Fig. 2), and it is assumed that this information is an input to the controller similar to the mass and inertia. If actuator failures occur in some EDFs, then the mean lever arm change accordingly. To sum up, the incremental INDI control input only depends on the mass/inertia and the lever arms (geometry) of the aircraft. The flight dynamics model information required for the dynamic inversion is replaced with the state derivative estimation based on the sensor measurements. This is a significant advantage regarding the improved robustness, especially designing a controller for a novel aircraft that is hard to model accurately. However, INDI becomes sensitive to the noise and delay introduced by the state derivative estimation using the sensor measurements. Filtering is required to reduce the noise level and to improve the INDI performance. It is crucial to test the INDI controller with a realistic sensor model due to the sensitivity. More information will be given about the sensor model and noise filtering in Section 3.2.2.

The overall INDI control input is represented as following, again ${U}_{INDI,0}$ represents the control input at the previous time step. It is important to note that the filter applied to estimate the state derivatives is also applied to the INDI control input to satisfy the data synchronisation (see Fig. 5).

(17) \begin{align}\begin{array}{l}{U}_{INDI}={U}_{INDI,0}+\mathrm{\Delta }{U}_{INDI}\end{array}\end{align}

Figure 5. High-level block diagram of the simulation model.

3.2.1 Linear controller: obtaining the required state derivatives

The aim is to control the aircraft’s attitude and velocity, and it is possible to generate the required state derivatives using a proportional-derivative (PD) type linear controller as follows. The linear controller gains are defined in Table 5.

(18) \begin{align} \begin{array}{r}\left.\begin{array}{r}{\dot{x}}_{req}={\left[\begin{array}{r}{\dot{\mathit{p}}}_{\mathit{r}\mathit{e}\mathit{q}}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\quad {\dot{\mathit{q}}}_{\mathit{r}\mathit{e}\mathit{q}}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\quad {\dot{\mathit{r}}}_{\mathit{r}\mathit{e}\mathit{q}}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\quad {\dot{\mathit{w}}}_{\mathit{r}\mathit{e}\mathit{q}}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\mathit{}\quad {\dot{\mathit{u}}}_{\mathit{r}\mathit{e}\mathit{q}}\end{array}\right]}^{T}\\[5pt] {\dot{p}}_{req}=({\phi }_{cmd}-\phi ){K}_{\phi }+({\dot{\phi }}_{cmd}-\dot{\phi }){K}_{\dot{\phi }}\\[5pt] {\dot{q}}_{req}=({\theta }_{cmd}-\theta ){K}_{\theta }+({\dot{\theta }}_{cmd}-\dot{\theta }){K}_{\dot{\theta }}\\[5pt] {\dot{r}}_{req}=({\psi }_{cmd}-\psi ){K}_{\psi }+({\dot{\psi }}_{cmd}-\dot{\psi }){K}_{\dot{\psi }}\\[5pt] {\dot{w}}_{req}=({w}_{cmd}-w){K}_{w}+({\dot{w}}_{cmd}-\dot{w}){K}_{\dot{w}}\\[5pt] {\dot{u}}_{req}=({u}_{cmd}-u){K}_{u}+({\dot{u}}_{cmd}-\dot{u}){K}_{\dot{u}}\end{array}\right\}\end{array}\end{align}

Table 5. Gains of the linear controller

As mentioned before, the velocity in body y-direction ( $v$ ) is not included in the equations since EDFs do not directly generate forces in this direction. To control the aircraft motion in body y-direction, an outer loop controller is used. Outer loop controller generates the required roll angle ${\phi }_{cmd,nav,hover}$ at low speeds to track the body y-velocity command ( ${v}_{cmd}$ ). ${\phi }_{cmd,nav,hover}$ is limited to $\pm $ 30 deg to avoid large bank angles considering the safety and limited thrust authority. Similarly, to control the altitude ( $h$ ) at low speeds, the vertical velocity command ${w}_{cmd,nav,hover}$ is generated by the outer loop altitude controller. To control the flight path angle ( $\gamma $ ) at high speeds, commands of pitch angle ( ${\theta }_{cmd,fpa}$ ) and body velocities, which corresponds to a specific angle-of-attack ( ${w}_{cmd,fpa}$ , ${u}_{cmd,fpa}$ ) and are applied simultaneously. It is also possible to control the angle-of-attack ( $\alpha $ ) by using a similar approach. At high speeds, coordinated turns are performed via applying a yaw angle command ( ${\psi }_{cmd,coord,turn}$ ) based on the coordinated turn equation defined in Ref. (Reference Stevens, Lewis and Johnson46). Details of the command generation logic and the transition speeds are explained in Section 3.2.3.

A manual command (with the subscript $cmd,manual$ ) is also defined for each channel to have a direct pilot command, which is primarily added to test the controller’s performance. To sum up, the attitude and velocity commands sent to the linear controller defined in Equation (18) are expressed as follows. ${\phi }_{cmd},{\theta }_{cmd},{\psi }_{cmd}$ represents the total roll, pitch and yaw angle commands; whereas, ${w}_{cmd}$ and ${u}_{cmd}$ are the total velocity commands in body z and x directions, respectively.

(19) \begin{align} \begin{array}{r}\left.\begin{array}{r}{\phi }_{cmd}={\phi }_{cmd,manual}+{\phi }_{cmd,nav,hover}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad with\quad\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\phi }_{cmd,nav,hover}=({v}_{cmd}-v){K}_{v}+({\dot{v}}_{cmd}-\dot{v}){K}_{\dot{v}}\\[5pt] {\theta }_{cmd}={\theta }_{cmd,manual}+{\theta }_{cmd,fpa}\\[5pt] {\psi }_{cmd}={\psi }_{cmd,manual}+{\psi }_{cmd,coord,turn}\\[5pt] {w}_{cmd}={w}_{cmd,manual}+{w}_{cmd,nav,hover}+{w}_{cmd,fpa}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad with\quad\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{w}_{cmd,nav,hover}=({h}_{cmd}-h){K}_{h}+({\dot{h}}_{cmd}-\dot{h}){K}_{\dot{h}}\\[5pt] {u}_{cmd}={u}_{cmd,manual}+{u}_{cmd,fpa}\end{array}\right\}\end{array}\end{align}

It is noted that the linear controller uses the kinematic relations and independent from the aircraft model. Therefore, tuning can be performed easily, and the gains of the linear controller are given in the following table.

3.3 Control allocation design and integration into the INDI controller

As mentioned in Section 3.1, the system is over-actuated and properly allocating the limited control authority in case of actuator saturation requires specific CA design. Based on the simulation results (Section 4.4), in the presence of severe disturbances (gusts/failures etc.), actuator saturation can cause stability problems if it is not handled properly. The INDI controller generates the incremental control input without considering the actuator limitations in the controller design. This could lead to large errors between the commanded and physically achieved forces and moments in case of severe actuator limitations. The CA becomes active in case of actuator limitations and allocates the limited control authority accordingly to guarantee stable flight.

If the INDI controller generates actuator commands beyond the minimum and maximum thrust and/or tilt angle, then the constrained optimisation problem given as follows is solved iteratively using the active set algorithm defined in Ref. (Reference Harkegard45).

(20) \begin{align} \begin{array}{r}\left.\begin{array}{r@{\quad}r@{\quad}r@{\quad}r} & \mathrm{\Delta }{U}_{CA}=[{\mathrm{\Delta }}_{{T}_{x,fl}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{x,fr}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{x,wl}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{x,wr}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{z,fl}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{z,fr}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{z,wl}}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\mathrm{\Delta }}_{{T}_{z,wr}}{]}^{T}\\[5pt] \underset{\mathrm{\Delta }{U}_{CA}}{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}} & J\left(\mathrm{\Delta }{U}_{CA}\right)=\parallel {W}_{u}\mathrm{}\mathrm{}(\mathrm{\Delta }{U}_{CA}-\mathrm{\Delta }{U}_{des}){\parallel }^{2}+\gamma \parallel {W}_{v}\mathrm{}\mathrm{}({T}_{v,INDI}\mathrm{}\mathrm{}\mathrm{\Delta }{U}_{CA}-\mathrm{\Delta }{v}_{INDI}){\parallel }^{2}\\[9pt] \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\;\mathrm{t}\mathrm{o} & \mathrm{\Delta }{U}_{CA,min}\lt \mathrm{\Delta }{U}_{CA}\lt \mathrm{\Delta }{U}_{CA,max}\\[5pt] \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\;\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\;\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s} & \mathrm{\Delta }{U}_{INDI}\\[5pt] \end{array}\!\!\!\right\}\end{array}\end{align}

The cost function is represented as a quadratic mixed optimisation problem with two parts named control minimisation and error minimisation [Reference Bodson42]. CA focus on error minimisation by prioritising the moment channels over the force channels in case of limited control authority (i.e. actuator saturation). This is achieved by selecting the corresponding weight matrix as ${W}_{v}=diag(\mathrm{1000},\mathrm{}\mathrm{}\mathrm{1000},\mathrm{}\mathrm{}100,\mathrm{}\mathrm{}50,\mathrm{}\mathrm{}50)$ . Note that the weights of ${W}_{v}$ matrix corresponds to $L,M,N,Fz,Fx$ channels, respectively (see ${v}_{INDI}$ in Equation (12)). The weights of each channel are tuned using simulation results, and they can be fine-tuned easily for particular cases since each weight corresponds to one specific channel. Weights of $L$ and $M$ are higher since in case of actuator saturation the limited control authority is allocated to track the roll and pitch moment commands in the first place to guarantee stable flight. If the roll and pitch moment commands are tracked well, the remaining control authority is allocated to first tracking the yaw moment commands and then the force commands. With this prioritisation order, CA makes sure that the aircraft can track the moment commands in extreme flight conditions, such as multiple axis commands with severe disturbances or failures. To illustrate this case, assume that the aircraft is at hover condition. The pilot gives a strong climb command, and at this moment a severe roll moment disturbance (due to the wind/gust or failures) also occurs. However, due to the strong climb command, EDFs are very close to their maximum rpm and there is not enough thrust to handle roll moment disturbance. In other words, the controller’s commands are out of actuator limits, and this leads to deficiency to handle roll disturbance. To resolve this problem, CA gives more priority to roll channel by allocating most of the remaining control authority to overcome the disturbance. This is an example case, but in most of the cases the first priority is always maintaining the stability in roll and pitch channels to avoid vital crashes. For specific cases, weights of ${W}_{v}$ can be adjusted to change the priority order. To conclude, CA has a major role in case of limited control authority. This case is simulated in Section 4.4, and it is observed that CA performance is highly critical and satisfactory. As mentioned, the first part of the cost function is related to the control minimisation. For the sake of simplicity, there is no aim to minimise specific control inputs. Therefore, ${W}_{u}$ matrix is chosen as ${W}_{u}=diag(1,\,1,\,\ldots,\,1)$ , whereas the desired incremental control input is zero, i.e. $\mathrm{\Delta }{U}_{des}=zeros\left(\mathrm{8,1}\right)$ . Finally, the tuning parameter $\gamma $ is taken as 0.0001.

Since the INDI controller generates the incremental control input, the CA needs to be also formulated incrementally. Therefore, evaluation of the minimum and maximum actuator limits (i.e. $\mathrm{\Delta }{U}_{CA,min}$ and $\mathrm{\Delta }{U}_{CA,max}$ in Equation (20)) at each time step is critical. Actuator limits are mapped from the physical control input ( $U$ ) to the INDI control input ( ${U}_{INDI}$ ) based on the thrust vector concept formulated in Equation (12). CA is adapted to the INDI controller formulation and the incremental limits are found as follows.

(21) \begin{align} \begin{array}{l}\left.\begin{array}{l}{\mathrm{\Delta }}_{{T}_{x,min}}={T}_{x,min}-{T}_{x,act}\mathrm{}\mathrm{}\mathrm{}\mathrm{},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad \mathrm{w}ith\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{T}_{x,min}=sin(\pi /2-{\delta }_{max})\sqrt{{T}_{x,act}^{2}+{T}_{z,act}^{2}}\\[5pt] {\mathrm{\Delta }}_{{T}_{z,min}}={T}_{z,min}-{T}_{z,act}\mathrm{}\mathrm{}\mathrm{}\mathrm{},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad \mathrm{w}ith\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{T}_{z,min}=cos(\pi /2-{\delta }_{min})\sqrt{{T}_{x,act}^{2}+{T}_{z,act}^{2}}\\[5pt] {\mathrm{\Delta }}_{{T}_{x,max}}={T}_{x,max}-{T}_{x,act}\mathrm{}\mathrm{}\mathrm{}\mathrm{},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad \mathrm{w}ith\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{T}_{x,max}=\sqrt{{T}_{max}^{2}-{T}_{z,act}^{2}}\\[5pt] {\mathrm{\Delta }}_{{T}_{z,max}}={T}_{z,max}-{T}_{z,act}\mathrm{}\mathrm{}\mathrm{}\mathrm{},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\quad \mathrm{w}ith\quad \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{T}_{z,max}=\sqrt{{T}_{max}^{2}-{T}_{x,act}^{2}}\\[5pt] \end{array}\right\}\end{array}\end{align}

In the above equations, ${\delta }_{min}$ , ${T}_{min}$ and ${\delta }_{max}$ , ${T}_{max}$ are the minimum and maximum actuator limits defined in Table 4, and ${T}_{x,act}$ , ${T}_{z,act}$ represent the x and z components of the thrust vector considering the actuator limits. Note that Equation (21) is given in generic form. It is calculated for each actuator state (i.e. front left, front right, wing left, wing right) specifically since limits are different (see Table 4).

In case of actuator limitations, the constrained optimisation problem is solved online at each time step within a sample time of 0.01 s. Therefore, the maximum number of iterations is set to 50 to limit the computational effort. Simulations are performed on a computer with a 1.8 GHz Intel core CPU. A detailed analysis is not made to check the required computational effort, but it is observed that the algorithm can be solved easily with the current computer technology. Detailed analysis of the active set algorithm’s computational properties compared to the other approaches can be found in Ref. (Reference Harkegard45).

Another critical point is the continuity of the actuator commands when the CA becomes active/inactive. Based on the results given in Section 4, discontinuous actuator commands are not observed. At this point, it is reminded that the initial point of the CA algorithm is the output of the INDI controller (i.e. $\mathrm{\Delta }{U}_{INDI}$ ). Moreover, $\mathrm{\Delta }{U}_{INDI}$ of the overall INDI controller output given in Equation (17) is replaced with the output of CA algorithm $\mathrm{\Delta }{U}_{CA}$ when CA algorithm is active (i.e. in case of actuator saturation), see Fig. 5.

Formulation of the unified INDI controller and integration of the CA method are described in this section. Performance of the overall controller is analysed in the next section via simulation results. A high-level block diagram of the closed loop simulation model is given in Fig. 5.

4.1 Takeoff & transition from hover to cruise flight

As can be seen in Fig. 6, the altitude command is applied at 3 s. Thrust is increased on the front and wing sections while the flap angles are kept at 90° as expected (Fig. 7). A heading rate is applied at 5 s while the aircraft is gaining altitude. Very small changes in the flap angles are observed to track the heading command. This is due to the highly sufficient heading control authority at hover.

Figure 6. Simulation results for takeoff and transition from hover to cruise flight.

Figure 7. Simulation results for takeoff and transition from hover to cruise flight, actuator states.

At 15 s, the aircraft reaches the desired altitude of 40 metres, and the speed command is applied between 15 and 35 s to perform the transition from hover to cruise flight. Aircraft gain speed by tilting the flaps to lower the flap angles and increase the thrust in body x-direction. Around 35 s, the desired cruise speed of 78 m/s is achieved. It is aimed to increase the angle-of-attack in cruise to take advantage of the aerodynamic lift. To achieve 4° cruise angle-of-attack, combined body velocity and pitch angle commands are applied between 35 and 40 s. A slight loss in altitude is observed during transients, but it is recovered at the steady state. The aircraft reaches the cruise trim condition with very small steady state errors.

According to Fig. 7, the CA becomes active frequently when the minimum thrust and flap angle limits are reached on the wing sections at cruise (i.e. after 40 s). The CA prioritises tracking roll and pitch moment commands over tracking yaw moment and force commands based on the weighting matrix ${W}_{v}$ defined in Equation (20). In this way, the limited control authority is adequately allocated to guarantee stable flight in the first place. The same transition test is performed without activating the CA algorithm and purely commanding the INDI controller outputs to the actuators. It is seen that the INDI controller violates the actuator limits significantly when the combined pitch angle and body velocity commands are applied after 40 s. Without proper prioritisation (i.e. without the CA), the INDI controller can not satisfy stable flight in case of actuator saturation. In Section 4.4, a comparison test is performed at hover to show the importance of the CA algorithm in case of limited control authority.

4.2 Climbing/descending and coordinated turn at cruise

As mentioned previously in Section 3.2.3, climbing/descending is achieved by adjusting the flight path angle via combined body velocity and pitch angle commands.

Based on Fig. 8, climb and descent manoeuvers are performed between 45 and 65 s by achieving positive and negative 5° flight path angles, respectively. Coordinated turn manoeuver is also tested at 75 s by applying 30° roll angle commands. Significant altitude loss is observed during the coordinated turns. This is an expected result since some of the total lift is used for turning. An outer loop controller can be designed to minimise the altitude loss similar to the one designed for hover (see ${w}_{cmd,nav,hover}$ in Equation (19)). The outer loop controller should increase total lift to overcome the lift deficiency during turns. In forward flight, this can be done by increasing angle-of-attack and/or airspeed. Some level of altitude loss might be inevitable due to the lack of control authority if bank angle is too high.

Figure 8. Simulation results for climbing/descending and coordinated turn at cruise.

According to Fig. 9, actuator saturations are observed several times during these manoeuvers. The CA works as expected to guarantee stable flight with the prioritisation of moment channels. It is noted that the CA algorithm can find a solution with a maximum number of iteration less than 10.

Figure 9. Simulation results for climbing/descending and coordinated turn at cruise, actuator states.

4.3 Transition from cruise to hover flight and landing

Results of transition from cruise to hover flight and vertical landing are shown in Figs. 10 and 11. At 45 s, combined body velocity and pitch angle commands are applied to achieve zero angle-of-attack. The body x-velocity command is then applied between 50 and 70 s to achieve hover flight while keeping the altitude. Vertical landing is started at 70 s, and also aircraft heading is adjusted during the landing.

Figure 10. Simulation results for transition from cruise to hover flight and landing.

Figure 11. Simulation results for transition from cruise to hover flight and landing, actuator states.

Based on Fig. 11, flap angles of both front and wing sections reach maximum limits during the transition to decelerate the aircraft. Therefore, the CA becomes active many times to prioritise moment commands over force commands. This is an expected result considering the significant deceleration command during the cruise to hover transition.

4.4 Importance of the CA: severe roll moment disturbance at hover

As mentioned previously, the CA has a very crucial role in case of actuator saturation. To show the CA’s effectiveness, a test case is simulated when the CA is on (i.e. INDI+CA) and off (i.e. only the INDI controller). As shown in Fig. 12, a strong roll moment disturbance is applied at hover (at 3 s) for both conditions (i.e. CA ON and CA OFF). Simultaneously, the pilot wants to gain altitude and apply body z-velocity command ( ${w}_{cmd}$ ).

Figure 12. Comparative simulation results that shows the importance of the CA.

According to Fig. 13, thrust authority is not enough to track the pilot command and reject the roll disturbance at the same time. The INDI controller generates thrust commands above the maximum physical limits (see dashed and red lines in Fig. 13). Therefore, disturbance rejection properties are not as desired due to the saturation related discrepancy between the commanded and physically achieved actuator states. When CA is OFF, approximately 25° roll angle is observed due to the disturbance (Fig. 13).

However, when the CA is active, disturbance rejection properties are much better since the CA prioritises the moment commands over the force commands considering the actuator saturation. There is no discrepancy between the physically achieved and commanded thrust when the CA is on (see dashed and blue lines in Fig. 13). The trade-off is the worse tracking of the body z-velocity command, as seen in the $w$ plot of Fig. 12. The deteriorated tracking performance of the body z-velocity command can be accepted since allocating the limited control authority to the moment channels is more important regarding the stability. When CA is ON, roll moment disturbance only cause 7–8 deg roll angle, which is much less than the case when CA is OFF (Fig. 12).

Figure 13. Comparative simulation results that shows the importance of the CA, actuator states.

To conclude, CA performs very well and allocates the limited control authority properly to guarantee stable flight in case of severe disturbance. More detailed analysis of the CA performance are given in Ref. (Reference Suicmez27).