Nomenclature
-
$L,M,N$
-
aerodynamic moments along the body axis,
$ft.lbs$
-
${{\Delta }}L,{{\Delta }}M,{{\Delta }}N\;$
-
incremental aerodynamic moments along the body axis,
$ft.lbs$
-
$\alpha, \beta $
-
angle of attack, side slip angle,
$rad$
-
$\dot \alpha, \dot \beta $
-
the derivative of angle of attack, side slip angle,
$rad/s$
-
$\phi, \theta, \psi $
-
roll, pitch and yaw angles along the body axis with respect to the Earth frame,
$rad$
-
${\phi _{ref}},{\theta _{ref}},{\psi _{ref}}$
-
reference roll, pitch and yaw angles along the body axis with respect to the Earth frame,
$rad$
-
$p,q,r$
-
pitch and yaw rate along the body axis,
$rad/s$
-
$\dot p,\dot q,\dot r$
-
the derivative of pitch and yaw rate along the body axis,
$rad/{s^2}$
-
${A_x},{A_y},{A_z}$
-
accelerations along the body axis,
$rad/{s^2}$
-
$\lambda $
-
the angle between path tangential angle and the aircraft heading, rad
-
${p_{ref}},{q_{ref}},{r_{ref}}$
-
reference pitch and yaw rate along the body axis,
$rad/s$
-
${I_{xx}},{I_{yy}},{I_{zz}},{I_{xz}},{I_{yz}},{I_{xy}}$
-
inertia parameters,
$slg - f{t^2}$
-
${v_{des}},u$
-
control variables, control effectors for control allocation
-
${B_c},S$
-
control effective matrix, steady-state distribution matrix
-
${W_v},{W_1},{W_2}$
-
weight matrices for control allocation
-
${\gamma _w}$
-
weighting factor for control allocation
-
${\gamma _{ref}}$
-
reference flight path angle,
$rad$
-
${\alpha _k}$
-
path tangential angle,
$rad$
-
$e\left( t \right)$
-
cross-track error,
$meter$
-
$\Delta $
-
look-ahead distance,
$meter$
-
${\gamma _{sensed}}$
-
measured flight path angle,
$rad$
-
$\dot \gamma $
-
the derivative of measured flight path angle,
$rad/s$
-
${\rm{\Gamma }}$
-
glide slope deviation angle,
$rad$
-
${\mu _{ref}}$
-
reference bank angle,
$rad$
-
${\dot \mu _{ref}}$
-
the derivative of the reference bank angle,
$rad/s$
-
${\mu _{sensed}}$
-
measured bank angle,
$rad$
-
${\gamma _c}$
-
convergence value for LOS
-
${\chi _{sensed}}$
-
measured azimuth angle,
$rad$
-
${\chi _{ref}}$
-
desired course angle,
$rad$
-
$\dot \chi $
-
the derivative of the course angle,
$rad/s$
-
${\chi _r}$
-
velocity path relative angle,
$rad$
-
${V_{ref}}$
-
reference airspeed,
$ft/s$
-
${V_{rs}}$
-
projected ground speed of the actual aircraft,
$ft/s$
-
${\delta _{{e_{lo}}}},{\delta _{{e_{ro}}}},{\delta _{{e_{li}}}},{\delta _{{e_{ri}}}}$
-
outer left, outer right, inner left, inner right elevators,
$rad$
-
${\delta _{{a_r}}},{\delta _{{a_l}}},{\delta _{{r_{up}}}},{\delta _{{r_{low}}}}$
-
right, left ailerons and upper, lower rudders,
$rad$
1.0 Introduction
During standard terminal arrival and approach in the terminal control area, path following or tracking is critical because these trajectories are constructed to avoid man-made and geographical obstacles. If the aircraft deviates too much from the stated path, the consequences can be fatal. For this reason, there are different types of algorithms for guidance and navigation. After, reference commands are transmitted to the flight control system. Also, for flight control systems or computers, linear and nonlinear algorithms can be used for highly nonlinear aircraft.
Nonlinear Dynamic Inversion (NDI) has taken a lot of attention from the aerospace community. In the literature for other fields, it is known as feedback linearisation (FBL) or input-output linearisation. However, for NDI, one differentiation of each control variable is required for inversion [Reference Acquatella, Falkena, van Kampen and Chu1]. Namely, a system is linearised inherently via an NDI closed loop. Nonlinearities in the system can be reduced by state feedback. As a disadvantage of NDI, it is an explicit control method. The controller still needs an explicit and accurate model of the nonlinear aircraft [Reference Acquatella, Falkena, van Kampen and Chu1]. If there are model mismatches and measurement errors, the controller performance is degraded and the consequence can be unstability. In other words, aerodynamic uncertainties have to be taken into account during the controller design due to the degradation of that controller [Reference Singh, Steinberg and Page2]. Also, even if model uncertainties are reduced with the help of online parameter identification, online self-excitation is required [Reference Lombaerts, Chu, Mulder and Joosten3]. The best way to solve model mismatch problems is to design a controller that is less dependent on the model [Reference Sieberling, Chu and Mulder4]. For this reason, Incremental Nonlinear Dynamic Inversion (INDI) is used.
The INDI method was initially introduced as modified NDI or simplified NDI [Reference Smith5]. For INDI, angular accelerations are utilised in place of model aerodynamic derivatives. INDI is also known as the implicit control method. This sensor-based approach compensates for wind gust disturbances and unmodelled dynamics inherently. As a disadvantage of INDI, there can be sensor errors, delays and biases during measurement or prediction of the angular acceleration.
There are different industry applications for NDI and INDI. The simplicity and robust design of this method have taken more attention year by year [Reference Balas6]. In the UK-VAAC (vectored thrust aircraft advanced control) program, Proportional–Integral (PI) is used with NDI for fly-by-wire (FBW) VAAC Harrier [Reference Smith and Berry7]. Boeing X-36 has a reconfigurable flight control that consists of dynamic inversion, and it is used with an online neural network to cancel inversion errors. Honeywell Research Center’s approach to NDI is MACH (multi-application control). It is modular and blends classical controller design with inversion. The outer loop is PI, and the inner is dynamic inversion [Reference Balas6]. NDI is used for Lockheed Martin JSF F-35. The F-35 is the first fighter using dynamic inversion from zero airspeed to supersonic speed, high angle-of-attack, during short takeoff and vertical takeoff and landing [Reference Harris and Stanford8]. For the full-scale advanced systems testbed (FAST) F/A-18 airplane, NDI has been developed as a baseline controller. Besides, the NDI controller is used in the hardware-in-the-loop simulation and is tested during an actual flight [Reference Miller9]. In [Reference Enns, Bugajskı, Hendrick and Stein10], NDI is used for F-18 high angle-of-attack research vehicle.
In [Reference Lu, van Kampen, de Visser and Chu11], INDI and NDI are used as fault-tolerant trajectory control during model uncertainties and actuator faults. In [Reference Lombaerts, Kaneshige, Schuet, Aponso, Shish and Hardy12], INDI and Proportional–Integral–Derivative (PID) controllers are used for the eVTOL vehicle, and different modes can be selected for hover, transition, and forward flight. In [Reference Wang, van Kampen, Chu and Lu13], the attitude controller is designed based on NDI, INDI and pseudo control hedging (PCH) in addition to the first-order linear model. In [Reference Acquatella, Falkena, van Kampen and Chu1], INDI is designed for the spacecraft attitude controller in the event of disturbance and uncertainty. In [Reference Sieberling, Chu and Mulder4], INDI is compared with NDI for aerodynamic uncertainties and inertia, centre of gravity (CG) mismatches. In [Reference Grondman, Looye, Kuchar, Chu and Kampen14], INDI is developed for the Cessna Citation II experimental aircraft. In addition, after conducting a test flight, INDI is compared with NDI. As a result, INDI outperforms NDI by far. In this work, the INDI is chosen as the rate controller.
In this paper, NASA GTM (generic transport model) (T-2) is used as a testbed [15]. It is a
$5.5{\rm{\% }}$
geometrically and dynamically scaled model of one of the transport aircraft [Reference Jordan, Langford and Hill16]. It is jet-powered and remotely controlled [Reference Cunningham, Cox, Murri and Riddick17]. The NASA team carried out comprehensive wind tunnel tests in the NASA Langley Tunnel to obtain model parameters (Fig. 1).
NASA Generic Transport Model (GTM) T-2 design. Reprinted by permission of the American Institute of Aeronautics and Astronautics, Inc., from [Reference Cunningham, Cox, Murri and Riddick17]. AIAA.

