Plant Model-Based Fault Detection during Aircraft Takeoff Using Non-Deterministic Finite-State Automata

: In this note, the application of a plant model-based fault detection method for nonlinear control systems on aircraft takeoff is introduced. This method utilizes non-deterministic ﬁnite-state automata, which approximate the fault-free dynamics of the plant. The aforementioned automaton is computed in a preliminary step while during evolution of the plant the automaton is continually evaluated to detect discrepancies between the actual and the nominal dynamics. In this way the fault detection module itself can be implemented on simpler hardware on board of the plant. Moreover, an implementation technique is presented that allows the use of the proposed fault detection method when the plant dynamics is given only by means of a graphical programming script. The great potential and practicality of the used method are demonstrated on a simulated takeoff manoeuvre of a battery-electrically driven aircraft.


Introduction
Fault detection is the first and most important component of the well-known scheme of Fault Detection, Isolation and Accomodation [1,2]. Approaches of fault detection can be model-free (based on simple "if-then" blocks) or model-based (e.g., based on models of signals, the plant dynamics itself or knowledge-based). There is a wide variety of techniques for fault diagnosis. Rather simple limit checking methods are applied in process automation systems [3]. For affine nonlinear systems the work [4] shows a differential geometry approach using residual generators. Non-standard fault detection methods include, for example, the use of artificial neural networks [5] or discrete event systems [6][7][8][9][10][11][12][13][14][15].
Roughly speaking, the present work takes up ideas relying on discrete event systems, yet differs in various aspects. To be specific, this work reformulates some known concepts of model-based fault detection. A substantial difference of this note to existing results is the utilization of a non-deterministic finite-state automaton (FSA) and the particular application on an aircraft takeoff fault detection scenario. The FSA is used to detect discrepancies ("residuals") between the model-predicted and the actual behaviour of the dynamical system under consideration. It approximates the fault-free dynamics in a discrete, non-deterministic form following the methodology of References [16,17]. The non-determinism allows to account for unknown parameters or (regular) disturbances to the plant dynamics.
Classically, this automaton can be obtained fully automated and computationally efficient under certain continuity assumptions on the plant dynamics. For the purpose of fault detection, we propose

Systems and Faults
In this paper we consider continuous-time continuous-state control systems governed by a differential inclusion of the formẋ (t) ∈ F(x(t), u(t)), where x : R → X and u : R → U is the state and input signal, respectively, and where X ⊆ R n , U ⊆ R m . The map F : X × U → P (X) in (1) is set-valued and strict. Disturbances or unknown parameters can be formalized by means of the non-determinism induced by F to the dynamics. For example, a parameter-depend control system governed bẏ where p : R → P, P ⊆ R l is a non-empty set of parameters and f : X × U × P → X is an ordinary map, follows dynamics (1) by setting F = f (·, ·, P).
The method that we are going to use formally applies to sampled versions of the previous dynamics. In fact, dynamics (1) transform through sampling to a discrete-time continuous-state difference inclusion x(k + 1) ∈ G(x(k), u(k)) with the state and input signals x and u, respectively. Here, the time range of x and u is Z + rather than R with discrete-time variable k. We assume the map G : X × U → P (X) to be strict throughout, thus a successor state always exists. G(x 0 , v) is defined as the image of the point x 0 under the flow of the continuous-time dynamics (1) for input signal v at a predefined sampling period τ > 0 [17] (Sec. VIII). The notion of a fault (in the dynamics) is conveniently modeled by a change from G in (3) to some other function so that x(k + 1) / ∈ G(x(k), u(k)) for at least one k ∈ Z + . For the description of the presented method we first formalize the notion of system and its behaviour [16,17].

Definition 1.
A system is a triple (X, U, G) where X and U are non-empty sets and G : X × U → P (X) is strict. The components of the triple are called state space, input space and transition function, respectively. A system whose state space equals R n for some n ≥ 1 is called plant. Definition 2. For a system S = (X, U, G) the behaviour of S is the set In words, the behaviour of a system consists of those pairs of input and output signal that are generated by the system according to (3). See

