Runway Safety Assistant Foreseeing Excursions: Calculating Means

1. Introduction

The most critical stages during the flight are the Take off and Landing, as extremely large changes in the kinetic energy of the aircraft happen during that time. The constraint that the length of the runway sets, as well as the existing environmental conditions make the task even more difficult for the operators, where they have to cope with them, occasionally without alternatives. One of the types of accidents that occur in these stages leads to exceeding the length of the runway, a type of accident with low mortality rate. However it leads to passenger confusion, temporary unavailability of the runway, damage to the aircraft and/or the runway itself, defamation of the airline, and more. Runway overrun incidents in civil aviation, occupied the first place in runway related incidents during the 2010-2014 period [1].

Since the human factor is the leading cause of aviation accidents [2], there should be a system that supports and warns the operators in cases of a potential runway overrun. That way, the risk and occurrences of such incidents can be decreased.

1.1. Calculation of Critical Speeds

Before starting any Take off or Landing procedure and based on a number of factors, some critical speeds (V1,V2 and VR for the Take off, Vapp and Vref for Landing) are derived. Their calculation is done through the flight manuals of the respective aircraft and recently through software solutions (e.g. Flysmart+ [3]). The first ones produce approximate calculations with large margins of error, while software calculations are more precise, including safety margins. The scope of these speed indicators is to complete each process with safety and prevent the airplane from overshooting the runway. At the same time, acceleration and deceleration rates that occur, prevent the straining of the aircraft’s mechanisms and provide comfort to the passengers.

1.2. Motivation

The initiative of this system’s development is EASA’s (European Union Aviation Safety Agency) research agenda for 2022-2024 . Reference LOC-03, named "Landing and Take off monitor" requests to "assess means to assist the flight crew in:

- preventing runway overrun and managing aircraft total energy

- monitoring the actual acceleration of the aircraft during the Take off run to detect mismatch between V1 and the actual remaining runway distance" [4]

Statistics of accidents in recent years, justify the inclusion of such a research interest from EASA.

1.3. Assumptions

The statically calculated critical speeds assume the smooth development of the procedures. In particular, it is assumed that during Take off the engines perform at 100%. Also, during Landing, the aircraft comes in contact with the ground before a certain point, and that the friction of the wheels with the runway is exactly as predicted. Although major deviation factors are easily observable by operators (e.g. engine loss, flat tire), others are not that obvious (e.g. reduced engine performance due to need of regular maintenance, worn tires), while in other cases operators misjudge the capabilities of the aircraft.

The is due to the inherent inability in predicting anomalies during the evolution of the process, such as those mentioned above. This also leads to the second cause, which is that values extracted from the static calculations correspond to velocities, while the main limiting factor is the length of the runway. Even though, the aircraft’s instruments indicate the aircraft’s speed accurately, there is no indication of its position relative to the runway, for which the crew can derive conclusions only by the visual contact they have with it. This is where the experience of the pilots must contribute to assess whether the remaining runway is sufficient to bring the aircraft to a full stop. However, such a thing is practically impossible to calculate precisely.

1.4. System Specifications

Based on the above, it is understandable that there is a need for an auxiliary decision support system that belongs to the aircraft’s cockpit and is able to calculate the point at which critical speeds are attained on the runway, and the most important being the point at which the aircraft is expected to come to a full stop. The system must be able to perform two different calculations, a static and dynamic one.

The static calculation is based excursively on user’s inputs and aims to calculate the distance needed in each procedure. For Take off the Accelerate Stop Distance Required (ASD) is calculated, while during Landing the corresponding distance is Landing Distance Required (LDR). Such calculations cannot take into account the real world interference, thus no matter how precise the calculation is, it is not able to predict unexpected events and/or the actual deviation. This is why real time/online algorithms are a necessity.

The dynamic calculation starts when the acceleration of the aircraft is initiated during Take off and at a certain height in the case of Landing. The calculations are based solely on factors other than user inputs, such as acceleration, velocity, height from the ground etc.

1.5. Obtaining Necessary Data

All the static calculations mentioned previously, require data regarding the aerodynamic characteristics of the aircraft, the performance of the engines and the interaction with the surrounding environment. Little to no data is available online concerning the aircraft and even if there are, data verification is required. For this purpose, the flight simulator X-Plane 11 is used, which provides both the necessary data and the verification framework for the correct operation of the system. This particular simulator is chosen based on its ability to represent the real world as closely as possible to the laws of physics.

To interact and retrieve these values from the simulator, suitable code has been developed in Python. The integrated UDP, sends the user-selected information in packets. These packets are sent on every new calculation the simulator makes, i.e. every new frame. The data is either processed or stored, depending on the use for which it is intended. The data fields sampled from the simulator can be found analytically at the end of the document.

2. State of the art

The specifications of a Take Off Performance Monitoring System (TOPMS) were defined for the first time in August of 1987, in the AS-8044 standard (Aerospace Standard 8044), while they were modified in 2007 and 2020 [5]. Three types of such systems are defined, depending on their function. These are:

• Type 1: non-predictive systems, only compare actual to expected performance
• Type 2: prediction the Take off distance
• Type 3: prediction Take off distance and braking distance in case of abort

Although a TOPMS prototype has been developed by NASA [6] and other organizations, there is no system widely used in aviation. EASA took the initiative to set up a study group (Working Group 88), exclusively for the study of TOPMS , which in 2017 concluded that such systems are unfeasible due to limitations in technology and data [7].

A recent proposal for a TOPMS came from Netherlands Aerospace Centre (NLR), in the Take off Performance Alerting Program (TOPAP) is proposed. It only detects and notifies for major errors in the expected and actual acceleration during Take off, which signifies difference in the actual and given weight, during the necessary pre-Take off calculations and the weight input at the Flight Management System (FMS) [8].

On the other hand, runway excursion alerting systems for Landing are being developed and used for more than a decade in commercial aviation. Since 2020, all new aircrafts must be equipped with a Runway Overrun Awareness and Alerting System (ROAAS), in order to comply with EASA’s CS-25 Amendment 24 (Certification Specifications for Large Aeroplanes). The two largest aircraft manufacturers, Airbus and Boeing, have developed their own solutions named Runway Overrun Prevention System (ROPS) and Runway Awareness and Advisory System (RAAS) respectively.

EASA explicitly defines that "The ROAAS shall reduce the risk of a longitudinal runway excursion during Landing by providing alert, in flight and on ground, to the flight crew when the aeroplane is at risk of not being able to stop within the available distance to the end of the runway.

1

During approach (from a given height above the selected runway) and Landing, the ROAAS shall perform real-time energy-based calculations of the predicted Landing stopping point, compare that point with the location of the end of the runway, and provide the flight crew with:

(a)

in-flight, timely, and unambiguous predictive alert(s) of a runway overrun risk, and

(b)

on-ground, timely, and unambiguous predictive alert(s) of a runway overrun risk. At the option of the applicant, the ROAAS may also provide an automated means of deceleration control that prevents or minimises runway overrun during Landing.

2

The ROAAS shall at least accommodate dry and wet runway conditions for normal Landing configurations." [9]

The Landing mode of the RUN.S.A.F.E. system, was developed in compliance with the ROAAS definition, apart from Item 1.b of the definition, hence new aircrafts can be equipped with the system, or it can be installed in older aircrafts in order to be up to date with the newer regulations. A full adaption to an aircraft will also be able to satisfy Item 1.b, accessing control to automated means of deceleration.

3. Take Off Procedure

Take off is defined as the phase during which the aircraft aims to abtain a certain speed that sets it capable of flying under its own lift. However, there are cases where the aircraft encounters problems during this process, leading to potential overruns, if handled with a rejected Take off (RTO). Factors that can lead to these situations include technical problems with the aircraft’s systems, adverse weather conditions, deficiencies in the execution of Take off procedures by the crew and factors affecting runway availability.

In the case of a RTO, the crew must follow certain procedures that aim to stop the aircraft with safety. Though, often times the overestimation of the crew’s and the aircraft’s abilities can lead to a runway overrun.

Currently, the standard term that prevents a runway overrun during a take- off is the V1 speed (also refereed to as decision speed). This critical speed marks the point beyond which the pilot no longer has the capability to cancel the Take off procedure, as the remaining length of the airport is not sufficient to safely abort the Take off [10]. Both static and dynamic calculations, executed by the system, estimate the distance needed to accelerate to V1 and then decelerate to a full stop.

3.1. Static Calculations