1.1 Main contribution
The cross-track error is not the sole error for establishing the predefined path. The error between the aircraft heading and the path course is required to describe how fast heading change occurs to re-establish the path again. In other words, after deviation, aggressive or gentle heading change to the predefined path has to be computed. We believe that this error has not been addressed in previous publications such as [Reference Fossen18–Reference Fossen20]. The relevant numerical simulations are presented to demonstrate its importance. Second, the mode control panel and realistic autopilot architecture have not been developed together in any single paper in the literature. It is advantageous since the related controllers can be evaluated efficiently when utilised collectively. We can provide a detailed explanation of the contribution:
-
1) This paper is an extension of [Reference Lu, van Kampen, de Visser and Chu11]. More precisely, the contribution of our work includes detailed lateral navigation in the terminal phase in addition to the approach and flare phases. Compared to the architecture in [Reference Lu, van Kampen, de Visser and Chu11], which has four loops, this work comprises two more loops. Furthermore, dynamic control allocation is implemented for the overactuated aircraft-NASA GTM design. The LOS algorithm is used for lateral navigation and complete instructions are provided to integrate it into the NASA GTM design. Even if LOS is used in [Reference Breivik, Hovstein and Fossen19], various design parameters are adjusted for NASA GTM design in this study, and the LOS algorithm has an adaptive feature. In addition, the virtual aircraft is used to compute cross-track error during lateral navigation. Virtual aircraft’s airspeed is computed with respect to actual aircraft. There is a sample case of lateral navigation. The aircraft makes an approach to the runway after leaving the terminal phase. During this process, the glide path is captured automatically. The approach and the flare phases are simulated. Detailed numerical values, block diagrams and results are provided.
-
2) A complete navigation and guidance flight is developed by combining the INDI, NDI, PID and LOS algorithms. Our paper contains a background of the general underlying control methodology, a more detailed description of the applied control architecture, and simulation results.
-
3) For distributing desired moments to control surface actuators, dynamic control allocation (CA) is used. INDI generates the desired control variables (desired moments), and these values are fed to CA. The weighted least square (WLS) equation is developed and solved with the active set method for CA [Reference Härkegård21]. In [Reference Lu, van Kampen, de Visser and Chu11], control effective matrix is 3 by 3, while in this work it is 3 by 8 because there are additional control surfaces to coordinate. Control issues become more challenging as there are more control surfaces.
-
4) In this work, a comprehensive mode control panel is developed for the NASA GTM model, and approximately all autopilot modes are proposed from a pilot’s point of view. A mode control panel is an interface between the pilot/operator and autopilot. This is not part of the flight control or navigation/guidance computers, but it directly affects aircraft control. In addition, the absence of a mode control panel makes the autopilot architecture unrealistic. It is not included in many papers in the literature. Therefore, constructing a mode control panel for simulations is essential for realistic analysis. With the mode control panel, the results of the given commands can be presented more realistically. In the future, it will be simple to achieve fully autonomous flight with the aid of the mode control panel modes.
-
5) In the LOS algorithm, virtual aircraft is used. The actual aircraft tracks the path of the virtual aircraft. At that time, virtual aircraft flies between waypoints. The equations are taken from [Reference Lekkas and Fossen22]. However, in this work, the virtual aircraft flies with the aircraft’s projection airspeed. It is a different approach. Equations to calculate this projected speed are provided. The purpose of this calculation is to ensure that, with the help of the virtual aircraft, the cross-track error line becomes perpendicular. If the virtual aircraft does not fly at the projected airspeed, the cross-track error line does not become perpendicular.
-
6) There is a lookahead distance (
$\Delta $
) value for the LOS algorithm. That distance varies between
${\Delta _{min}}$
and
${\Delta _{max}}$
. These minimum and maximum parameters are used to either drive the aircraft gently if it is too close to the path or aggressively if it is too far from the path [Reference Lekkas and Fossen23]. However, in addition to distance, the angle between the aircraft and the path is also important. If this angle is too much, but the distance of the aircraft is too close to the path, logic has to be developed to improve lateral navigation. As a result, LOS
${\Delta _{max}}$
value changes in an adaptive manner regarding the angle between the aircraft heading and the path tangential angle. This adjustment is also beneficial for path following during control surface faults. We believe that this method has never been applied before. -
7) Also, single-engine operations, control surface lock-in-place and turbulence/wind conditions are considered in the simulation. Additionally, a variety of faults and failures are examined to demonstrate the robustness of the autopilot architecture.
The paper is organised as follows: In Section 2.0, background information is provided on the controllers, which are NDI, INDI and Dynamic CA, respectively. Next, Section 3.0 presents the simulation architecture overview. It contains the mode control panel and different types of controllers. In Section 4.0, numerical simulations are provided to illustrate the results. Finally, conclusions are presented.
2.0 Background for controllers
2.1 Nonlinear dynamic inversion (NDI)
In this section, a summary of the formulation for NDI is given. Consider the nonlinear system as:
In this equation,
$f\left( x \right)$
represents nonlinear dynamics of the aircraft,
$B\left( x \right)$
is the control matrix,
$u$
is the control effectors or commands. Finding a direct relationship between the input and the intended output and inverting it is the idea behind NDI. For control purposes, with the help of the NDI given in [Reference Lombaerts, Kaneshige, Schuet, Aponso, Shish and Hardy12],
$u$
becomes:
$v$
is a virtual input vector. It can be desired accelerations. In Equation (2),
$f\left( x \right)$
consists of aerodynamic derivatives of the aircraft. In other words, aircraft model is required. The inverse of the
$B\left( x \right)$
matrix must be calculated. If it is not a square matrix, the inversion of this matrix can be computed by control allocation (CA) methods. NDI method depends on the aerodynamic model due to
$f\left( x \right)$
.
2.2 Incremental nonlinear dynamic inversion (INDI)
The nonlinear system dynamics is given in Equation (1). Taylor series expansion of the system can be taken at the last sensor sampling moment (
$t - {{\Delta }}t$
) [Reference Wang, van Kampen, Chu and Lu13]:
\begin{align}\dot x &\approx {\dot x_0} + {\left. {\frac{{\partial \left[ {f\left( x \right) + B\left( x \right)u} \right]}}{{\partial x}}} \right|_{x = {x_0},u = {u_0}}}\left( {x - {x_0}} \right)\nonumber\\ &\quad+ {\left. {\frac{{\partial \left[ {f\left( x \right) + B\left( x \right)u} \right]}}{{\partial u}}} \right|_{x = {x_0},u = {u_0}}}\left( {u - {u_0}} \right) + {\mathcal{O}}\left( {{{\Delta }}{x^2}} \right)\end{align}
\begin{align}\dot x \approx {\dot x_0} + {\left. {\frac{{\partial \left[ {f\left( x \right) + B\left( x \right)u} \right]}}{{\partial x}}} \right|_{x = {x_0},u = {u_0}}}{{\Delta }}x + B\left( {{x_0}} \right){{\Delta }}u\\[8pt]\nonumber\end{align}
where
${{\Delta }}x = \left( {x - {x_0}} \right)$
and
${{\Delta }}u = \left( {u - {u_0}} \right)$
are the state and control increments in one sampling time step (
${{\Delta }}t$
). In Equation (4),
${\dot x_0}$
can be computed by direct measurement or state estimation. Also,
${u_0}$
can be computed directly. The variation
$\left( {x - {x_0}} \right)$
can be neglected in principle. Furthermore Equation (4) can be modified as:
After mathematical manipulation of Equation (5) and
${{\Delta }}u$
is represented as
${{\Delta }}{u_{ref}}$
, the control law becomes:
INDI is the incremental form of the NDI. In Equation (6), the derivative of the state vector at the last time step (
${\dot x_0}$
) is used instead of
$f\left( x \right)$
[Reference Lu, van Kampen, de Visser and Chu11]. For example, angular accelerations (
${\dot p_0}$
,
${\dot q_0}$
,
${\dot r_0}$
) can be used.
The INDI is shown in Fig. 2. Also, we can write
${v_{des}} = v - {\dot x_0}$
. The increment input is
${{\Delta }}{u_{ref}}$
and
${u_0}$
is the input from the previous time step. The total input has to be calculated:
Incremental nonlinear dynamic inversion.