Fault Detection Method
Below, we introduce discrete approximations [28], also known as discrete abstractions, with which the presented method detects faults. We will use the notion from Reference [17] in a simplified and modified form that is suitable for our purposes.
The key property of a discrete approximation is that it quantizes the continuous state space of a plant into a discrete set and that it "approximates" the behaviour of the plant in an appropriate sense. We continue these descriptive explanations after the precise definition. Definition 3. Let S = (X, U, G) be a plant, let X and U be a cover of X and U, respectively, by non-empty sets. A system S = (X , U , G ) is called discrete approximation of S if holds whenever Ω 1 , Ω 2 ∈ X , Υ ∈ U . Elements of X and U are called cells and input symbols, respectively.
Property (4) is crucial for the "approximation" of the plant behaviour: If for two cells Ω 1 , Ω 2 ∈ X and an input symbol Υ ∈ U it holds Ω 2 / ∈ G (Ω 1 , Υ) then the contrapositive of (4) implies that the plant does not produce signals (u, x) ∈ B(S) such that for some k ∈ Z + . In turn, the observation of a signal possessing the property (5) gives evidence that a fault in the dynamics is present. The latter remark, which is a rigorous mathematical statement, results in a straightforward algorithm for fault detection as follows. Firstly, the components of the algorithm are, besides the discrete approximation (X , U , G ) of the plant (X, U, G), the two quantizers ("analog-to-digital converters") Q state : X → P (X ) and Q input : U → P (U ), (6) which are induced by the covers of X and U, respectively. (I.e., Ω ∈ Q state (a) if and only if a ∈ Ω, and analogously for Q input ). Then, for every point in time k ∈ Z + the values u(k), x(k) and x(k + 1) are measured and converted to Υ(k), Ω(k) and Ω(k + 1), respectively, by means of the quantizers. More concretely, Υ(k) ∈ Q input (u(k)) and Ω(s) ∈ Q state (x(s)) for s ∈ {k, k + 1} is fulfilled. (If the image of a quantizer contains more than one element then an element may be picked randomly for the conversion.) Finally, the test is performed. A fault is signaled if the latter set membership relation is violated. The algorithm is illustrated in Figure 2.
False? Figure 2. Illustration of the fault detection algorithm applied to a plant (X, U, G) with discrete approximation (X , U , G ) and quantizers (6).

Computing Discrete Approximations
Consequently, the presented fault detection method requires a discrete approximation of the plant to monitor.

Computational Method
Methods to compute discrete approximations (for the various existing notions of discrete approximations) have been developed since several years, for example, References [17,19,29]. Typically, the computation of a discrete approximation (X , U , G ) for a given plant (X, U, G) is done in three steps [29], which are outlined subsequently.
Firstly, assuming for simplicity that U is finite, U = U is set and a finite cover X of X is chosen. The cover typically contains two distinguished subsetsX and X \X , which cover the "operating range" of the plant and its complement ("overflows"), respectively. For example, in Figure 3 the cover is formed by the shown squares and outer enclosing rectangles. Secondly, for every cell Ω ∈X the set G(Ω, Υ) is overapproximated. Specifically, a set G(Ω, Υ) that satisfies G(Ω, Υ) ⊆ G(Ω, Υ) is calculated. Here, a proper superset is typically needed; the reachable set G(Ω, Υ) is usually not explicitly representable, whereas an explicitly representable superset might be found. The overapproximation property is also enough in order to establish (4) in the following last step: The automaton is obtained, that is, G is defined, by establishing (4) through intersection tests for every Ω 1 ∈X , Ω 2 ∈ X and Υ ∈ U . To be specific, if (8) is non-empty then Ω 2 ∈ G (Ω 1 , Υ). The trivial choice G (Ω, Υ) = X is made for Ω ∈ X \X ("overflow" cell).