The fundamental principle for statically calculating distance during the accelerating and the decelerating phase of a rejected Take off is Newton’s second law, in which the sum of the forces acting on a body is equal to the rate of change of body’s momentum:

$$\sum \overrightarrow{F}=m\overrightarrow{a}$$

(1)

Therefore, by modeling all the forces that act on the aircraft parallel to its longitudinal axis, with sufficient accuracy, its acceleration can be calculated at any given moment. Acceleration’s first integral is velocity, hence the bounds of the integral (time frame) of the accelerating and decelerating part can be found. By integrating once more, within the same bounds, we find the distance covered.

The algorithm that implements this double integral is based on the Riemann sum, in which the definite integral S of a function f is computed as a sum of area of parallelograms, with variable of integration, in this case, being time t:

$$S=\sum _{n=1}^{n}f\left(t_i\right)\Delta t_i$$

(2)

with $\Delta t_i=t_i-t_{i-1}$.

Figure 1 exhibits the algorithm that calculates the distance required to attain a certain speed, that being V1 for acceleration and 0 for deceleration. The Velocity and Distance variables are initialized to v_init (0 for acceleration and V1 for deceleration) and 0 respectively. For each repetition of the loop, we firstly calculate the sum of forces acting on the axis parallel to the runway $\sum \overrightarrow{F}$, and we divide with the mass to find acceleration $\overrightarrow{a}$. The following two commands essentially implement a right Riemann sum, and update the values for speed and distance. The loop exits when the final speed is reached.

The division (time step) used directly affects the accuracy of the result. By reducing the time step, the accuracy increases as well as the operations needed. For this particular implementation, the golden medium between the two, was found to be at 0.05 seconds, as it leads to instant calculations, without sacrificing the accuracy (deviation less than 1 metre).

As far as the right Riemann sum is concerned, we found that the time step is of such a scale that switching between the right, left and middle Riemann sum does not lead to a notable difference to the final result. Hence, the right is chosen, since it needs the less calculations and the algorithm will be embedded into a resource limited device.

3.1.1. Forces Acting on the Aircraft

The forces acting on the aircraft, that are parallel to the runway, are those that change the speed of the aircraft. These are the thrust of the engines $F_{Thrust}$ , the aerodynamic drag $F_{Drag}$, the friction between the tires and the runway $F_{Friction}$ and the component of the weight on the axis of the runway $m\ast g\ast sin\left(\theta \right)$ which depends on its slope $\theta$, as shown in Figure 2.

Each of the forces is being studied separately in the following paragraphs.

3.1.1.1. Engine Thrust

The Boeing 737-800 has two turbofan CFM56–7B26 type engines. Turbofan engines have expected performance characteristics that depend on altitude and airspeed [11]. The diagram for the aircraft’s engines is available by the simulator, as depicted in Figure 3. The produced thrust in the diagram refers to thrust production at 100% N1.

In order to make use of the information in the diagram, the thrust that is in pounds has to be converted into Newtons, and the altitude to the corresponding air density. Taking into account that the X-Plane simulator uses the International Standard Atmosphere (ISA), the air density on the surface of the sea is exactly 1.225kg/m³ with the temperature being at 15^oC and the pressure at 1013.25 hPa. The corresponding densities at the heights of 5000ft, 10000ft and 15000ft are 1.055 kg/m³, 0.904 kg/m³ and 0.771 kg/m³ respectively.

Only a certain range of data is extracted from the diagram of Figure 3 which contains the first four altitude lines of 0ft, 5000ft, 10000ft and 15000ft and velocities from 0 up to 0.3 Mach (198 Knots) with a step of 0.01 mach, with a total of 124 data points being collected. This certain range was selected because the calculations should be applicable to any airport of the world, with the highest being Daocheng Yadingat in China at an altitude of 14471ft, while the speed of 0.3 mach (198 Knots) is way greater than the V1 speed for this certain aircraft, in any possible configuration.

Based on these data and a polynomial fit, a function is generated, that exports the thrust at any density and airspeed, within range. The following polynomial is used to fit the data to.

$$\begin{array}{c}\hfill f(x,y)=A+Bx+Cy+D{x}^2+E{x}^{2}y+F{x}^2{y}^2+G{y}^2+Hx{y}^2+Ixy\end{array}$$

(3)

The selection of this quadratic function in two variables, with cross terms is done for three reasons. Each altitude line can be described by a second degree polynomial. Additionally, the distance between each altitude line decreases as the height increases, which also signifies a second degree relation between the thrust and the altitude. Finally, the cross terms (those who contain both x and y variables), are used as they are necessary for an accurate fit of this complex correlation.

Subsequently, the selected velocities and densities construct the meshgrid A, containing 124 rows (all the combinations of 31 velocities in 4 different densities) and 9 columns (as many as the terms of Equation 3). Each row represents a combination of x’s and y’s, i.e. velocities and densities.

Columns include the result of the variable related part of the terms of Equation 3. The first column contains all 1’s, as there it is a non-variable related term. The second column contains only the x (velocity) related to that row. The eighth column contains the product of the row’s velocity multiplied by the square of the row’s density and so forth.

Column B contains the samples of all 124 points of the diagram. The problem can be solved as a linear system $Ax=B$, in which x are the coefficients A,B,...I of Equation 3. The solution is calculated by solving through the least squares method. The Mean Absolute Error for the solution found is 7.005 Newtons.

The optical representation of the retrieved data and the fitted polynomial is depicted in Figure 4.

3.1.1.2. Friction between the Runway and the Tires

When the aircraft is on the runway, the only point that interacts with it are the tires, between which friction forces are being developed. They may be almost negligible during the acceleration, especially when comparing them with the thrust produced by the engines, but the aircraft’s braking system is the main stopping mean when decelerating. This makes its modeling vital.

The friction produced is calculated as the total force acting vertically on the tire multiplied by a friction coefficient $\mu$, as shown in Equation 4.

$$F_{Friction}=F_{Vertical}\ast \mu$$

(4)

For simplicity reasons and as no significant deviation emerges from it, the force was considered to be a point force. Hence, based on Figure 2, the friction force is expressed as:

$$F_{Friction}=v\ast (m\ast g\ast cos\left(\theta \right)-Lift)\ast \mu$$

(5)

in which Lift is the aerodynamic lift produced by the wings of the aircraft. The above equation makes the modeling of the Lift produced by the wings mandatory.

3.1.1.3. Lift and Drag

The two forces produced by the interaction of the body with the surrounding air are the Lift and Drag. Calculated by similar equations, these can be expressed as:

$$Lift={\displaystyle \frac{1}{2}}{C}_L\rho A{V}^2$$

(6)

$$Drag={\displaystyle \frac{1}{2}}{C}_D\rho A{V}^2$$

(7)

in which:

• Lift : Aerodynamic lift
• Drag : Aerodynamic drag
• C_L: Lift Coefficient
• C_D: Drag Coefficient
• $\rho$ : Air density
• A: The projected area of the body vertical to the force
• V: The airspeed of the aircraft

It is well established that the lift and drag coefficients are heavily depended on the Angle of Attack (AoA), which is the angle difference between the airflow and the lateral axis, that passes through the nose and the tail of the aircraft. Both lift and drag coefficients have typical curves, similar to these in Figure 5.

While on the runway, the AoA is close to zero with deviations of ±1.5 degree as found through various tests in the simulator. Based on the typical curves, a strongly positive correlation is expected between the lift coefficient and the AoA, while on the other hand no correlation between the drag correlation and the AoA.

To model the coefficients, Equations 6 and 7 are solved for them:

$$C_L={\displaystyle \frac{2Lift}{\rho A{V}^2}}$$

(8)

$$C_D={\displaystyle \frac{2Drag}{\rho A{V}^2}}$$

(9)

Necessary code is developed, that retrieves the Lift and Drag forces, the air density and the velocity of the aircraft of each frame of the simulation, and calculates the coefficients through reverse engineering.

The crew of the aircraft have a set of flap settings to choose from. This aircraft’s set includes the settings notated as 0,1,5,10,15,25,30,40. For each of these settings, flaps extend progressively more and achieve higher lift coefficients and therefore higher lift forces in the same speed. The counterbalance of these higher lift coefficients are also progressively higher drag coefficients for each of the settings. The five most common settings during Take off and those modeled for the system are 1,2,5,10 and 15. Similar methods can be used to model all the flap settings.

3.1.1.4. Weight Component

One of the characteristics of the runway is its slope. This is described as the height difference between the start and the end of the runway divided by its length, and is expressed as percentage (Equation 10).