As stated before, INDI replaces the nonlinear dynamics
$f\left( x \right)$
with a simplified estimation, typically based on partial derivatives or measurements, to achieve control without full knowledge of the dynamics. Mostly, angular accelerations are utilised for INDI, and it is an implicit control method. It compensates for wind gust disturbances and unmodelled dynamics inherently. Measured angular accelerations (
${\dot x_0}$
) are fed to the INDI controller, so aerodynamic derivatives are not required to obtain these incremental changes. Different approaches can be used for obtaining angular accelerations (
${\dot x_0}$
). For example, in [Reference Sieberling, Chu and Mulder4], a filter with a simple prediction polynomial is used to obtain the angular accelerations. In [Reference Bacon, Ostroff and Joshi24], second-order low pass filters are used for the same purpose. In our work, we use the same filter given in Equation (8) for obtaining angular acceleration [Reference Bacon, Ostroff and Joshi24].
with
${\omega _n} = 25$
rad/s and
${\zeta _n} = 0.8$
[Reference Lu, van Kampen, de Visser and Chu11].
As a disadvantage of INDI, there can be sensor errors, delays and biases during measurement or prediction of the angular acceleration. Prediction filters can compensate for measurement delays, but the model is required for these filters. It is assumed that there are no sensor errors in our work.
In [Reference Wang, van Kampen, Chu and Lu13], stability and robustness analysis of INDI is done by Lyapunov methods and nonlinear system perturbation theory.
2.3 Dynamic control allocation
The control allocation problem can be given as
${v_{des}} = {B_c}{{\Delta }}{u_{ref}}$
. Instead of
$B\left( x \right)$
matrix given in Equation (1),
${B_c}$
matrix is used. The control effective matrix
${B_c}$
used for the relation between the moments and the control surfaces, is given in Equation (13).
${v_{des}}$
is modified as
${v_{des}} = I\left[ {v - {{\dot x}_0}} \right]$
and
$I$
is the inertia matrix. In this work, the
${B_c}$
matrix is not square (3x8 matrix). The formulation can be rewritten to exclude the inertia matrix (
$I$
) from the
$B$
matrix because we choose
${v_{des}}$
as a desired moment.
Since the
${B_c}$
matrix cannot be inverted directly, the dynamic control allocation is used. The controller scheme shown in Fig. 2 can be modified using Fig. 3 due to
${B_c}$
matrix. The incremental input, given in Equation (9) is used in Fig. 3.
Modified incremental nonlinear dynamic inversion.

${B_c}$
matrix gives the relation between desired moments (
${{\Delta }}L$
,
${{\Delta }}M$
,
${{\Delta }}N$
) and control effector (
${\delta _{{e_{lo}}}},{\delta _{{e_{li}}}}$
etc.).
${v_{des}}$
is the desired control variables (desired moments).
${{\Delta }}{u_{ref}}$
is the control effector vector.
In [Reference Berge and Fossen25], the term ‘dynamic’ is introduced. In [Reference Härkegård26], a novel systematic approach is proposed to dynamic control allocation by including an additional term for the actuator rates in the optimisation criterion. In this paper, the same approach is employed for the control allocation problem. Namely, actuator rates are penalised with rate limiting during the real-time calculation. For additional readings on dynamic control allocation, see [Reference Härkegård26] for the equations, and [Reference Ergöçmen27] for the fault tolerant dynamic control allocation.
Weighted least squares (WLS) is chosen for the dynamic control allocation due to fast convergence for calculating
${{\Delta }}{u_{ref}}$
given in Equation (10). The parameters given in Equations (11), (12) are used in Equation (10).
\begin{align}\Delta u_{ref}&= arg{ \min_ {\underline{u }\leq\Delta u \leq\bar{u}}}\|{{W_u}\left( {{{\Delta }}u\left( t \right) - {u_p}\left( t \right)} \right)} \|^2 \nonumber\\ &\quad + {\gamma _w}{\kern.2pt}\| {{W_v}\left( {{B_c}{{\Delta }}u\left( t \right) - {v_{des}}\left( t \right)} \right)} \|{^2} \end{align}
In these equations,
${{\Delta }}u\left( t \right)$
is the control effector,
${u_p}\left( t \right)$
desired steady-state control input and
${{\Delta }}{u_{prev}}$
is the control effector from the previous sampling time.
${W_1}$
,
${W_2}$
and
${W_v}$
are weight matrices. While greater
${W_2}$
values slow down (penalise) actuator movement, bigger
${W_1}$
values are employed to quickly converge to the desired value.
${W_v}$
matrix is a weight matrix for control variables. Greater values mean significance and enable rapid convergence.
$\left( {{{\Delta }}u\left( t \right) - {\rm{\;\;}}{u_p}\left( t \right)} \right)$
is the control error and
$\left( {{B_c}{{\Delta }}u\left( t \right) - {v_{des}}\left( t \right)} \right)$
moment error.
$S$
matrix is a steady-state distribution matrix. It is for desired control purposes, for example, minimum deflection. For instance, the corresponding row of the S matrix is set to
$0$
if some of the control surfaces are required to move less during normal flight.
${u_p}\left( t \right)$
is the desired steady-state control input.
${\gamma _w}$
is weighting factor, which is used to give importance between control and moment error.
${B_c}$
is the control effective matrix. By
${B_c}$
matrix, the efforts are transformed from desired control variables (
${v_{des}}$
) (moments increments:
${{\Delta }}{L_{control}}$
,
${{\Delta }}{M_{control}}$
,
${{\Delta }}{N_{control}}$
) to control effectors (
${{\Delta }}{u_{ref}}$
) shown in Fig. 3. Control effectors are respectively given as: left outer, left inner, right outer, right inner elevators (
${\delta _{{e_{lo}}}}$
,
${\delta _{{e_{li}}}}$
,
${\delta _{{e_{ro}}}}$
,
${\delta _{{e_{ri}}}}$
); right and left ailerons (
${\delta _{{a_r}}}$
,
${\delta _{{a_l}}}$
); upper and lower rudders (
${\delta _{{r_u}}}$
,
${\delta _{{r_l}}}$
).
${B_c}$
matrix is given in Equation (13).
\begin{align}{B_c} = \left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}}0 {}&0 {}&0 {}&0 {}&{\frac{{\partial L}}{{\partial {\delta _{{a_r}}}}}} {}&{\frac{{\partial L}}{{\partial {\delta _{{a_l}}}}}} {}&{\frac{{\partial L}}{{\partial {\delta _{{r_{up}}}}}}} {}&{\frac{{\partial L}}{{\partial {\delta _{{r_{low}}}}}}}\\{\frac{{\partial M}}{{\partial {\delta _{{e_{lo}}}}}}} {}&{\frac{{\partial M}}{{\partial {\delta _{{e_{li}}}}}}} {}&{\frac{{\partial M}}{{\partial {\delta _{{e_{ro}}}}}}} {}&{\frac{{\partial M}}{{\partial {\delta _{{e_{ri}}}}}}} {}&0 {}&0 {}&0 {}&0\\0 {}&0 {}&0 {}&0 {}&{\frac{{\partial N}}{{\partial {\delta _{{a_r}}}}}} {}&{\frac{{\partial N}}{{\partial {\delta _{{a_l}}}}}} {}&{\frac{{\partial N}}{{\partial {\delta _{{r_{up}}}}}}} {}&{\frac{{\partial N}}{{\partial {\delta _{{r_{low}}}}}}}\end{array}} \right]\end{align}
The related derivatives in the
${B_c}$
matrix are given as:
The weight is
$49.6$
$lb$
. Significant lengths are
$\bar c = 0.9153ft$
is the the mean aerodynamic chord,
${S_a} = 5.9018f{t^2}$
is the wing area.
$b = 6.8488ft$
is the wing span.
$q$
is the dynamic pressure. The inertia parameters are given as (
$slug - f{t^2}$
):
The dimensionless aerodynamic control derivatives are given as:
\begin{align*}{C_{{l_{{\delta _{{a_r}}}}}}} &= - 0.0760{\rm{\;\;\;\;\;\;\;\;}}{C_{{n_{{\delta _{{a_r}}}}}}} = - 0.0025\\ {C_{{l_{{\delta _{{a_l}}}}}}} &= 0.0760{\rm{\;\;\;\;\;\;\;\;}}{C_{{n_{{\delta _{{a_l}}}}}}} = 0.0025\\ {C_{{l_{{\delta _{{r_{low}}}}}}}} &= 0.0145{\rm{\;\;\;\;\;\;\;\;}}{C_{{n_{{\delta _{{r_{low}}}}}}}} = - 0.0846\\ {C_{{l_{{\delta _{{r_{up}}}}}}}} &= 0.0145{\rm{\;\;\;\;\;\;\;\;}}{C_{{n_{{\delta _{{r_{up}}}}}}}} = - 0.0846\\ {C_{{m_{{\delta _{{e_{lo}}}}}}}} &= - 0.4436{\rm{\;\;\;\;\;\;\;\;}}{C_{{m_{{\delta _{{e_{ro}}}}}}}} = - 0.4436\\ {C_{{m_{{\delta _{{e_{li}}}}}}}} &= - 0.4436{\rm{\;\;\;\;\;\;\;\;}}{C_{{m_{{\delta _{{e_{ri}}}}}}}} = - 0.4436 \end{align*}
All parameters are given to compute
${B_c}$
, the control effective matrix. Other numerical values for
$S$
matrix (steady-state distribution matrix),
${W_1}$
,
${W_2}$
and
${W_v}$
weight matrices, which are heuristically tuned in the control allocation are given as:
\begin{align}S = \left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}0 {}&1 {}&0\\0 {}&1 {}&0\\0 {}&1 {}&0\\0 {}&1 {}&0\\{3.8} {}&0 {}&{ - 1.7}\\{ - 3.8} {}&0 {}&{1.7}\\{ - 0.1} {}&0 {}&{12}\\{ - 0.1} {}&0 {}&{12}\end{array}} \right]\begin{array}{*{20}{l}}{{\rm{\;\;\;\;\;\;\;\;}}} {}{{W_1} = diag\left( {3,3,3,3,1,1,1,1} \right)}\\{{\rm{\;\;}}} {}{}\\{{\rm{\;\;\;\;\;\;\;\;}}} {}{{W_2} = diag\left( {3,3,3,3,1,1,1,1} \right)}\\{{\rm{\;\;}}} {}{}\\{{\rm{\;\;\;\;\;\;\;\;}}} {}{{W_v} = diag\left( {2,2,2} \right)}\\{} {}{}\end{array}\end{align}
The actuator limits are given in Table 1.
Control effectors position and rate limits