Implementation
In order to allow for efficient computation, the cover in the first step of the scheme is formed by translated copies of a closed hyper-rectangle and a few unbounded hyper-rectangles (aforementioned overflow cells) to fulfill the cover property. When choosing the overapproximating superset G(Ω 1 , Υ) in (8) also as a hyper-rectangle then (8) turns into an intersection test between hyper-rectangles, which can be efficiently evaluated. So the difficult part when implementing the scheme described in Section 2.3.1 for control systems with continuous-time dynamics (1) is to establish a computationally efficient overapproximation method by hyper-rectangles.
System Knowledge Simulation Data

Theory Versus Reality
Such an overapproximating method exists [30] and can be efficiently applied [17,19]. However, one prerequisite is that F in (1) fulfills some smoothness conditions. More generally, the assumption that F in (1) is given as a mathematical formula constitutes a rather ideal situation ("white box" model, see Figure 3) which in applications is rarely present. In contrast, common engineering practice for modelling complex control systems (like aircraft) is the assembling of the dynamics from various subsystems to one huge "simulation model". Quite common is the use of graphical programming tools for it [31]. In this way, however, the function in (1) is hidden behind the assembled interconnection of the subsystems. It is typically unpractical to extract the mathematical formula for further processing, for example, for computing discrete approximations. Therefore, it suggests itself that, with abandoning formal guarantees from theory, empirical techniques may allow to generate an "empirical" discrete approximation from the simulation model. The details to this idea are given below.
Besides the aforementioned case of a "white box" model we further distinguish two more situations. First, the plant dynamics are given bẏ where f is known and possesses the previous continuity properties and where W ⊆ R n is unknown ("gray box" model), but is encoded in the simulation model and 0 ∈ W. (Works like References [17,32,33] model uncertainties or disturbances by W in (9)). In this case, a discrete approximation (X , U , G ) for the nominal dynamics, that is, for the plant (X, U, G) obtained from sampling (9) for W = {0}, is computed in a first step in the convential way. In a second step, transitions are supplemented to the obtained automaton by evaluating the simulation model for a chosen set of test points in X and U. More concretely, for every (Ω, Υ) ∈ X × U two sets of test points x i ∈ Ω, u j ∈ Υ are chosen, where (i, j) ∈ I × J and I and J are finite index sets. Next, the images y i,j of x i under u j obtained from evaluating the simulation model are computed. Then, using the quantizers (6) the corresponding triple T = (Υ, Ω 1 , In this way, transitions due to the initially neglected set W are added to the automaton. Hence, it approximates more accurately the actual dynamics (9) of the plant. See Figure 3b. In the last situation ("black box" model), the simulation model is exclusively used to build up the automaton. As before, for every (Ω, Υ) ∈ X × U two sets of test points x i ∈ Ω, u j ∈ Υ are chosen, where (i, j) ∈ I × J (I and J as before) and the images y i,j of x i under u j obtained from evaluating the simulation model are computed. Finally, the automaton is defined by virtue of See Figure 3c. In the case that the model additionally depends on parameters an additional set of test points d k ∈ P, k ∈ K is chosen, where K is a finite index set. Then, the images y i,j,k of x i under u j for parameter d k are computed and the elements Q state (y i,j,k ) for all (i, j, k) ∈ I × J × K are used in the right hand side of (10) accordingly.
Finally, a remark on input inaccuracies is given. Formally, the control input u(t) in (1) is assumed to be constant during the sampling period. In reality, those control values may vary. These disturbances can be compensated by simply letting neighboring cells of U be strictly overlapping by an amount accounting for the expected input inaccuracies. The previously described computational techniques remain unchanged in this case.

Application: Takeoff Monitoring
To investigate the presented fault detection method in application, we consider a model of a small battery-electrically driven aircraft accelerating on runway until V 1 -speed. At this pre-calculated speed pilots must decide to continue takeoff and lift off or to start braking to abort takeoff. In case of continuation the pilot commands rotation at speed V R ≥ V 1 , that is, initiates the lift-off (cf. Figure 4). We will consider the "black box" situation described in Section 2.4. Nevertheless, for convenience the dynamics of the aircraft are outlined in the first part of this section. Three illustrative fault scenarios will be considered. In addition, this section includes remarks on implementation and finally the experimental evaluation of the presented fault detection algorithm on the chosen example is presented.

Regular Relations during Takeoff
The dynamic model of the battery-electrically driven airplane aims at a detailed energy efficiency modelling of the involved components (Propeller, electric motor including motor controller and battery). Figure 5 shows the components of the propulsion model. Inputs of the system are the aerodynamic velocity V aero , the commanded rotational speed of the propeller N cmd and the air temperature T Air , which is assumed to be constant. Outputs of the system are the thrust T, produced by the propeller, the motor temperature T Mot and the battery state of charge SoC.
The propeller, which produces a thrust T ∈ R, is modelled by means of lookup tables of dimensionless thrust and power coefficients C T (V aero , N) and C P (V aero , N). With these, provided thrust and required power are calculated depending on the aerodynamic velocity V aero and rotational speed N [20]. The motor (including motor controlling) translates the mechanical power P, required by the propeller, into electrical power, provided by the battery, containing motor and controller efficiency η Mot+Cont (N, P) [20]: The battery voltage U Bat is assumed to be constant during acceleration phase. Increasing of the motor temperature T Mot is modelled bẏ with the thermal conductivity of the motor windings. The cooling fluxQ = α cool · S cool (T Mot − T Air ) is determined with heat transfer coefficient α cool , the cooling surface S cool , the thermal conductivity k T Mot and the air temperature T Air , where all terms except the motor temperature T Mot are assumed to be constant during the acceleration phase of the takeoff. Furthermore, the motor dynamics contain an internal saturation, which may reduce the resulting propeller RPM with respect to limitations I Mot (t) ≤ I, P(t) ≤ P of the motor. Here, I Mot and P denotes the motor current and the motor power, respectively, and some bounding constants I ∈ R and P ∈ R. The resulting propeller RPM N(t) equals Finally, the decrease of the state of charge of the battery SoC is calculated with the nominal capacity C nom of the battery:Ṡ For more detailed descriptions of the components the reader is referred to Reference [20]. The aircraft dynamics are given by a parametrized differential equation of the form (2), where f : R 3 × R 2 × P → R 3 . The parameter vector p = [µ V Wind ] in (2) includes the surface friction coefficient µ and the headwind velocity V Wind . Regarding the later fault scenarios, the surface friction coefficient is realistically bounded by µ := [0.005, 0.04] [24,25] and the headwind amounts to V Wind := [0, 2] [ m /s]. More realistic values of the friction coefficient can be determined by friction tests. It is assumed that real friction parameter values do not spread that wide like given in our paper. Furthermore, wind is assumed to hit the aircraft straightly from the front (Aircraft takeoff is typically performed with headwind best possible with regard to the runway direction).
The state vector is given by where the components represent the aircraft kinematic velocity V = V aero + V Wind , the motor temperature T Mot and the state of charge of the battery (SoC), respectively. The control vector is given by where C L and N cmd , respectively, is the commanded lift coefficient and commanded propeller RPM (rotational speed of the propeller).
Aerodynamic lift and drag forces L and D are computed using standard formulas, for example, (15) with drag coefficient C D depending on C L by a polar, constant air density ρ and reference surface of the wing S ref . The kinematic accelerationV is given bẏ with gravitational acceleration g, the mass m of the aircraft and the friction coefficient µ. Therefore, it is assumed that (1) thrust T and drag D are co-linear withV, (2) lift L and weight force mg are co-linear among each other and perpendicular to thrust and drag, and (3) the runway is exactly horizontal. The dynamics for T Mot and SoC taken from Reference [20].

Fault Modelling
Three different types of faults are modelled in the simulation of the takeoff run. For fault 1, which is an exceptionally increasing motor temperature due to an assumed error in the cooling system, factor e 2 ∈ R is added to the state derivativeṪ Mot from time t 2 . This fault will cause an unusual faster increase of the motor temperature T Mot with (11). Fault 2 is an increasing friction due to a problem with the tyres, so a factor e 1 ∈ R is multiplied to the value of the drag D in (15) beginning from time t 1 of the takeoff run simulation time t. As a consequence, the acceleration (16) will decrease. Fault 3 illustrates a short circuit in the powertrain, due to which the battery current I Bat increases. This is realized by a factor e 3 ∈ R, which is multiplied to I Bat from time t 3 . This fault will lead to an unusual steeper decrease of the battery state of charge SoC (12). See Table 1 for the values of the aforementioned constants.

Implementation of the Takeoff Monitoring
The takeoff monitoring system is divided into two parts according to the presented fault detection method. In the first part, a discrete approximation for the model defined in Section 3.1 is computed. The second part contains the "real-time" algorithm for fault detection during takeoff.

Computing the Discrete Approximation
The discrete approximation S = (X , U , G ) for the given plant (aircraft) is determined according to Section 2, where the "black box" model situation is assumed. A detailed description of the construction of S follows. To begin with, all states, controls and parameter ranges are given in  To obtain the transition function G of the discrete approximation S the sampling time for evaluating the black box model is chosen as τ = 1 s and about 1.67 × 10 8 test points are used. These points form a uniform grid on the Cartesian product of the six intervals defined in Table 2.
The calculation is performed using MATLAB [34]. The latter environment provides easy-to-use grid organization using the function ndgrid on the one hand. On the other hand, it provides user-friendly data accessing functions, for example, linear indexing by means of the function sub2ind. Furthermore, using the function parfor from the Parallel Computing Toolbox [34] the computation of the discrete approximation can be parallelized easily.
The computation of the discrete approximation takes 43 h using 6 threads at 3.7 GHz-CPUs. The required memory for the discrete approximation amounts to 330 MB, where indices are saved as 4-byte integers. Finally, approximately 500 MB are needed for the on board implementation of the introduced algorithm.
Of course, the price to pay for using user-friendly programming environments is an increased runtime in comparison with compiled binaries. Nevertheless, we would like to point out that the computation of the discrete approximation is done once per plant, so a higher runtime does not limit applicability.

Simulation of the Takeoff Run
For a simulation of the takeoff run the initial state is assumed with the constant control u(t) = u 0 with u 0 = 0.95 2500 is applied. The friction parameter µ is assumed to be randomly distributed within its bounds given in Table 2 over the runway to allow an investigation of the robustness of the monitoring module with regard to deviation of parameters. For the simulation of fault 2 (drag increase, see Table 1), additionally wind was added to test the systems reliability in presence of an unknown parameter. Here, headwind is modelled in form of a sine function between values 0.1 and 1.9 m /s and a period time of about 15 s.

Implementation of the Fault Detection Algorithm
The algorithm in Figure 2 is implemented as follows. From time 0 s the state x(t) and input u(t) are measured and converted to Q state (x(t)) and Q input (u(t)), respectively, in a time raster of 0.05 s width and saved in a memory. The time raster is defined according to t = k 1 τ + k 2 · 0.05 with k 1 ∈ N and k 2 ∈ {0, . . . , 19}. From time t = τ = 1 s, fault detection begins with evaluating (7) according to Figure 2. (Before the chosen sampling time τ fault detection is not feasible by concept since the state at time τ is located in the future.) The actual states x(t) are converted similarly by (6) and the corresponding images at time t − τ are taken from memory and used in (7) (cf. Figure 2). Cells and input symbols older than t − τ are discarded. The program sequence for automated takeoff monitoring is depicted in Figure 6.
True? FAULT yes no no no yes yes Figure 6. Program Scheme.
The core function (highlighted block) brings about a decision to continue or abort takeoff. With the hardware described in Section 3.3 one automaton evaluation takes 0.034 ms. This comparison in the real-time module can be implemented as a simple index search. The MATLAB function is member can be used to determine whether the index Q state (x(t)) is saved in the image of the transition function for the pair (Q state (x(t − τ)), Q input (u(t − τ))) or not.

Experimental Evaluation
Below, the aforementioned three fault scenarios are simulated and discussed. The Figures 7-9 include the progression of the affected state (thick black line) together with the error flag of the fault detection program. The state progression without fault (gray dotted) is added to the plot as well as the bounds allowed by the discrete approximation (thin black stair lines).
State progression and the fault detection of Fault 1 is shown in Figure 7. The constant e 2 leads to a steeper increasing of the motor temperature beginning at 6 s compared to the fault free progression (gray dotted line).  Fault 2 scenario can be seen in Figure 8. The increased friction slows down the acceleration during takeoff. Consequently, the lines for state progression with and without fault diverge beginning at t 1 = 6 s. At 8.1 s a violation of property (7) occurs for the first time and the fault detection module signals a fault. The resolutions chosen to compute the discrete approximation and the wide ranges of both parameter µ and V Wind do not allow an ealier detection of the fault. For the same reason an alternating 'true-false' fault signal with a gap around 9 and 10 s appears first. From 12.5 s a stable error is detected. A reason for this late fault detection is the rougher resolution of state V used in the discrete approximation in comparison to the state T Mot (cf. Table 2).
As described before, a sinusoidal headwind is added to stress the presented method's reliability. Therefore, the phase shift of the sinusoidal headwind was varied manually, so that the duration between appearance and detection of the fault is nearly maximal. By using this adverse constellation of fault and environment, the reliability of the presented system can be tested. Figure 9. Increasing battery current leads to a steeper decrease of the state of charge in comparison to the non-fault progression beginning at 2 s. Here, no fault can be detected in contrast to the cases of faults 1 and 2. Despite the fault, the SoC stays within the limits recorded in the discrete approximation. The reasons for this is the poor resolution of the cover of the discrete approximation in this component (cf. Table 2) and the comparatively 'slow' reacting state SoC. The scenario of fault 3 shows that the presented fault detection scheme has to be tested extensively for every application and simulation times of the discrete approximation are to be adjusted for proper work. Furthermore, the presented method can not guarantee complete reliability in case of uncertainties in model and measurement.

Discussion and Outlook
A plant model-based fault detection method was presented and applied to an aircraft takeoff scenario. Discrepancies between nominal and actual dynamics are detected utilizing a discrete approximation of the plant dynamics in form of a finite-state automaton. This discrete approximation is computed in a preliminary step. The fault detection module itself is less complex and can be executed on simple hardware.
To show an application of the presented fault detection method, faults during the takeoff run of a battery-electrically driven airplane have been simulated and detection capabilities have been analysed. The application example is based on a "black box" simulation model with which the discrete approximation is generated. It could be shown that with reasonable choices of the covers of the discrete approximation and the resolution of the test points the fault detection system performs reliably. Fault scenario 3 demonstrates that especially for 'slow reacting' states resolution must be enhanced to allow proper fault detection. Furthermore, it could be shown that the separation of the presented fault detection system in two segments simplifies its applicability. In spite of the long preliminary calculation period, the on board fault detection module provides real time capability. For applications with known dynamics (e.g., 'white box'), the presented method is perhaps not the most efficient moreover, but it can be utilized nevertheless.
The presented fault detection system offers further development potential. With small customization of the presented theory the system can be extended to include plant outputs or observable states in the fault monitoring instead of states. Therefore, more system information can be observed and the reliability of the system can be improved. For the demonstrated application, a plant output could be the resulting propeller RPM speed or the battery current. Moreover, the choice of the sampling time (denoted by τ in this note) should be investigated to reduce general computing costs on the one hand and to improve reliability of fault detection on the other hand. A decoupling of different cause-and-effect relationships between faults and states is explicitly not suggested. By using only one FSA for the complete plant dynamics approximation, even faults with an unknown or unexpected cause-and-effect relationship can be detected. In time-and security critical processes, it is not necessary to know the exact cause of a failure. If any system parameter was not regular during the process, it should be aborted and the failure should be reconstructed under laboratory conditions. Author Contributions: Conceptualization, A.W. and F.S.; methodology, A.W.; software, F.S.; validation, F.S.; investigation, A.W. and F.S.; resources, A.K; data curation, A.W. and F.S.; writing-original draft preparation, A.W. and F.S.; writing-review and editing, all authors; visualization, A.W. and F.S.; project administration, A.K.; funding acquisition, A.K. All authors have read and agreed to the published version of the manuscript.