$$Slope={\displaystyle \frac{\mathrm{Start}\mathrm{point}\mathrm{altitude}-\mathrm{End}\mathrm{point}\mathrm{altitude}}{Length}}\ast 100\%$$

(10)

In order to transform the slope percentage to degrees ($\theta$), the inverse tangent is used, as in Equation 11.

$$\theta =arctan\left({\displaystyle \frac{Slope}{100}}\right)$$

(11)

Given the angle $\theta$, the gravitational force can be analyzed in x and y axes, with one being perpendicular and the other being vertical to the aircraft’s axis of movement:

$$F_x=mgsin\left(\theta \right)\mathrm{and}{F}_y=mgcos\left(\theta \right)$$

(12)

Depending on the runway slope, $F_x$ either contributes or resists the achievement of a desired velocity, while $F_y$ is needed to calculate the total weight acting on the wheels.

3.1.2. Calculating Forces While Acceleration

There are two subsections about calculating the forces during a Take off. Since the system aims to calculate the distance that will be covered by the aircraft in a occasion of a Rejected Take Off (RTO) at speed V1, there are different calculations for the accelerating and the decelerating part of the procedure. The discrepancy between the two aims to the minimize the computations, and hence the time, needed to provide a result, without sacrificing its accuracy.

Engine Thrust

To calculate the thrust of the engines, a function is used to calculate the thrust based on density and airspeed. This output refers to thrust produced by one engine at 100% N1, and is hence multiplied by 2, that equals the number of engines.

Depending on the selected by the pilots N1 (primary thrust indicator), that may occasionally exceed 100%, the result has to be multiplied by the equivalent fraction. However that the relation between the N1 setting and the thrust produced is not linear, hence the result cannot be multiplied by the N1’s percentage. The X-Plane simulator provides the mapping between N1 and produced thrust, that can be used to calculate the actual thrust in relation to N1. For simplicity reasons and since this paper aims to prove the feasibility of such a system, it is considered that a constant N1 of 104% was used. This N1 setting corresponds to 101.5% of the maximum thrust produced.

Additionally, the engines have a spoolup time that is needed to reach the desired maximum thrust. This time is calculated to be 7.5 seconds for this type of engine. The modeling of this transition utilizes a rational function with the general form of:

$$F\left(t\right)={\displaystyle \frac{a{t}^2+bt+c}{d{t}^2+et+f}}$$

(13)

Then, the appropriate coefficients a,b,c,d,e,f are found, that describe the thrust for these 7.5 seconds as a fraction of the engine’s final thrust, in relation to the elapsed time (seconds), since the setting of the engine levers. The results of the actual and calculated spoolup thrust are presented in Figure 6 and Figure 7. The mean difference is at approximately 891 Newtons, which is negligible and not able to affect the final result.

Lift and Drag Coefficients, and Friction through Lift

The horizontal stabilizer is set, based on the weight and the center of gravity of the aircraft. This provides aerodynamic stability around the lateral axis, both on ground and in the air.

The above mentioned setting ensures that the AoA of the aircraft throughout the Take off will be around 0 degrees, with minor fluctuations of ±1 degree, as it is observed through testing. The lift coefficient faces big changes in this range. On the other hand, the friction force, that depends on the lift force, is multiplied by the rolling friction coefficient that is equal to or less than 0.025, on any runway condition, throughout acceleration. To this end, for simplification reasons and only for the accelerating part of the Take off, it is concluded that a mean lift coefficient can be used for each of the flap settings. On the other hand, the drag coefficient faces little to no changes within this range, as seen through the typical graphs, hence the mean value for each flap setting can also be used, without heavily affecting the final result as well.

Weight Component

The weight component is calculated once, through the $\theta$ angle (as described before), and is used as a steady force in the calculations.

3.1.3. Calculating Forces during Deceleration

Equally, if not more, important than the aircraft’s ability to attain speed, is the ability to reduce it. As a structure that weighs between 52000kg and 78000kg [13], it is quipped with a number of means that can reduce the kinetic energy, those being brakes, spoilers and thrust reversers.

The aircraft’s braking system works under the same principles as in cars, meaning that it transforms kinetic energy to thermal energy through friction. The overheating of the braking system leads to reduced performance and this is the reason that more than one stopping means are mandatory for modern medium and large size aircrafts. The second stopping mean is the wing’s spoiler. Placed on top of the wing, it differs from speedbrakes, as not only it increases the drag coefficient, it also disrupts lift production. Last are the thrust reversers, moving surfaces that divert the engine’s exhaust airflow, resulting in forces opposite to the aircraft’s moving direction.

The aircraft also has a number of options for the autobrakes, those being named RTO, 1,2,3,MAX. The last four selections concern the Landing of the aircraft and aim to achieve certain deceleration rates. The RTO selection though is used in the case of a Rejected Take Off and applies full pressure to the aircraft’s braking system, enabling the anti-skid system, until the full stop.

The calculation principle in the case of deceleration is common to that of the acceleration, meaning that the distance is calculated as a sum of parallelograms, that derive from acceleration’s double integral.

Engine Thrust

When decelerating there are practically three modes for the engines to work in. The one is the idle thrust, in which the engine produces the minimum thrust pushing the plane forward, the reversers in idle thrust, that produce the minimum reverse thrust and the reversers in the maximum setting. Apart from calculating the thrust produced in each setting, as a fraction of the maximum produced thrust, the transitions to each one of the abovementioned states is also modeled using the rational function used to model the spoolup transition (Equation 13). Indicatively, the modeled and actual transition from full throttle to full reverse throttle can be seen in Figure 8.

Lift and Drag Coefficient

Throughout acceleration, the friction force is orders of magnitude smaller than those of deceleration. Using mean lift and drag coefficients in this case, leads to fault-propagation and inevitably to unrealistic deviations in terms of distance.

A more elaborate modeling of the lift and drag coefficients is done by applying a least squares first degree polynomial fit to the collected data. In Figure 9, the data for Flaps setting 1 are presented, along with the resulting fitted line. The linear behaviour of the lift and drag coefficients in this range of Angles of Attack is already established.

The lift coefficient for each of the Flap settings, can now be expressed through this general form equation:

$$C_L\left(AoA\right)=a\ast AoA+b$$

(14)

In Figure 10 there are the fitted lines for all five of the modeled flap settings. It is clearly visible and expected that higher Flap settings, offer higher lift coefficients.

Each of the flap settings is also modeled with the spoilers deployed. Thus, any possible configuration can be predicted. The exact same process for both deployed and undeployed spoilers is used to model the drag coefficient in each one of the flap settings, as presented in Figure 11.

Angle of Attack

The modeling of AoA is now required in order to calculate the lift and drag coefficients. The factors affecting the AoA are divided in two categories, named primary and secondary. The primary factors include:

• Flap setting
• Runway condition
• Total aircraft weight
• Ground speed

While the secondary factors include:

• Ambient temperature
• Usage of speedbrakes
• Thrust mode

Regarding the primary set of factors, the first two are described by discrete values. The flap setting can be 1, 2, 5, 10 or 15 and the runway condition can be dry, wet or damp, as in the simulator. The last two though are described by continuous values. In order to model the AoA, the weight spectrum was sampled in seven values. Two of them are the minimum and maximum weight of the aircraft (51860kg and 79091kg). The range in between is sampled in 56662kg (5000 kg payload, 2 hours in normal cruise worth of fuel), 60463kg (7500kg payload and two and a half hours autonomy), 64495kg (10000kg payload and 3 hours of flight autonomy), 70641kg (15000kg payload and 3 hours of fuel autonomy) and 75826 (17000kg payload and 4 hours of fuel autonomy). These weights are selected to represent possible flight configurations for short or longer range routes. The interval between each two weights is about 5000kg, as it was noticed that there are no notable differences between smaller intervals.

Two runs are recorded for all the possible combinations between aircraft weight, flap setting and runway condition. A common V1 speed of 70m/s (136 knots) is selected among all runs, on which the RTO is be initiated, as all runs must contain similar data, so they can be further analysed. This speed is the minimum allowable speed among the configurations, that will not lift the aircraft off the ground. During the deceleration, the speed and AoA are being recorded. A total of 210 decelerations are executed in this experiment, for 105 different combinations of conditions.

Then, the processing procedure initially includes the fitting of the data for each couple of runs to a second degree function, using the least squares method. Figure 12 presents the data of AoA, as well as the fitted second degree function. The figure refers to a experiment configuration of 51860kg, Flap setting 1 and a dry runway.