3.0 Simulation architecture overview
The overview of the simulation is shown in Fig. 4. The NASA GTM T2, the nonlinear aircraft model, is used as an aircraft model.
Components of the simulation.

The Fault Injection block is used to initiate an emergency. These situations include control surface lock-in-place and engine failure or fault. With the help of this block, the controller’s robustness can be evaluated. The Atmospheric Environment block consists of atmospheric values, and is used for adjusting wind and turbulence. The Mode Control block is the interface between the pilot/operator and the autopilot. There are different modes with respect to lateral (roll channel) and vertical plane (pitch channel). In addition, there are approach and airspeed hold modes.
Navigation and Guidance block consists of path following and approach controller (PIDs and line-of-sight algorithm), flight path controller (PIDs), and attitude controller (NDI). The path following and approach controller includes both pitch and roll channels. In roll channel-lateral guidance, the LOS algorithm is utilised. In the pitch channel, PIDs are utilised. Flight path controller consists of PIDs. NDI is employed in attitude controller. An aircraft model is not necessary for this NDI as explained in the attitude controller section.
The flight control computer consists of a rate controller (INDI) and control allocation (dynamic CA). Dynamic CA takes into account the maximum and minimum rate and angle limits of the aircraft. In this part, aerodynamic control derivatives are required for allocating desired moments to control surfaces. Each of these blocks and their algorithms is described individually.
3.1 Mode control block
In this work, almost all modes are simulated for the autopilot. Pitch channel, roll channel, approach (APP), and airspeed hold are the major channels for mode control block, as shown in Fig. 5.
Mode control block.

In roll channel, (1) heading hold (HDG HLD) mode and (2) NAV/LNAV mode can be activated. After pushing the HDG HLD button, the reference value can be adjusted and labelled as
${\psi _{ref}}$
. If the NAV/LNAV mode is activated, predetermined waypoints are retrieved from the database. A path is constructed with the help of these waypoints. It is assumed that this path was uploaded before, and the pilot selects this arrival and approach path from the flight management computer (FMC).
In pitch channel, (1) altitude hold (ALT HOLD) mode, (2) vertical navigation (VNAV) mode, and (3) flight path angle (FPA) mode can be activated. After pushing the ALT HOLD button, the value can be selected. The value is labelled as
${h_{ref}}$
. VNAV mode can be activated by pushing this button, but when the aircraft reaches a predetermined waypoint, it begins descent. This path information is in the FMC and can be retrieved by the pilot. After pushing the FPA (another adjustable mode) button, the value becomes adjustable (labelled as
${\gamma _{ref}}$
).
APP button is adjusted and the criteria are fulfilled, NAV/ LNAV and VNAV are activated automatically. For example, if the aircraft is in ALT HOLD mode, but the APP button is pushed, then the mode will be triggered. Once the aircraft reaches the specific waypoint for descent activation, the system will automatically activate VNAV mode. As a result, ALT HOLD is not functional, and the aircraft starts to approach a runway in this example. The same process is valid for HDG HOLD. After a specific waypoint for activation, it changes to NAV/LNAV mode. In summary, when the APP is pushed, it is in an armed position. However, if the criteria are met, it becomes the activated one.
3.2 Path following and approach controller
After the mode control block, these reference signals are transmitted to the roll and pitch channel blocks, which are shown in Fig. 6.
Path following and approach controller.

In pitch channel, there is a vertical activation logic command or approach activation signal route. With the help of this signal, different modes can be activated and these modes can be selected by the mode control block. Finally, a reference flight path angle (
${\gamma _{ref}}$
) activation logic command or approach activation signal route. In roll channel, there is a lateral activation logic command or approach activation signal route. After that, depending on the mode which is selected from the mode control block,
${\psi _{ref}}$
is generated. In the next sections, each of these channels is described individually.
3.2.1 Roll channel-lateral navigation
LNAV/NAV: This part is mainly about guidance and LOS. The guidance system generates a command for the flight control computer (FCC).
In this work, the main objective of the guidance system is path following, which only has spatial constraints but no temporal constraints. However, for path tracking, both constraints are employed. Different techniques for creating a reference trajectory or path have been described in the literature. For example, some of them are developing a reference model, optimal trajectory generation, generating a time-varying reference trajectory, etc. given in [Reference Fossen18]. For path following, the main objective is to track a vehicle that moves along a predefined path [Reference Breivik, Hovstein and Fossen19]. In this work, virtual aircraft is used to generate a path. Namely, the aircraft tries to follow the virtual aircraft, which virtually flies along the predefined path. There are different types of guidance techniques for path following, such as carrot following algorithm, non-linear guidance law (NLGL), vector field, optimal-based, learning-based, geometric, etc. [Reference Rubí, Pérez and Morcego28, Reference Sujit, Saripalli and Sousa29].
Geometric guidance techniques are (1) pure pursuit (PP), (2) LOS, and (3) constant bearing (CB). In this work, LOS is used as a guidance technique for minimising cross-track error problem [Reference Papoulias30]. Besides, two guidance principles can be used for the LOS technique: (1) enclosure-based steering, (2) lookahead-based steering. Due to its easy implementation and several advantages, lookahead-based steering is used in this work. There are different modifications and applications for LOS [Reference Curry, Lizarraga, Mairs and Elkaim31].
The measured azimuth angle (
${\chi _{sensed}}$
) of the aircraft is computed by
${\dot x_{sensed}}$
and
${\dot y_{sensed}}$
:
Defining a straight-line path with two waypoints is given as
$p_k^n = \left( {{x_0},{y_0}} \right)$
and
$p_{k + 1}^n = \left( {{x_1},{y_1}} \right)$
which are shown in Fig. 7. Also, with respect to the path-fixed reference frame,
${x_n}$
is rotated by an angle whose name is path tangential angle (
${\alpha _k}$
) and corresponding equation is given:
Line of sight (LOS) navigation.

There are errors if the aircraft position deviates from the predefined path. These errors (cross-track error and along-track distance, respectively) are
$\varepsilon \left( t \right) = {[e\left( t \right){\rm{\;\;\;\;}}s\left( t \right)]^T}$
and corresponding equations are given:
From Equations (23), (24) (vir. means virtual):
For path following, minimising cross-track error is enough. Path tracking, on the other hand, is used when the problem is both minimising cross-track error
$e\left( t \right)$
and along-track distance
$s\left( t \right)$
. In this work, only the path following is preferred which means cross-track error is used. The along track error does not increase due to the same ground speed between virtual and the real aircraft. Namely, along track error is always zero.
Finally, the primary goal of the guidance is
${\rm{li}}{{\rm{m}}_{t \to \infty }}e\left( t \right) = 0$
.
There are different solutions for calculating cross-track error
$e\left( t \right)$
. For example, one of them is calculated with the geometric formulation. Instead, in this paper, virtual aircraft is used for the calculation. Virtual aircraft position is
${P_{virtual}}$
shown in Fig. 7, and virtual aircraft kinematics for path generation is given:
The virtual aircraft flies with the NASA GTM aircraft’s projected ground speed (
${V_{rs}}$
). In Equations (27), (28),
${V_{rs}}$
is the projected ground speed of the actual aircraft along the path. For example, in Fig. 7 the aircraft does not fly through the path heading, but the virtual aircraft continues to fly through the path with a path heading. As a result, if there is a difference between the aircraft airspeed and the virtual aircraft airspeed throughout the path, the cross-track error cannot be perpendicular to the path. For this reason, the projected ground speed of the actual aircraft (
${V_{rs}}$
) through the path is calculated with Equation (29).
${V_{gs}}$
is the ground speed of the actual aircraft.
Namely, the aircraft’s ground speed along the path is denoted by
${V_{rs}}$
. Wherever the airplane is flying,
${V_{rs}}$
is dependent on the path angle
${\alpha _k}$
. Specifically, the actual aircraft’s groundspeed is divided into two vectors that are perpendicular to one another. One is perpendicular to the path, while the other (
${V_{rs}}$
) is parallel to it and it is shown in Fig. 8
The ground speed (
${{\rm{V}}_{{\rm{gs}}}}$
) and the projected ground speed (
${{\rm{V}}_{{\rm{rs}}}}$
) of the aircraft.

Now, we can use lookahead-based steering to generate the desired course angle (
${\chi _{ref}}$
):
The given equation can be used to calculate the velocity path relative angle (
${\chi _r}\left( e \right)$
):
Lookahead distance (
$\Delta $
) is shown in Fig. 7. The value of the lookahead distance increases or decreases the convergence because it is a guidance parameter [Reference Breivik and Fossen32]. Namely, this value has an effect like a proportional gain, and we can rewrite the equation again:
${\chi _r}\left( e \right)\,:\!= arctan\left( { - {K_p}e} \right)$
. Lookahead distance values can be taken as constant values. However, the same value can be a problem if the cross-track error is too large or too small. For this reason, the equation given in [Reference Lekkas and Fossen23] is used:
For lookahead distance, the selected maximum value is
${\Delta _{max}}$
, and the minimum is
${\Delta _{min}}$
. This equation is useful for forcing the aircraft to follow the predefined path aggressively when too far, whereas gently when too close. That is to say,
$\Delta \left( {e\left( t \right)} \right)$
value changes from
${\Delta _{min}}$
to
${\Delta _{max}}$
depending on the aircraft’s cross-track error.
In Figs. 9 and 10, the results for different
${\Delta _{min}}$
and
${\Delta _{max}}$
values are shown. After 45 degrees heading change shown in Fig. 9, aircraft is established on-course faster with the values of
${\Delta _{max}} = 250$
. However, after 90 degrees heading change shown in Fig. 10, aircraft overshoots the path with the value of
${\Delta _{max}} = 250$
. If
${\Delta _{max}} = 500$
is chosen, the aircraft is established on-course smoothly. As a result, logic has to be followed for
${\Delta _{min}}$
and
${\Delta _{max}}$
values to improve lateral navigation.
Lateral navigation regarding different
${{{\Delta }}_{{\rm{min}}}}$
and
${{{\Delta }}_{{\rm{max}}}}$
values after waypoint-2.

Lateral navigation regarding different
${{{\Delta }}_{{\rm{min}}}}$
and
${{{\Delta }}_{{\rm{max}}}}$
values after waypoint-3.

In Fig. 11, the right aileron deflections are shown during lateral navigation with respect to different
${\Delta _{min}}$
and
${\Delta _{max}}$
values. These right aileron deflections are for the lateral navigation shown in Figs. 9 and 10. There are small changes in deflections.
Right aileron actuator position during lateral navigation regarding different
${{{\Delta }}_{{\rm{min}}}}$
and
${{{\Delta }}_{{\rm{max}}}}$
values.

In Fig. 12, the whole lateral navigation is shown with respect to
${\Delta _{max}} = 250$
and
${\Delta _{max}} = 500$
.
Lateral navigation regarding different
${{{\Delta }}_{{\rm{min}}}}$
and
${{{\Delta }}_{{\rm{max}}}}$
values.

We include another logic for minimising cross-track error. First, the angle (
$\lambda $
) between path tangential angle (
${\alpha _k}$
) and aircraft heading (
${\psi _{sensed}}$
) is given in Equation (33) shown in Fig. 13. If the absolute value of
$\lambda $
is between
$0$
and
$50$
degrees,
${\Delta _{min}} = 20$
and
${\Delta _{max}} = 250$
. However, if the absolute value of
$\lambda $
is between
$50$
and
$120$
degrees,
${\Delta _{min}} = 20$
and
${\Delta _{max}} = 500$
. The result of the lateral navigation with
$\lambda $
logic is shown in Fig. 14.
${\rm{\lambda }}$
logic, the ground speed and the projected ground speed of the aircraft.

Lateral navigation regarding
${\rm{\lambda }}$
logic.

Besides,
${\gamma _c}$
in the Equation (32) is the convergence value, which is a design parameter. For NASA GTM model,
${\gamma _c} = 0.01$
. With the help of Equation (32),
$\Delta \left( {e\left( t \right)} \right)$
can be computed and fed to the Equation (31). The desired course angle is
${\chi _{ref}}$
and given in Equation (30). For computing
${\psi _{ref}}$
, the equation is given:
In Fig. 15, there is a lateral activation logic. It decides on which mode path will be used. For example, if LNAV/NAV is adjusted from the mode control block, it activates the waypoint database path for generating
${x_{ref}}$
and
${y_{ref}}$
which is shown in Fig. 15. There are two types of waypoints which are fly-by and fly-over. In fly-over waypoints, the next leg is activated when the aircraft is over the waypoint. In fly-by waypoints, there is turn anticipation, and the aircraft begins the turn before reaching the waypoint, smoothly establishing the next leg. In our work, fly-over waypoints are used. However, in [Reference Ergöçmen and Tilki33], fly-by waypoints are employed and the switching mechanism is given in detail.
Roll channel-lateral guidance.

HDG (Heading): Also, other than LNAV/NAV mode, heading mode can be activated and adjusted. After that, the value that is adjusted from the mode control block is transmitted to the lateral guidance shown in Fig. 15. In this block, this signal enters as a manual
${\psi _{ref}}$
and directly goes to the lateral activation switch. If HDG mode is pushed, lateral activation logic activates this path and sends a signal to the lateral activation switch. Finally, the adjusted heading value (manual
${\psi _{ref}}$
) can pass.
Logic for approach (lateral): While descending and approaching the runway, the lateral navigation is vital. Instead of the path following before approaching, heading can be adjusted and the air traffic controller (ATC) can vector the aircraft to the final approach course. For vectoring, HDG HOLD mode has to be adjusted from the mode control block. To make an approach, the APP has to be activated. After pushing the APP button from the mode control block, the aircraft continues to fly in HDG HOLD mode until final approach fix (FAF) is reached. FAF is the point where the aircraft starts the approach. If the aircraft arrives at the FAF, the mode automatically changes to LNAV mode for lateral path following (LOS guidance law) with the help of the lateral activation logic.
3.2.2 Pitch channel-vertical navigation
Vertical guidance is shown in Fig. 16. There are PIDs for generating flight path angles (
${\gamma _{ref}}$
). However, the selection of the block for generating
${\gamma _{ref}}$
depends on the mode control block, which activates the vertical activation logic.
Pitch channel-vertical guidance.