The following step includes the scaling of the modeling on the collected data. As previously mentioned, the primary factors of weight and velocity are continuous, thus the same practice as in thrust calculation can be used. The data combination and results of Flap setting 15 and a Dry runway can be seen in Figure 13. More precisely, the blue lines represent the fitted second degree functions on the data of all 7 different weight settings, run with a Flap setting of 15 and on a dry runway. Then, these results are fitted to form a plane (seen in red dots), as in the thrust modeling. A total of 15 planes (for all possible combinations of 5 flap settings and 3 runway conditions) are constructed this way.

The secondary factors that affect the AoA are put in this category as they affect the whole plane equally on the vertical axis. Thus, a correction factor has to be added to the output of the corresponding plane. For example, negative temperatures affect the friction between the tires and the runway, as water and vapors turn into ice, which leads to reduced maximum braking capability and thus Angles of Attack closer to 0. Deployed spoilers are found to have the least affect on the AoA, but not at a negligible scale. The third and last secondary factor, that concerns the thrust mode (idle, reverse idle and reverse max), is found to have quite a noticeable impact on the AoA while decelerating, but still only on the vertical axis.

A total of 36 correcting factors are calculated for each of the combinations between temperature (positive or negative), the three possible conditions of the runway and the three possible thrust states. They are organised in two dimension arrays and integrated into the AoA calculating function.

Tire Friction While Braking

The RTO setting at the autobrakes system applies full pressure to the pads, activating the anti-skid braking system of the aircraft. This leads to fluctuations of the actual friction force between the tires and the runway. The maximum force that can be applied before the tire starts skidding, actually depends on two factors: the vertical force acting on it and the friction coefficient $\mu$. This coefficient describes the maximum ratio between the friction force, T, and the vertical force acting on the wheel, N, without loosing grip.

$$\mu ={\displaystyle \frac{T}{N}}$$

(15)

Friction coefficient $\mu$ depends on many factors, such as the type of the runway and its condition, the temperature, the tire condition, pressure and deformation. Since it is hard to determine the tire related factors in real world, the factors considered to affect friction coefficient, are the temperature (above or below $0^{o}C$, the freezing point of water), and the runway condition (dry, damp or wet). The only type of runway that commercial aircrafts use are pavement runways, thus this is the only type of runway examined.

Friction coefficient on pavement usually has values between 0.1 and 0.8 [14], depending on the condition of the runway and the ambient temperature. In order to take into consideration the anti-skid braking system of the aircraft, meaning the time intervals during which the pressure is released for the tire to roll and regain grip, the mean values of the coefficients were calculated.

More precisely, a plethora of tests are conducted, in the whole range of available temperatures of the simulator ($-{24}^{o}C-+{52}^{o}C$), in all available runway conditions (Dry, Damp, Wet). We concluded that temperature only affects the friction coefficient depending on whether it higher or lower than $0^{o}C$, which is the freezing point of water. Thus, it is expected to observe smaller friction coefficients in negative temperatures, as well as when the runway is wet.

The total combinations for the three runway conditions and the two temperature states (positive or negative) are six, as many as the output friction coefficients. They can be seen collectively in Figure 14. Hence, these are the friction coefficients that are used while braking in the RTO auto brakes setting, that take into account the anti skid braking system.

Weight component

The weight component is calculated once through $\theta$ as described before, and used as a steady force for the calculations

3.1.4. Additional Necessary Calculations/ Conversions

Air Density

Air density plays a significant role over the aerodynamic calculations as well as the engine thrust. To calculate the air density, only two variables are necessary, those being the ambient temperature and the barometric pressure on the runway, named QFE.

The input variables are in degrees Celsius, regarding the temperature and in hecto Pascals for the pressure. The way to convert them into air density is through the ideal gas law:

$${\rho }_{air}={\displaystyle \frac{P_{air}}{R_D\ast Temperature}}$$

(16)

In which R_D = 287.05 JKg^-1K^-1 is the gas constant for dry air.

Headwind Calculation and Boundary Layer

All Take off and Landing procedures are preferably and usually executed with a headwind instead of a tailwind for safety reasons, as it allows for a greater air speeds in comparison to ground speed.

To find the relation between the two speeds, the wind is analysed in horizontal and vertical components to the runway axis. The only data needed, are the heading of the wind and the heading of the runway. The forming angle between the two, noted as $\alpha$ in Figure 15, is used to find the headwind component. In particular:

$$\alpha =Headin{g}_{Wind}-Headin{g}_{Runway}$$

(17)

Then the wind speed can be multiplied by the cosine of $\alpha$, to find the headwind component:

$$Headwin{d}_{speed}=Win{d}_{speed}\ast cos\left(\alpha \right)$$

(18)

The input and the output share the same measurement unit. Such data are available to the pilots, by the control tower.

The boundary layer that is formed close to the ground, which reduces the speed of the air, has to be taken into account as well. Various tests with different wind speeds and directions, unveil that the wind’s speed can be calculated if multiplied by coefficient $CBoundary=0.6617$. Hence, the headwind on the runway is expressed as:

$$Headwin{d}_{onrunway}=Headwin{d}_{speed}\ast C_{Boundary}$$

(19)

Velocity Conversion

Avionics offer a variety of information including the Indicated Air Speed (IAS), the True Air Speed (TAS) and the Ground Speed (GS). The IAS is displayed in the primary flight display of the aircraft and the speed that V1, VR and V2 refer to, as well as V_App and V_Ref for the Landing procedure. The conversion from IAS to GS i required to find the covered distance.

The different states followed to convert IAS to GS can be seen in Figure 16. The first step converts IAS to Calibrated Air Speed (CAS). Systems onboard the aircraft add an error to the actual speed of the aircraft. A collection of values for both IAS and CAS indicates that the relation between the two is expresses as in Equation 20.

$$CAS=0.989\ast IAS+5.7\left(Knots\right)$$

(20)

Calibrated Air Speed refers to speed in ISA (International Standard Atmoshpere) conditions, which means that the ambient temperature and the barometric pressure are not taken into consideration. For given conditions of temperature and pressure, the air density can be found, and the conversion to True Air Speed is done through Equation 21.

$$TAS=CAS\ast \sqrt{{\displaystyle \frac{{\rho }_0}{{\rho }_{air}}}}$$

(21)

in which ${\rho }_0$ is the ISA conditions air density of 1.225 kg/m³, and ${\rho }_{air}$ is the actual air density.

The last step is to subtract the headwind from the TAS as in Equation 22, in which all speeds refer to knots.

$$GS=TAS-Headwind$$

(22)

3.1.5. Overall Calculation of ASD Required

Since all the forces are modeled, it is possible to calculate the Accelerate Stop Distance Required (ASDR). The aggregate of inputs needed, includes:

• The externally calculated V1 speed, in knots
• Wind direction, in degrees
• Wind speed, in knots
• Ambient temperature, in degrees Celsius
• Ambient Atmospheric Pressure on Elevation, QFE, in hPa
• Flap setting
• Total aircraft mass in kg
• Runway slope (as percentage)
• Runway condition

There is also the capability of calculating the distances when either using or not the spoilers, as well as when selecting different modes in engines during an RTO. When aborting a Take off though, the pilots use all the braking means available, which means that they make use of the spoilers and the engines are set to the maximum reverse thrust. This is why these options do not necessarily have to be given as inputs.

Initially, the steady factors are being calculated, those being: the air density, the lift and drag coefficients during accelerating, the weight component, the headwind and the friction coefficient while decelerating. Then the algorithm executes the two while loop, to calculate the acceleration and deceleration distances, with a top speed of V1, respectively. In each iteration of the two loops, the rest of the calculations (such as aerodynamic forces) are done based on the current airspeed or groundspeed of the aircraft.

As mentioned, the V1 speed is given by the user. A minor alteration in the calculations and the additional input of the runway’s length, would allow for the calculation of the V1 speed by the system. In reality though, the majority of the runways, allow the attainment of $V_R$ speed (speed in which the nose of the aircraft is lifted and the actual Take off takes place), long before the aircraft is in danger of overrunning the aircraft. Additionally, the relations between the critical speeds of the aircraft are explicitly stated in EASA’s CS-25 file [9], and declare that $V_R$ speed may not be less than V1 speed.

In this condition, the strict definition of V1 speed is abolished and the system is unable to calculate the correct speed. Furthermore, discrepancies between the system’s and other systems’ outputs (those that are exclusively made to calculate the critical speeds), would create confusion to the pilots, which should be avoided.