Vertical navigation (alt): If ALT HOLD is selected and the value (
${h_{ref}}$
) is adjusted, then
${h_{ref}}$
value is transmitted to the vertical navigation (alt) block shown in Fig. 16 and
${h_{sensed}}$
value is subtracted from this value. After that, with the help of the PID, (
${\gamma _{ref}}$
) is generated. If the flight path angle is manually adjusted from the mode control block, it activates the vertical activation logic and switch again. Finally, the manual (
${\gamma _{ref}}$
) value goes directly.
Vertical navigation (GS): This block is used during the descent to the runway. Other than the approach, this VNAV mode can be used for descending from an enroute flight level. However, in this work, it is mainly used for approaching the runway. As stated for the roll channel, APP is not activated until arrival to FAF, even if the APP button is pushed. The same idea for the pitch channel, APP is in the armed position before FAF, and ALT HOLD is currently activated on this flight. If the aircraft arrives at FAF and also APP button is pushed before, with the help of the vertical activation logic, the aircraft’s vertical navigation changes from the vertical navigation (alt) block to the vertical navigation (GS) block. After that, the aircraft starts to descent with minimising deviation error (
${\rm{\Gamma }}$
) from the glide slope as shown in Fig. 17. Basically, for the calculation of
${\rm{\Gamma }}$
, this formulation can be used:
${\rm{\Gamma }} \simeq 57.3d/R$
in degrees. Distance from the glide slope (
$d$
) is shown in Fig. 17. However, it is used mainly for glide slope, which uses a beam, for example, the instrument landing system (ILS). In this work, a deviation from the glide slope (
${\rm{\Gamma }}$
) is calculated differently.
${\rm{\Gamma }}$
is found by this equation:
Glide slope.

The angle between the beginning of the runway threshold (
${x_0}$
,
${h_0}$
) and the aircraft (
${x_a}$
,
${h_a}$
) is
${\gamma _{deviate}}$
and is given as:
It is assumed that, the height (
${h_{true}}$
) between the aircraft and the runway elevation can be measured by radar altimeter. The distance (
${x_{true}}$
) between the aircraft and the runway threshold by GPS can be measured easily shown in Fig. 17.
For the glide slope, the angle is
${\gamma _{path}}$
(the value is 3 degrees). The equation is given as:
Finally, we can find
${\rm{\Gamma }}$
and minimise this value with the PID algorithm. As a result, the aircraft approaches the runway with the predefined glide slope angle (
${\gamma _{path}}$
). After passing through PID,
${\gamma _{ref}}$
is generated shown in Fig. 16.
Flare: The landing can be subdivided into four phases (glide slope tracking, flare, touchdown and landing roll). Glideslope tracking is introduced in vertical navigation. During a flare, the pilot performs an exponential curve, and autoland system does the same:
${h_{flare}}$
(ft) is the height for the flare initiation.
$\tau $
is a time constant. After differentiating Equation (38):
At the initiation of the flare (
${h_{flare}}$
,
$t = 0$
), the rate of descent for 3 degrees glideslope is given:
${U_0}$
is the velocity (ft/sec) of the aircraft, and it is assumed to constant during flare path. From Equation (40), if the rate of descent is taken
$ - 5.4188{\rm{\;\;}}ft/sec.$
during the initiation of the flare, it becomes:
\begin{align}{\dot h_{flare}} = - \frac{{{h_{flare}}}}{\tau }{\rm{\;\;}}& \Rightarrow {\rm{\;\;\;\;}} - 5.4188 = - \frac{{{h_{flare}}}}{\tau }\nonumber\\&\Rightarrow {\rm{\;\;\;\;}}{h_{flare}} = \left( {5.4188} \right)\left( \tau \right)\end{align}
This rate of descent (
$\dot h$
) has to be less than
$ - 2{\rm{\;\;}}ft/sec.$
during touchdown. As a result, we have to control the rate of descent. The flare control block diagram is shown in Fig. 18.
${h_{sensed}}$
is multiplied with
$ - 1/\tau $
and becomes
${\dot h_{ref}}$
. It is subtracted from
${\dot h_{sensed}}$
. Error dynamic is
$e = \left( {\left( { - 1/\tau } \right){h_{sensed}}} \right) - {\dot h_{sensed}}$
. After that, this error is multiplied by the gain and transfer function [Reference Blakelock34]. Finally, it becomes
${\gamma _{ref}}$
for the flight path controller.
Flare control.

In the flare control problem,
$\tau $
and
${h_{flare}}$
have to be calculated.
$\tau $
can be calculated from the distance to the touchdown point (
$D$
) shown in Fig. 18. The value of
$D$
is
$600ft$
. It is assumed that flare is finished and touchdown occurs approximately in 4 time constant (4
$\tau $
). The distance of the flare path is given:
In Equation (42), there are two unknowns (
$\tau $
and
$R$
). From the trigonometric relation:
After substituting Equation (43) in Equation (42), we obtain
$\tau = 1.93$
and
$R = 200ft$
and
${h_{flare}} = 10.5ft$
. The duration of the flare is
$4\tau $
and it equals to 7.72 sec.
NASA GTM Model descent with −0.32 ft/min along the 3 degrees glideslope. The flare manoeuver begins at 510.5 ft (+10.5 ft from threshold height) shown in Fig. 19. The beginning of the runway is 7,744 meters from the Waypoint-1 (0; 0). The altitude of the runway is 500 ft. However, the touchdown point is approximately 7,850 meters. At this point, vertical descent rate is −0.03 ft/min (−0.016ft/sec), which is acceptable for a smooth landing.
Flare manoeuver.

In Fig. 20, the whole process for the mode control panel and path following and approach controller (roll and pitch channels) are summarised. There are three inputs for the roll/pitch channels and autothrottle. The autothrottle is a separate one. The most complicated one is APP mode. After pushing APP, also if one of the modes for roll (NAV/LNAV or HDG), and one of the modes for pitch (ALT, FPA or VNAV) are activated, the logic becomes true (1). However, there is a second check for whether the aircraft is at the final approach fix (FAF). If it is not at the FAF, the previous mode continues and APP is not activated. If the aircraft is at the FAF, the approach logic is activated for the pitch and roll. NASA GTM tracks the path from FAF to the runway laterally and vertically. For roll channel-lateral navigation, the LOS algorithm generates reference heading
${\psi _{ref}}$
by comparing the actual position (
${x_{sensed}},{y_{sensed}}$
) and the runway position (
${x_0},{y_0}$
). For this computation, the related equations in Section 3.2.1 roll channel-lateral navigation are used. For the pitch channel-vertical navigation, the reference path angle (
${\gamma _{ref}}$
) is generated with the help of the deviate path angle (
${\gamma _{deviate}}$
) given in Equation (36) and the runway path angle (
${\gamma _{path}}$
) given in Equation (37). After
${\gamma _{ref}}$
, if the aircraft altitude reaches the
${h_{flare}}$
altitude, the flare mechanism is activated. All equations for the pitch channel are given in Section 3.2.2 pitch channel-vertical navigation.
Finite state machine: mode control panel with pitch/roll channel of the path following approach controller.

Finally, there are three outputs, which are reference heading (
${\psi _{ref}}$
), reference path angle (
${\gamma _{ref}}$
), and reference throttle (
${\delta _t}$
) from the state flow shown in Fig. 20.
3.3 Flight path controller
After path following and approach controller,
${\psi _{ref}}$
and
${\gamma _{ref}}$
are transmitted to the flight path controller shown in Fig. 21.
${\psi _{ref}}$
is generated by the roll channel and
${\gamma _{ref}}$
by pitch channel before. These values are subtracted from respectively
${\psi _{sensed}}$
and
${\gamma _{sensed}}$
values, which are measured by sensors. After that, these errors pass through PID controllers. Finally,
${\mu _{ref}}$
(reference bank angle) and
${\alpha _{ref}}$
(reference angle of attack/alpha) are generated, respectively.
${\mu _{ref}}$
is limited with
$ + / - {50^ \circ }$
and
${\alpha _{ref}}$
with
$ + / - {20^ \circ }$
. Also,
${\beta _{ref}}$
is adjusted to zero due to acquiring a coordinated turn during the heading change. It is possible to compute this reference side slip angle (
${\beta _{ref}}$
) in other complex yaw damper systems.
Flight path control.

3.4 Attitude controller
Flight path control loop sends reference values (
${\mu _{ref}}$
,
${\alpha _{ref}}$
and
${\beta _{ref}}$
) to attitude controller shown in Fig. 22. These values are subtracted from respectively
${\mu _{sensed}}$
,
${\alpha _{sensed}}$
and
${\beta _{sensed}}$
values, which are measured by sensors. Then, they pass through individual 3 PID controllers to generate virtual controls (
${\dot \mu _{ref}}$
,
${\dot \alpha _{ref}}$
and
${\dot \beta _{ref}}$
). These values are used in the NDI function block (shown in Fig. 22). There is no need to use aerodynamic derivatives (aircraft model) for generating reference rate values (
${p_{ref}}$
,
${q_{ref}}$
,
${r_{ref}}$
) in this block. In other words, dynamic inversion is achieved, but the aircraft aerodynamic derivatives are not required. The formulation of the function [Reference Lu, van Kampen, de Visser and Chu11] is given as:
\begin{align}\left[ {\begin{array}{*{20}{l}}{{p_{ref}}}\\{{q_{ref}}}\\{{r_{ref}}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}{{\rm{cos}}\alpha {\rm{cos}}\beta } {}&0 {}&{{\rm{sin}}\alpha }\\{{\rm{sin}}\beta } {}&1 {}&0\\{{\rm{sin}}\alpha {\rm{cos}}\beta } {}&0 {}&{ - {\rm{cos}}\alpha }\end{array}} \right]\left[ {\begin{array}{*{20}{l}}{{{\dot \mu }_{ref}}}\\{{{\dot \alpha }_{ref}}}\\{{{\dot \beta }_{ref}}}\end{array}} \right] + {T_{BV}}\left[ {\begin{array}{*{20}{c}}{ - \dot \chi {\rm{sin}}\gamma }\\{\dot \gamma }\\{\dot \chi {\rm{cos}}\gamma }\end{array}} \right]\end{align}
Attitude control.

There is a velocity to body frame transformation matrix (
${T_{BV}}$
) and
${T_{BV}}$
is the transpose of the body to velocity frame transformation matrix (
${T_{VB}}$
). Transformation matrix (
${T_{BV}} = T_{VB}^T$
) is given:
\begin{align}{\left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}{{cos\alpha {\rm{\;\;\;\;}}cos\beta }} {}&{{\rm{\;\;\;\;}}sin\beta {\rm{\;\;\;\;}}} {}&{{sin\alpha {\rm{\;\;\;\;}}cos\beta }}\\{r21} {}&{\left( {cos\beta {\rm{\;\;\;\;}}cos\mu } \right)} {}&{r23}\\{r31} {}&{\left( {cos\beta {\rm{\;\;\;\;}}sin\mu } \right)} {}&{r33}\end{array}} \right]^T}\end{align}
\begin{align*} r{21}& = \left( { - cos{\rm{\;\;\;\;}}sin\beta {\rm{\;\;\;\;}}cos\mu + sin\alpha {\rm{\;\;\;\;}}sin\mu } \right) \\ r{23} &= \left( { - sin\alpha {\rm{\;\;\;\;}}sin\beta {\rm{\;\;\;\;}}cos\mu - cos\alpha {\rm{\;\;\;\;}}sin\mu } \right)\\r{31} &= \left( { - cos\alpha {\rm{\;\;\;\;}}sin\beta {\rm{\;\;\;\;}}sin\mu - sin\alpha {\rm{\;\;\;\;}}cos\mu } \right)\\r{33} &= \left( { - sin\alpha {\rm{\;\;\;\;}}sin\beta {\rm{\;\;\;\;}}sin\mu + cos\alpha {\rm{\;\;\;\;}}cos\mu } \right) \end{align*}
Other unknowns are derivatives of azimuth angle (
$\dot \chi $
) and flight path angle (
$\dot \gamma $
), which are given in Equations (46), (47). Also, in Equations (46), (47), there are accelerations along the body axis (
${A_x}$
,
${A_y}$
and
${A_z}$
). They are measured directly by accelerometers.
$V$
is the total velocity of the aircraft. Bank angle (
$\mu $
) and alpha (
$\alpha $
) are measured directly.
3.5 Rate controller
After attitude controller, the signal (
${p_{ref}}$
,
${q_{ref}}$
,
${r_{ref}}$
) reaches to rate controller shown in Fig. 23. It is subtracted from the measured signal (
${p_{sensed}}$
,
${q_{sensed}}$
,
${r_{sensed}}$
), and then passes through the individual three PIDs. Virtual controls (
$v = [{\dot p_{ref}},$
${\dot q_{ref}},{\dot r_{ref}}{]^T}$
) are calculated by this process. After that, this signal reaches the INDI block. In this part, angular accelerations in the last time step (
${\dot{x}_{0}} = [\dot{p}_{{0}_{sensed}},\dot{q}_{{0}_{sensed}},\dot{r}_{{0}_{sensed}}]^{T}$
) are calculated by the filter given in Equation (8). As stated before, with the help of the virtual controls and angular accelerations, the aircraft’s main dynamics, and aerodynamic derivatives are not required. By subtracting virtual controls from angular accelerations and multiplied with the inertia matrix, incremental moments (
${{\Delta }}{L_{control}}$
,
${{\Delta }}{M_{control}}$
,
${{\Delta }}{N_{control}}$
) is calculated for control allocation (CA).
Rate control.

We can simplify and write the formulation without aircraft dynamics. Instead of
$A\left( x \right)$
matrix given in [Reference Ergöçmen and Yavrucuk35], which contains aircraft aerodynamic derivatives (
${L_p}$
,
${N_p}$
,
${L_r}$
,
${N_r}$
,
${M_q}$
,
${L_p}$
,
${N_p}$
,
${L_r}$
,
${N_r}$
) and states, angular accelerations are used. The formulation [Reference Sieberling, Chu and Mulder4] for the INDI is given in Equation (48):
\begin{align}{\left[ {\begin{array}{*{20}{l}}{\dot p}\\{\dot q}\\{\dot r}\end{array}} \right]_{ref}} = {\left[ {\begin{array}{*{20}{l}}{{{\dot p}_0}}\\{{{\dot q}_0}}\\{{{\dot r}_0}}\end{array}} \right]_{sensed}} + {\left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}{{I_{xx}}} {}&0 {}&{ - {I_{xz}}}\\0 {}&{{I_{yy}}} {}&0\\{ - {I_{xz}}} {}&0 {}&{{I_{zz}}}\end{array}} \right]^{ - 1}}\left( {{B_c}} \right){\rm{\;\;\;\;\Delta }}u\end{align}
Before the signal goes to CA block, virtual input (desired moments) becomes
${v_{des}} = {[{{\Delta }}{L_{control}}{\rm{\;\Delta }}{M_{control}}{\rm{\;\Delta }}{N_{control}}]^T}$
[Reference Lombaerts, Kaneshige, Schuet, Aponso, Shish and Hardy12] shown in Fig. 23:
\begin{align}\left[ {\begin{array}{*{20}{l}}{{\rm{\;\;\Delta }}{L_{control}}}\\{{\rm{\;\;\Delta }}{M_{control}}}\\{{\rm{\;\;\Delta }}{N_{control}}}\end{array}} \right] = \left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}{{I_{xx}}} {}&0 {}&{ - {I_{xz}}}\\0& {}{{I_{yy}}} {}&0\\{ - {I_{xz}}} {}&0 {}&{{I_{zz}}}\end{array}} \right]\left[ {\begin{array}{*{20}{l}}{{{\dot p}_{ref}} - \dot{p}_{{0}_{sensed}}}\\{{{\dot q}_{ref}} - \dot{q}_{{0}_{sensed}}}\\{{{\dot r}_{ref}} - \dot{r}_{{0}_{sensed}}}\end{array}} \right]\end{align}
Finally, increment input (
${{\Delta }}{u_{ref}}$
) shown in Fig. 23 is given in Equation (51).
${B_c}$
is not a square matrix, it is a 3 by 8 matrix. For this reason, dynamic control allocation is used.
\begin{align} &= B_c^{ - 1}\left[ {\begin{array}{*{20}{c@{\quad}c@{\quad}c}}{{I_{xx}}} {}&0 {}&{ - {I_{xz}}}\\0 {}&{{I_{yy}}} {}&0\\{ - {I_{xz}}} {}&0 {}&{{I_{zz}}}\end{array}} \right]\left[ {\begin{array}{*{20}{l}}{{{\dot p}_{ref}} - \dot{p}_{{0}_{sensed}}}\\{{{\dot q}_{ref}} - \dot{q}_{{0}_{sensed}}}\\{{{\dot r}_{ref}} - \dot{r}_{{0}_{sensed}}}\end{array}} \right] \\[8pt]\nonumber\end{align}
In this equation, the incremental control command (
${{\Delta }}u$
) is calculated.
${u_0}$
is the input at the last time step. However, we need the total:
Finally, the control command (
$u$
) is computed and ready to be sent to the actuator model. This entire process is summarised in Fig. 23.
4.0 Results
4.1 Case-1: Control surfaces lock-in-place
Lateral guidance is shown in Fig. 24. There are five fly-over waypoints. It means the NASA GTM model has to fly over it. There is no turn anticipation.
Lateral guidance regarding lock-in-place ailerons and elevators.