In addition to the previous calculations, a distance equivalent of two seconds at speed V1 is added on top of the previous result. This is considered to be the reaction time in case that an engine failure happens at speed V1. It is an indivisible part of the Accelerate Stop Distance calculations, as CS-25 [9] states.

3.2. Dynamic Distance Calculations

Static calculations depend on many different factors and assume a smooth run in each process. However, one of the main and essential functions of the system is the prediction of distances in real time, in the event that the Take off or Landing of the aircraft does not run as predicted.

The dynamic calculations performed during Take off are not based on any of the inputs mentioned in the static calculations. The only input they rely on is the acceleration (or deceleration) of the aircraft. In this context, it is necessary to analyze the data as well as the theory, to find a computationally efficient method, capable of coping in real time, taking measurements of the acceleration at regular time intervals and performing the appropriate calculations.

3.2.1. Acceleration

Figure 17 shows a typical acceleration graph of a Boeing 737-100, which was modified and used by NASA (National Aeronautics and Space Administration). It is clearly visible that after spoolup (that usually lasts less than ten seconds, depending on the engine), there is a linear decrease in acceleration. This is because as the speed increases, the produced thrust of the engines decreases and the aerodynamic resistance increases. At the same time, the friction between the wheels with the ground is reduced, but is the force with the least impact.

After collecting data for the Boeing 737-800 on three random configuration, a corresponding graph was constructed. The expected linear relation after the spoolup, is visible in all three configurations of Figure 18. More precisely, after the 7.5 second period that the engines need to reach full power, the acceleration rate decreases linearly in time.

3.2.1.1. Calculation of the remaining distance

To calculate the distance needed to achieve a certain velocity, we find the $\alpha$ and $\beta$ coefficients that apply to the linear part of acceleration, when expressing it as one:

$$Acceleration\left(t\right)=\alpha \ast t+\beta$$

(23)

Then, by integrating Equation 23, the velocity equation derives:

$$Velocity\left(t\right)=\int Acceleration\left(t\right)={\displaystyle \frac{\alpha }{2}}{t}^2+\beta t+\gamma$$

(24)

$\alpha$, $\beta$ are the same coefficients as in Equation 23 and $\gamma$ is the already attained speed. Practically, it is the speed that the aircraft has achieved during the spoolup, right up to the start of the linear part. Since the speed that needs to be achieved is already known (speed V1), the time needed can be found by the solutions to this quadratic equation:

$${\displaystyle \frac{\alpha }{2}}{t}^2+\beta t+\gamma -V1=0$$

(25)

Taking into consideration Equation 25 and its coefficients, its discriminant is:

$$\Delta ={\beta }^2-4{\displaystyle \frac{a}{2}}(\gamma -V1)$$

(26)

While the solutions are:

$$t_{1,2}={\displaystyle \frac{-\beta \pm \sqrt{\Delta }}{2\frac{a}{2}}}$$

(27)

The $\alpha$ coefficient always has a negative value, as the acceleration rate tends to decrease during the procedure, while the $\beta$ coefficient always has a positive value, representing the maximum achieved acceleration rate, immediately after the completion of the spoolup. Hence, projecting the accelerating function in time, will eventually attain negative values. Then, the velocity starts to decrease due to the negative acceleration. This means that, the solutions to the quadratic function, if any, will refer to the moments that V1 has been achieved, with positive acceleration values and negative acceleration values. The only real solution though, is always the first of the two.

It is easily understood that there will always be solutions to satisfy the equation. Non existence of solutions would mean that the aircraft is not able to attain a speed equal to V1. In addition the linearity of the acceleration function ceases to be in effect, when it gets close to 0, as the aircraft gets closer to its terminal velocity.

Based on the above, the solution will always lie in the first of the two solutions. More precisely, based on Equation 27, the solution is:

$$t_{V1attainment}={\displaystyle \frac{-\beta +\sqrt{\Delta }}{\alpha }}$$

(28)

Having found t_{V1 attainment}, it is also possible to find the distance that will have been covered during that time as well. By integrating the velocity function, we get the distance function in time:

$$Distance\left(t\right)={\displaystyle \frac{\alpha }{6}}{t}^3+{\displaystyle \frac{\beta }{2}}{t}^2+\gamma t+\delta$$

(29)

in which $\alpha$ and $\beta$ the coefficients from Equation 23, $\gamma$ the speed immediately after the completion of the spoolup and $\delta$ the distance already covered at t=0. By replacing t with t_{V1 attainment}, the distance is calculated.

3.2.1.2. Approximation of Acceleration’s $\alpha$ and $\beta$ Factors

An equally important feature is the approximation of the $\alpha$ and $\beta$ coefficients for the accelration equation. Firstly, the non linear data of the spoolup’s 8 seconds (Figure 18) are not taken into consideration.

The next 3.5 seconds worth of data (time and acceleration) are used in a linear regression in order to determine the gradient of the function. The intercept is determined by the first collected value. These two values are used as seeds for the next algorithm.

Theoretically, all the samples that we get later than the 11th second can be saved in constantly increasing arrays and then processed by a linear regression algorithm. However, this solution is not efficient for a real time system, due to the exponential complexity of the algorithm. Instead, the Recursive Least Squares filter (RLS) is used, that uses the weighted least squares method to approximate coefficients $\alpha$ and $\beta$, without having to execute increasingly more calculations for every incoming measurement of acceleration. RLS is used to approximate optimal coefficients of linear equations, structured like:

$$y=Ax$$

(30)

The algorithm also has an exponential complexity but the simplicity of the system, consisting of only two coefficients, keeps the total calculations at a consistent low number for each incoming data package. The linear equation that describes the acceleration in time can be also expressed as:

$$acceleration\left(t\right)=\left[\begin{array}{cc}a& b\end{array}\right]\ast \left[\begin{array}{c}t\\ 1\end{array}\right]$$

(31)

The RLS algorithm can be tuned through the forgetting factor $\lambda$, which sets the significance of older samples. In this configuration, the forgetting factor was set to 1, as all the samples need to be of the same significance in order to output an accurate result and in order to converge faster.

3.2.2. Deceleration

The decelerating part of a Rejected Take Off also has to be calculated dynamically. The RTO setting of the autobrakes though leads fluctuating braking force and thus fluctuating accelerations. The RLS filter can also be used to calculate the distance, but not with an identical structure. The alterations in measurements set the convergence impossible.

3.2.2.1. Calculation of the Remaining Distance

Due to the above reasons, the speed and the quadratic relationship approach that describes it are used, instead of the acceleration. The quadratic equation that describes velocity in time is:

$$V\left(t\right)=\alpha \ast t^2+\beta \ast t+\gamma$$

(32)

in which $\alpha ,\beta ,\gamma$ the coefficients of the equation and t the time. Similarly as before, the solutions to the quadratic equation can determine the moment in time at which the aircraft has zero velocity. The two solutions to the equation being zero, the final desirable speed and in correspondence to the coefficients of the previous equation, are:

$$t_{1,2}={\displaystyle \frac{-\beta \pm \sqrt{{\beta }^2-4\alpha \gamma }}{2\alpha }}$$

(33)

The deceleration rate during this procedure is not steady. The more the velocity drops, the more lift force is reduced. This leads to greater vertical forces on the wheels, and thus capability to apply more baking force before the tire starts skidding. Based on this observation, it can be determined that $\alpha$ coefficient of the velocity function will always be negative and $\beta$ will certainly be negative.

Based on the above ($\alpha <0,\beta <0,\gamma >0$ being equal to to the speed of the aircraft), the equation has an inverted bell shape, whose maximum value lies above the horizontal axis. The graph meets the axis in two points one in negative t values and one in positive. The solutions that is of interest to us is the latter one.

Since, the signs of $\alpha$ and $\beta$ coefficients are both negative, the solution is found when subtracting the square root of the derivative:

$$t_{tillstop}={\displaystyle \frac{-\beta -\sqrt{{\beta }^2-4\alpha \gamma }}{2\alpha }}$$

(34)

Then, the t_{till stop} can be applied to the integral of the velocity function and the distance can be found. Thus, the distance is found through:

$$Distanc{e}_{tillstop}={\displaystyle \frac{a}{3}}{t}_{tillsto{p}^3}+{\displaystyle \frac{b}{2}}{t}_{tillsto{p}^2}+c{t}_{tillstop}+\delta$$

(35)

The constant term $\delta$ that results from the integral, is equal to the distance that has been covered before the start of the dynamic calculations and after the initiation of the Rejected Take Off.