NASA GTM starts the first leg between waypoint-1 (0; 0) and waypoint-2 (1,000; 0). During flying this leg, a 15-degree lock-in-place is simulated for the right or left elevator and the right or left aileron. After waypoint-3 (initial approach fix), (3,000; 2,000), it changes heading from 45 degrees to 135 degrees for establishing this leg. Waypoint-4 (4,000; 1,000) is an intermediate fix, and it is used for level flight before the approach. Waypoint-5 (6,000; 1,000) is a final approach fix, and NASA GTM begins its descent and approaches the runway with a 3 degree glide slope. 40 meters left and 33 meters right deviation with respect to
$X$
distance occurs for a 15-degree left or right aileron lock-in-place, 10 meters left/right deviation with respect to
$X$
distance occurs for 15 degrees left or right elevator lock-in-place shown in Fig. 25. It is an acceptable result. Whenever the mentioned deviations occur, NASA GTM can re-establish the route every time.
Recovery from control surface lock-in-place.

During aileron lock-in-place, if the LOS has an adaptive manner concerning angle (
$\lambda $
) given in Equation (33), the cross-track error is minimised fast shown with the blue line in Fig. 26. After lock-in-place occurs, the aircraft deviates from the path approximately 40 meters with adaptive LOS (with
$\lambda $
logic). After that, compared with the brown line which is the result of the LOS without an adaptive manner (without
$\lambda $
logic), the aircraft establishes the path again quickly before reaching the waypoint. Namely, without an adaptive feature, LOS algorithm regulates cross-track error ineffectively compared with the adaptive LOS case shown in Fig. 26.
Recovery from 15 degrees right aileron stuck. Blue line: with adaptive LOS (with
${\rm{\lambda }}$
logic). Brown line: without adaptive LOS (without
${\rm{\lambda }}$
logic).

For Required Navigation Performance (RNP) approaches (also known as GPS approaches), the GPS receiver should have a receiver autonomous integrity monitoring (RAIM) capability. It monitors errors during navigation and if the aircraft deviates from a specific limit from the path, the RAIM alerts the pilot. This total system error (TSE) is the sum of these errors:
-
• Path definition error (PDE): Errors caused by incorrect or ambiguous representation of the desired flight path in the navigation database.
-
• Navigation system error (NSE): Errors in the aircraft’s navigation system, including inaccuracies in the position estimation from GPS, inertial systems or other sensors.
-
• Flight technical error (FTE): Errors caused by the pilot or autopilot in maintaining the aircraft’s position along the defined flight path.
For this TSE, limits are 2 nm (3.7 km) left and 2nm right of the path for the enroute phase, 1 nm (1.9 km) for the terminal phase, and 0.3 nm (0.6 km) for the approach and it is getting narrow to the runway threshold [36]. The lock-in-place occurs at the terminal phase. The maximum deviation is 40 meters shown in Fig. 26. The limit for the terminal phase is 1.9 km. As a result, if we compare the deviation with the limit, the autopilot can fly the aircraft with a high level of safety.
In Fig. 27 ail stuck15, a 15-degree left aileron lock-in-place occurs approximately at the 4th second of the simulation. The right aileron compensates undesired roll due to an emergency. After that, NASA GTM Model follows the path with a 25 and −40 degrees roll angle (
$\phi $
) maximum.
$\beta $
angle is 1 degree which is an acceptable value. The right aileron, the left and right elevators, and rudders are not saturated.
15 degrees left aileron lock-in-place.

In Fig. 28, a 15-degree right aileron lock-in-place occurs. In this case, the heading (
$\psi $
) oscillates slightly at the beginning of the heading (
$\psi $
) change. However, during a level flight after a heading change, there is no oscillation. The left aileron, the left and right elevators, and rudders are not saturated.
15 degrees right aileron lock-in-place.

In Figs. 29 and 30, right (inboard and outboard) and left elevator (inboard and outboard) respectively lock-in-place occurs at the 4th second. Small deviation roll moment due to emergency does not induce too much heading change. As a result, path following for lateral navigation can be done easily. The right aileron, the left aileron and rudders are not saturated.
15 degrees right elevator lock-in-place.

15 degrees left elevator lock-in-place.

4.2 Case-2: Left or right engine out
Left or right engine out is studied, and lateral and vertical guidance for them are shown in Fig. 31. At the 4th second of the simulation, left or right engine out occurs. This point is between waypoint-1 (0; 0) and waypoint-2 (1,000; 0). However, the NASA GTM model can fly safely and land on the runway easily (7,744; 0). The proposed controller without reconfiguration is robust enough for these types of emergencies as well as control surfaces lock-in-place.
3D Guidance for one engine out without wind and turbulence.

4.3 Case-3: Moderate turbulence and wind
In this case, the controller is put to the test in windy and turbulent conditions. As a turbulence model, NASA-TM-4,511 and MIL-STD-1,797 Appendix-A are used. Moderate turbulence and 10 kt (16.87 ft/sec) at 70 degrees wind are shown in Fig. 32. The NASA GTM flies at 75.88 kt (128.07 ft/sec) and after 1 second it decreases the airspeed to 65 kt (109.70 ft/sec) for an approach. Moderate turbulence affects the aircraft from 3 axes between −9 kt (−15 ft/sec) to +9 kt (15 ft/sec). If the aircraft is exposed to this turbulence and wind, its airspeed oscillates shown in Fig. 33. However, the controller compensates for these undesired effects. Flight parameters and control surface positions are shown in Fig. 34. The elevators command is between 7 and −3 degrees, ailerons command is between −15 and 15 degrees, rudder command is between −10 and 10 degrees.
10 kt at 70 degrees wind and moderate turbulence.

Aircraft airspeed during 10 kt at 70 degrees wind and moderate turbulence.

Flight parameters during wind and turbulence.

3D Guidance at 800 ft altitude with wind (10 kt and from 70 degrees) and moderate turbulence are shown in Fig. 35. In a different case, 3D guidance with 10 kt at 100 degrees wind and moderate turbulence are shown in Fig. 36.
3D guidance with 10 kt at 70 degrees wind and moderate turbulence.

3D guidance with 10 kt at 100 degrees wind and moderate turbulence.

As a result, the controller is robust enough for turbulence and wind effects.
5.0 Conclusion
In this work, the highly nonlinear NASA GTM design is used. In terms of controllers, the full autopilot design combines PID, INDI, NDI and CA algorithms. Path following is executed using the adaptive LOS algorithm. The lookahead distance’s minimum and maximum values changes concerning the angle between the aircraft and the path. Virtual aircraft airspeed (
${V_{rs}}$
) is calculated from the actual aircraft for the cross-track error calculation. It is a different approach.
The parameters of the LOS algorithm are adaptive. The logic (
$\lambda $
) behind it modifies this value in relation to the error between the path angle and the aircraft heading. In addition to enabling quick path entry, adaptive LOS also provides fault tolerance in the event of engine-out issues or control surface damage.
A mode control panel is used for realistic simulations. The finite state machine diagram for the mode control panel is provided. In this work, the whole autopilot framework has been studied extensively. Detailed instructions are given for it. The designed controller has fault-tolerant capabilities. To evaluate the controller’s robustness, control surface issues for various control surfaces lock in place, and engine failures are studied. The terminal, approach, and flare phases are simulated. Additionally, the controller can easily compensate for undesired effects, such as moderate turbulence and wind.
The autopilot design can be changed to a completely autonomous one for future work. Instead of using database waypoints for navigation, a rapidly exploring random tree (RRT) can be used for generating waypoints. The autopilot’s mode switches to heading mode automatically in order to avoid the threat. Also, an emergency flight planner (EFP) can be developed for different kinds of emergencies. This system is directly applicable to autonomous flight, incorporating trajectory generation, path planning and path following.
Competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.