3.2.2.2. Approximation of Velocity’s Coefficients

A rejected Take off is considered to start when the throttle levers are pulled back. The duration of the pull back is approximately 0.8 seconds. The next four seconds, many changes happen to the aircraft. The deployment of the spoilers and the spooldown of the engines affect the deceleration rate in non linear ways, thus these four seconds of data are being discarded. The next two seconds worth of velocity and time data are gathered and used in a least squares polynomial fit, in order to calculate the seeds of the last algorithm used. These seeds consist of three coefficients, $\alpha ,\beta$ and $\gamma$, in correspondence to Equation 32.

Then, once again the Recursive Least Squares filter is used, with a slightly different configuration. The filter now tries to converge to three coefficients, instead of two, while the problem can be described as a linear equation like:

$$acceleration\left(t\right)=\left[\begin{array}{ccc}\alpha & \beta & \gamma \end{array}\right]\ast \left[\begin{array}{c}{t}^2\\ t\\ 1\end{array}\right]$$

(36)

Once again, the forgetting factor $\lambda$ is set to 1, as all the samples are of equal importance.

4. Landing Procedure

Landing is the most critical part of a flight and it includes heavy work load by the pilots as technical and aerodynamic configurations have to be taken into account, as well as the ambient conditions at the Landing site. The accuracy and efficacy of the crew’s actions directly affect the passengers’ convenience and safety.

The most common type of accident during Landing is the runway overrun [1], with causing factors including high approach and touch down speeds, insufficient use of the reversers etc. [16].

4.1. Velocities and Parts of Landing

Before the Landing procedure starts, two critical velocities are calculated, those being V_App and V_Ref, based on a number of factors similar to those of Take off. V_App is the speed at which the aircraft approaches the airport, while for the final part, this speed is reduced down to V_Ref. This is practically the stall speed multiplied by a safety factor of 1.23, as CS-25 defines. This ensures a safe Landing along with a minimum Landing distance requirement.

The Landing procedure can be divided into four stages: Approach, Flare Free roll and Braking, as Figure 19 indicates.

The approach phase starts at the threshold of the runway at which supposedly the plane is at a height of 50ft and a speed equivalent to V_Ref. The glide path followed may vary depending on limitations in the landscape around the runway, but it is 3 degrees (or less) for the vast majority of cases. Then, the flare phase starts at 20-30ft above the runway, depending on the aircraft itself. This maneuver, aims to reduce the vertical speed from ≈700ft/min, down to 60-180ft/min, so that the aircraft can touch down smoothly. Elongated or early initiation of the maneuver can set the procedure in danger of a runway overrun, as the aircraft may be found hovering above the runway, wasting meaningful metres.

Immediately after touching the runway, the free roll stage begins. It lasts about 2 to 4 seconds, and aims to lower the nose of the aircraft, so that all wheels are in contact with the ground, before applying any of the braking means. Last but not least, the braking stage utilizes at least one of the decelerating means, till it comes to a full stop, or it attains a taxiing speed of 20-30knots. Since the system must be able to foresee an overrun, it calculates the distance needed to come to a full stop.

4.2. Static Calculations

The static calculations are performed after the user has given all the inputs needed to predict the total Landing Distance Required (LDR). Each part is estimated separately, in ways described in the following paragraphs.

All tests and modeling for the Landing phase, are executed using the autopilot function of the aircraft, on a Instrument Landing System (ILS) CAT-III equipped runway. In each part, there is a tendency for a slight overestimation (5-15%) by using mean to upper factors where applicable, as a safety measure, in cases that the system fails to estimate LDR accurately.

Approach

The usual approach glide path for commercial aircrafts stands at 3 degrees. Additionally, through simulations, it is concluded that the optimal height to initiate a flare for the Boeing 737-800 is at 30ft. This means that the height difference between the two stages is at 20ft.

In Figure 20, the velocity component analysis for this specific glide path is shown, separating V_ref into its vertical and horizontal components.

The time it takes for the aircraft from 50ft to 30ft is:

$$t_{50ft-30ft}={\displaystyle \frac{50ft-30ft}{V_{Ref}\ast sin\left(3^o\right)}}$$

(37)

where t refers to seconds and V_Ref is converted to feet per second. The distance covered during that time is:

$$D_{50ft-30ft}=t_{50ft-30ft}\ast V_{Ref}\ast cos\left(3^o\right)$$

(38)

where D_50ft-30ft refers to feet and V_Ref is converted to feet per second.

The overestimation alterations for this stage included the reduction of the theoretical glide path to 2.86^o and the disposal of cos(3${ }^0$) for calculating the horizontal speed. Thus:

$$t_{50ft-30ft}={\displaystyle \frac{50ft-30ft}{V_{Ref}\ast sin(2.{86}^o)}}$$

(39)

$$D_{50ft-30ft}=t_{50ft-30ft}\ast V_{Ref}$$

(40)

Flare

Modeling the flare maneuver through a rational function is considered to be meaningless, since two flares may differ vastly. Instead, numerous tests are executed in order to find the mean values that described the outcome of the process. It is concluded that a flare lasts 4 to 8 seconds, with a tendency towards the lower limit. For the static calculations it is considered to last 5.5 seconds, slightly more than the average time. The average velocity maintained during flare is found to be 2.1% less that V_ref. Hence, the distance covered during Flare is:

$$D_{Flare}=5.5\ast 0.979\ast V_{Ref}$$

(41)

where D_Flare and V_Ref are in feet per second.

Free Roll

Through similar processes, it is concluded that there is an average loss of 7.5% off the V_ref speed at the moment that the aircraft touches the ground and is considered to remain steady for the duration of this phase. Free roll lasts for 2 to 4 seconds, with the mean to upper value of 3.5 seconds being used for the static calculation. The free roll distance is calculation is:

$$D_{FreeRoll}=3.5\ast 0.925\ast V_{Ref}$$

(42)

where D_{Free Roll} and V_Ref are in feet per second.

Braking

For the fourth and last part of Landing, the pilots have the option to either use the braking system manually or by using the auto brakes in one of the four available settings, 1, 2, 3, MAX. The static calculations have to be done based on an expected behaviour, so there is an input for the Landing procedure with these settings available.

Each one of them achieves a certain deceleration rate during the braking stage, ranging from 1.2 m/s${ }^2$ to 4.1 m/s${ }^2$[17]. Additionally, ranging temperatures have a minor impact to this rate, while thrust reversers add 0.05m/s${ }^2$ to 0.22m/s${ }^2$ to this rate.

Based on the expected deceleration rate, the braking distance is found. Velocity while braking can be expressed though this time function:

$$V\left(t\right)=V_{SB}+acceleration\ast t$$

(43)

in which V(t) is the velocity in time t, V_SB is the velocity at t=0, at which braking started and acceleration isthe expected acceleration rate.

Setting the final desired speed to 0 and solving for t, we get the time needed till the aircraft stops:

$$T_{tillstop}={\displaystyle \frac{V_{SB}}{acceleration}}$$

(44)

By integrating the velocity function and replacing t with T_{till stop}, we get the distance covered while braking:

$$D_{covered}={\displaystyle \frac{acceleration\ast T_{tillsto{p}^2}}{2}}+V_{SB}\ast T_{tillstop}$$

(45)

4.3. Dynamic Calculations

As in the Take off procedure, dynamic algorithms are also developed for the Landing procedure.

Pre-Approach

As EASA’s Certification Specifications describe, a ROAAS system must be able to inform the crew for a longitudinal excursion in flight (from a given height) and on ground. For this reason, an extra stage named Pre-Approach is added to the dynamic calculations, starting at 330 feet above the runway level and ending at 50ft.

The system in this stage saves the last 100 values of vertical and horizontal speed and calculates the moving average of both of them. The sensitivity of the system can be tuned by changing the number of the last values kept. Based on the data output rate of about 30 samples per second that the simulator produces, this setup leads to a moderate sensitivity. Embedding the system to a real aircraft, would require to tune this number in regard to the data output rate that the aircraft’s avionics.

Based on the average vertical speed and the current altitude, the time needed to reach 50ft is calculated. This time is multiplied by the average horizontal speed to find the distance that will be covered till the aircraft reaches 50ft of altitude above the runway.

Approach

Similarly, the moving average of both horizontal and vertical speeds is calculated, but for this part, the algorithm is more sensitive keeping the 80 last samples. The new threshold is now 30ft above the ground, that the flare stage starts. Other than that, the calculations for the distance prediction are identical.

Flare

As previously mentioned, the flare maneuver, cannot be modeled or predicted, thus the same configuration is used with an even higher sensitivity, keeping only 50 of the last values. The threshold for the last airborne part of the Landing is 0ft. The calculations are identical to those of the previous two parts.

Free Roll

Immediately after touching down, the free roll stage starts. For the static and dynamic calculations, a duration of 3.5 seconds is considered.

Initially, the start time is saved to a variable, while the moving average of the last 40 values of horizontal speed in feet per second is calculated.

For every new sample, the past time from the start of the free roll is calculated:

$$t_{freeroll}=t-t_{freerollstart}$$

(46)

Then the time left to the 3.5s of the considered duration is computed:

$$t_{freerollleft}=3.5s-t_{freeroll}$$

(47)

Using the moving average of the horizontal speed, the distance left till the end of the free roll is determined:

$$D_{freerollleft}=V_{meanhorizontal}\ast t_{freerollleft}$$

(48)

In case that the braking distance is not initiated, before the end of the 3.5 seconds, then the time left is considered to be 0.5 seconds and the distance left is given by:

$$D_{freerollleft}=V_{meanhorizontal}\ast 0.5s$$

(49)

The reaction time of 0.5 seconds is selected because it sets an additional safety distance to the calculations and thereby gives time for the pilots to initiate the braking. At the same time it does not trigger the system unintentionally with an excessive overestimation.

Braking

The last part of Landing, braking, is triggered when pressure is applied to the brakes either by the auto brake system or by the pilots. The values of initial time and velocity are saved for further use.

Again the Recursive Least Squares algorithm that is utilized as in the deceleration of the aircraft at during a Rejected Take Off. The algorithm aims to minimize the sum of the squares of the difference between the approximated and the actual coefficients that describe the velocity as a second degree function:

$$V\left(t\right)=\alpha \ast t^2+\beta \ast t+\gamma$$

(50)

The difference is that the deceleration in this case is overall steady with minor fluctuations. The seed values in this case is 0.005 for coefficient $\alpha$, the anticipated deceleration rate for $\beta$ and the initial velocity at the start of deceleration for $\gamma$.

The reason for using an $\alpha$ coefficient different than 0, as it would be expected for a steady deceleration, lies at the way that the remaining distance is calculated. As a second degree function, the time at which the speed is eliminated can be found through the discriminant:

$$\Delta ={\beta }^2-4\alpha \gamma$$

(51)

And the actual time:

$$t_{tillstop}={\displaystyle \frac{-\beta +\sqrt{\Delta }}{2\alpha }}$$

(52)

By integrating the velocity function (Equation 50) and using the time found through Equation 52, the distance needed till the aircraft comes to a full stop can be found:

$$D_{tillstop}={\displaystyle \frac{\alpha \ast t_{tillsto{p}^3}}{3}}+{\displaystyle \frac{\beta \ast t_{tillsto{p}^2}}{2}}+\gamma \ast t_{tillstop}$$

(53)

The constant of integration that arises from the integral is practically zero, because the dynamic algorithm of this state starts after the initiation of braking.

Dynamic Calculation of LDR While LANDING

The principle behind dynamically calculating LDR, is that the remaining distance for the current stage is calculated in real time and the static estimations for the rest of the process are summed to output a single result.

In Table 1, the current stage is presented in each line of the table, while the type of result (statically calculated or dynamically calculated) used in the sum, is presented in the columns of the table.

The dynamic result has meaning only when compared to the position of the aircraft relative to the runway. For this reason, a data structure is made, including the coordinates of the thresholds of all runways, based on the simulator’s data. The input given by the user includes the ICAO code of the airport and the runway marking. The returned data are the elevation of the runway, its length and the threshold’s geographical location. Equivalent aeronautical data, such as the FAA’s National Airspace System Resource (NASR) System [18] can be utilized for a real world application.

Having the threshold’s and the aircraft’s geographical data, through the GPS system, their relative distance is calculated through the Haversine algorithm.

Having the absolute distance of the aircraft from the runway, let it be $D_n$ and the distance from the previous calculation, let it be $D_{n-1}$, we can find whether the aircraft has passed the threshold. If the sign of the subtraction $D_{n-1}-D_n$ is positive, this means that the newly calculated distance is smaller than the previous, thus the aircraft is still approaching the threshold. In the occasion that the sign is negative, then the aircraft is above the runway.

Let LDRD be the dynamically calculated LDR. If the aircraft has not passed the threshold, then:

$$LDR=LDRD-D_n$$

(54)

If the aircraft is past the threshold, then:

$$LDR=LDRD+D_n$$

(55)

Finally, since the length of the runway is already known, the comparison between LDR and the length can determine whether the aircraft is in danger of overrunning the runway.

5. Performance Evaluation

The performance evaluation of safety systems is a crucial part of a safety system’s development. This is why numerous tests were designed and executed, on the majority of the input parameters that a user can input on the system.

The only parameter that was not taken into consideration, even though the system fully supports its input, is the runway gradient. This was done for two reasons.

Firstly, the simulator has the option to either render the runway as a completely flat surface, or to add the gradient based on the surface characteristics. The latter though, creates unrealistic runway gradients, with steep bumps, that often times cause the aircraft to bounce on them. This creates misleading scenarios and the data collected from such runs could not be of any use. This is why the former option was selected.

The second reason is that, even if the slope was to be rendered as realistic as possible, the gradient of a runway is calculated as the height difference between the two thresholds. Any height differences in between are practically not taken into consideration. This is an already discovered problem with such systems [19], while the best promising method would be to collect analytical data for the elevation across the runways.

The modeFRONTIER automation and design software [20] is used to design the experiments. More precisely, the software has a Design of Experiments (DoE) mode, in which all the independent parameters, their range and the total desired number of experiments are entered. A set of experiments is exported based on a quasi random Sobol sequence, exploring the widest possible range of combinations.

A total of 200 experiments are designed for each procedure (Take off and Landing) and the corresponding data are collected and used to reproduce the results that the system would output, in both static and dynamic modes.

5.1. Take off performance evaluation

All the experimental runs, consisted of Rejected Take Offs exactly at speed V1. The Variables and the corresponding value range (or available options) can be seen in Table 2.

Velocity V1 is practically a dependent variable, thus a simplified equation is constructed to calculate the appropriate velocity for each of the runs. This equation is:

$$V1=⌊{\displaystyle \frac{Weight}{1000}}⌋+60Knots$$

(56)

This simplified equation outputs a range of V1’s from 112 knots to 138 knots depending solely on the weight of the aircraft. This aims to approximate the real V1 speed without getting close to the speed that the plane would get airborne.

Static Take off Calculation Performance

The static calculations, that the comparison is done with, do not include the 2-second at V1 distance, that EASA’s regulations set as a reaction time and safety parameter. This aims to show the exact difference between the static and the dynamic calculations, without taking into consideration the human factor.

Three comparisons are done to determine the accuracy of the system. The first is about the acceleration distance, seen in Figure 21. For each of the runs, the percentile deviation is calculated, while the results are presented in a histogram, along with the standard deviation and the mean value of these results. It is found that the average difference is -0.18% with a standard deviation of 0.55%. This means that in 95.5% of occasions, the difference is between [-1.28%, +0.92%].

Similarly, regarding the deceleration distance, seen in Figure 22, the mean difference is +0.07%, with a standard deviation of 1.15%, due to the higher complexity of the calculations executed.

Cumulatively, the ASDR distance is slightly underestimated by 0.02%, with a standard deviation of 0.73%, as Figure 23 indicates. This means that in 99.7% of times, the static calculations differ by a percentage within [-2.21%, 2.17%].

Dynamic Take off calculation performance

The performance of the dynamic calculations is evaluated in two ways, namely the time needed by the algorithm to converge and the deviation of the last prediction.

As seen in Figure 24, the first predictions of the algorithm, at 11 seconds, that derives from the seeds, underestimates the distance needed to achieve V1. Half a second later the mean prediction is at +2% and 12 seconds after the initiation, the mean prediction is at +1%. It comes though with a high standard deviation (light blue background that equals to $\pm 1\sigma$), that sets the result non utilizable.

Instead, the 15 second time mark has been set as the time after which, the RLS’s mean result is close enough to a 0% deviation and the standard deviation is around 1% with a tendancy to decrease even more in time. Thus, the algorithm has converged and the result can be safely compared with the statically calculated one.

Throughout the 200 runs, the last prediction of the dynamic algorithm has a mean difference of +0.29% when compare to the actual result, with a standard deviation of 0.26%. The histogram can be seen in Figure 25

A similar evaluation is also done for the decelerating part as well. In Figure 26, the mean estimation throughout time can be seen. At about 6.8 seconds after the initiation of a Rejected Take Off, the first dynamic result is available. Initially there is an underestimation of about 10%, but this changes rapidly. The mean overestimation peaks at about 2.5%, 8.5 seconds after the initiation, and then continuously converges to a 0% deviation.

Since there is no heavy overestimation that can falsely trigger the system with a warning, there is not a time mark equivalent to the previous procedure’s. All the dynamic results produced by RLS are treated equally and they are being compared to the statically calculated ASDR, that also includes a distance equal to two seconds at speed V1.

The algorithm’s last prediction can be seen in Figure 27. The mean deviation stands at -0.02%, while the standard deviation is 1.28%. This proves that till the end of the deceleration the algorithm has converged realistically close to the actual distance.

5.2. Landing Performance Evaluation

Another 200 experiments are designed for the Landing process through the same DOE generating software.

In Table 3, the corresponding variables and their ranges (or available options) can be seen. The maximum weight is adjusted to the maximum Landing weight suggested by the manufacturer [13]. The settings for flaps used during Landing are 30 and 40 and the Autobrakes can be set at setting 1,2,3 or MAX, and so it is in the experiments.

The output file is adjusted in certain ways. Specifically, all experiments that are executed with a tailwind greater than 10 knots, are manipulated in their wind direction or wind speed. Landings with tailwinds are also avoided in real world operations, for they enhance the danger of a runway excursion.

The Vref speed is calculated through a simplified algorithm, presented in Equation 57. It outputs a range of 129 to 143 knots in Flap setting 40 and 136 to 150 knot in Flap setting 30, based solely on the aircraft’s weight. This provides a rough approximation of the Vref speeds used in real world procedures. Nevertheless, the output speeds are greater then the V_s0 (stall) speed, in all cases.

$$Vref=\left\{\begin{array}{c}108+⌊{\displaystyle \frac{Weight-40000}{5}}⌋\ast 7,flaps=40\\ 115+⌊{\displaystyle \frac{Weight-40000}{5}}⌋\ast 7,flaps=30\end{array}\right.$$

(57)

5.2.0.15. Static Landing Calculation Performance

The results for the comparison between the static calculations and the actual distances per stage of Landing and for the total Landing Distance Required can be seen in Figure 28. All stages and the overall distance, are slightly overestimated, except for the Flare and Free Roll distance. Their mean deviation is around 25% - 30%, while they also have much larger standard deviations.

These data may at first sight be an issue for the performance of the system. However, as previously mentioned, flare is a hard to model and almost impossible to predict maneuver. An analysis between this stage and the free roll, that also has an even greater mean value, is conducted to reveal more insights about the two stages. Executing a linear regression between the data and visualizing it, as in Figure 29, it is obvious that there is a negative correlation between the deviations of the two stages. This means that the bigger the deviation is in the stage of Flare, the smaller it is in the next stage, that of the Free Roll, and vice versa.

Finally, a more detailed look in the results of the statically calculated LDR distance, reveals an average overestimation of 8.98% and a standard deviation of 5.05%. These results can be seen in Figure 30.

As depicted in the diagram, there are instances where the overall LDR distance is underestimated (Deviation < 0). However, upon individual examination of these instances, it becomes apparent that the disparity arises from the Flare stage, which is underestimated by an average of 27% (Table 4).

Similarly, corresponding results arise from the examination of cases where the estimation exceeded the actual result by at least 20%. The majority of the overestimation stems from the horizontal component (Table 5, study of three of the cases).

5.2.0.16. Dynamic Landing Calculation Performance

The performance of most of the Landing procedure’s dynamic stages cannot be represented collectively, due to the lack of a common basis for the data to be presented at. A common basis that could be used is time, but the varying and usually short duration of the stages sets such a comparison devoid of revealing useful and insightful information.

The part that could be analysed though, due to its greater duration, is that of braking, that utilizes the RLS algorithm. From the moment braking is initiated, the RLS is able to produce results, due to the almost complete absence of slope in the deceleration rate. Therefore, after processing the data from the 200 experiments, the average values and variances, during the first 8 seconds of braking, are obtained. As shown in Figure 31, initially there is a negative deviation in the result of approximately -7.5%. This changes at 4 seconds, reaching a maximum later, at 5.2%. From there on, the algorithm converges to 0%, with the prediction at 8 seconds having an average value of 0.13% and a standard deviation of 2.4%.

6. Conclusion

The aim of the system is to predict runway overruns in real time for both Take off and Landing procedures. Two types of algorithms were developed, including static and dynamic calculations of the remaining distance needed by the aircraft to come to a stop. The algorithms are developed frugally in order to minimize complexity, without compromising in terms of accuracy. Additionally, the vast majority of affecting factors are taken into consideration in the widest range possible, so that the system is able to perform in any probable scenario, with the same reliability.

All the system’s inputs are carefully selected, such that they are either available to the pilots, or accessible through the aircraft’s avionics system. This will lead to a less complex embedding process, along with the reliability of the already developed systems and processes.

The results of the evaluation process provided sufficient evidence that the algorithms can be used in a real time overrun avoiding system. Both procedures, Take off and Landing, can be predicted with adequate accuracy statically, while the dynamic algorithms are able to predict the distance needed based solely on the development of the process. Their comparison can provide a reliable result, evaluating the danger of an overrun, in the case that the Take off or Landing procedure is continued as is.

The algorithms have also been embedded to a RaspberryPi 3b, with a 7 inch touch display and speakers for interaction purposes, to form a prototype system. The Raspbian operating system is used, while the application is developed through the Qt application development framework. All the algorithms are developed in C++, due to adaptability and performance reasons. The application has two modes, one for Take-off and one for Landing, with the corresponding configuration screens. The procedure screen of both modes has a common structure, and the system displays and announces messages depending on the real time result [21].

Future visions of the system, include the adoption of the aircraft’s center of gravity in the static calculations, so that all factors are taken into consideration when calculating the Accelerate Stop Distance Required and the Landing Distance Required. Furthermore, it is vital that the system is also evaluated for the runway’s slope, an unweighted factor for reasons described in previous sections.

Another considered approach is to reproduce real runway overrun events, in order to assess whether the system and algorithms would be able to notify the pilots in sufficient timing, for the accident to be prevented.

Appendix A.1

The following datarefs were used from the X-Plane 11 simulator, during the development and validation of the system.

Table A1. Datarefs used from the X-Plane 11 simulator.

Field Num.	Field name	Field Description
1	Time	Time since last reset (seconds)
3	Ground Speed	Ground speed (knots)
3	IAS	Indicated Air Speed (knots)
3	EAS	Equivalent Air Speed (knots)
3	TAS	True Air Speed (knots)
4	Speed	Speed (Mach)
4	VSI	Vertical Speed Indicator (ft/min)
4	Acceleration	Longitudinal Acceleration (g)
14	Brakes	Brake usage (ratio 0-1)
18	AoA	AoA (degrees)
20	Longitude	Geographic coordinates (degrees)
20	Latitude	Geographic coordinates (degrees)
21	Distance	Distance covered since last reset (ft)
25	Throttle	Throttle usage ratio (0-1)
35	Thrust	Thrust produced by engines (lbs)
64	Lift	Total lift produced (lbs)
64	Drag	Total drag produced (lbs)
137	Friction force	Friction force with runway (lbs)

Table A1. Datarefs used from the X-Plane 11 simulator.

Field Num.	Field name	Field Description
1	Time	Time since last reset (seconds)
3	Ground Speed	Ground speed (knots)
3	IAS	Indicated Air Speed (knots)
3	EAS	Equivalent Air Speed (knots)
3	TAS	True Air Speed (knots)
4	Speed	Speed (Mach)
4	VSI	Vertical Speed Indicator (ft/min)
4	Acceleration	Longitudinal Acceleration (g)
14	Brakes	Brake usage (ratio 0-1)
18	AoA	AoA (degrees)
20	Longitude	Geographic coordinates (degrees)
20	Latitude	Geographic coordinates (degrees)
21	Distance	Distance covered since last reset (ft)
25	Throttle	Throttle usage ratio (0-1)
35	Thrust	Thrust produced by engines (lbs)
64	Lift	Total lift produced (lbs)
64	Drag	Total drag produced (lbs)
137	Friction force	Friction force with runway (lbs)