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

Abstract

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 finite-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.

1. 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 a technique for generating an “empirical” approximation in case when only a simulation model, for example, in form of a graphical programming script, is available. Clearly, this is not an exceptional case in industrial engineering environments.

Though our method can be formally used in a wide range of situations (from plants with known dynamics to the aforementioned case of simulation models) it is not the purpose of this paper to compare our method in each application case with well-established fault detection techniques. The goal of this work is to present a fully automated and easy-to-implement technique and we are going to highlight the case of simulation models.

The differences to existing works are various. Firstly, “empirical” approximations have been proposed previously, yet we extend this idea to generate approximations given by non-deterministic finite-state automata rather than directed graphs (deterministic automata) [18]. As for fault-detection based on discrete event systems [6,7,8,9,10,11,12,13,14,15], References [6,7,8,9] build up automata by hand on the basis of discrete events. The authors of References [10,11,12] model the used discrete automaton based on the continuous dynamics of the plant under consideration, as we do. In Reference [15] the original system is mapped to a simplified system and conditions are given under which it is sufficient to perform fault-detection on the simplified model only. The authors of References [13,14] add focus on stochastic automata. In contrast, our approach relies on a fully automated sampling-based methodology involving a non-determinism already mentioned before. The latter has been proven to be applicable in real applications in the context of symbolic controller synthesis, for example, Reference [19].

To demonstrate the practicality this work investigates the application of the novel fault detection method to aeronautics. It is applied to the acceleration phase during aircraft takeoff using a detailed propulsion model of a battery-electrically driven aircraft from Reference [20]. Our motivation origins from the fact that a huge amount of accidents in aviation is related to the first phases of flight. For example, between  1999 and 2002 about $25\%$ of all accidents of airliners happened during takeoff and climb phase [21,22]. The accident of a Boeing 737–800 in Goa is a recent example, where an insufficient stabilization phase for the two engines caused a runway excursion to the side [23]. To minimize risks during acceleration phase, system states have to be monitored carefully by the pilot. Automated fault detection, supporting pilot assistance systems can reduce workload for the pilot, risk potential can be identified earlier and therefore accidents can be avoided.

Several work concentrates on fault detection during the takeoff. References [24,25] show a takeoff monitoring concept, which contains two segments: In the first (pre-takeoff) segment, scheduled performance data is generated, which is compared with measured data in the second (real-time) segment. In contrast to the approach we will present, data measurement is done manually and specific to the problem at hand. In References [26,27] a so-called Flight Safety Assessment and Management is investigated which uses the theory of discrete event systems.

The rest of the paper is organized as follows. First, Section 2 explains the presented fault detection method. It includes details to the concept of discrete approximations of plant dynamics. Section 3 discusses the application of the main contributions to aeronautics in simulation. Section 4 contains a short conclusion and outlook.

2. Modification of Common Methodology

Before presenting the details of the novel method, the basic notation and terminology used subsequently is introduced.

The set of real numbers, non-negative real numbers, integers and non-negative integers is denoted by $\mathbb{R}$, ${\mathbb{R}}_{+}$, $\mathbb{Z}$ and ${\mathbb{Z}}_{+}$, respectively. For a set A we denote by $\mathcal{P}\left(A\right)$ the power set of A. $A\backslash B$ stands for the set difference. A cover of a set A is a subset of $\mathcal{P}\left(A\right)$ such that the union of its elements is equal to A. All elements of the cover represent A in turn. For non-empty sets A and B we denote by $B^A$ the set of all maps $A\to B$. If $A\subseteq \mathbb{R}$ then an element of $B^A$ is called signal with time range A. A set-valued map $f:A\to B$ is called strict if $\varnothing \ne f\left(a\right)$ for all $a\in A$.

2.1. Systems and Faults

In this paper we consider continuous-time continuous-state control systems governed by a differential inclusion of the form

$$\dot{x}\left(t\right)\in F(x\left(t\right),u\left(t\right)),$$

(1)

where $x:\mathbb{R}\to X$ and $u:\mathbb{R}\to U$ is the state and input signal, respectively, and where $X\subseteq {\mathbb{R}}^n$, $U\subseteq {\mathbb{R}}^m$. The map $F:X\times U\to \mathcal{P}\left(X\right)$ 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 by

$$\dot{x}\left(t\right)=f(x\left(t\right),u\left(t\right),p\left(t\right)),$$

(2)

where $p:\mathbb{R}\to P$, $P\subseteq {\mathbb{R}}^l$ is a non-empty set of parameters and $f:X\times U\times P\to 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)\in G\left(x\right(k),u(k\left)\right)$$

(3)

with the state and input signals x and u, respectively. Here, the time range of x and u is ${\mathbb{Z}}_{+}$ rather than $\mathbb{R}$ with discrete-time variable k. We assume the map $G:X\times U\to \mathcal{P}\left(X\right)$ 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 $\tau >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)\notin G\left(x\right(k),u(k\left)\right)$ for at least one $k\in {\mathbb{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\times U\to \mathcal{P}\left(X\right)$ is strict. The components of the triple are called state space, input space and transition function, respectively. A system whose state space equals ${\mathbb{R}}^n$ for some $n\ge 1$ is called plant.

Definition 2.

For a system $S=(X,U,G)$ the behaviour of S is the set

$$\mathcal{B}\left(S\right):=\{(u,x)\in {(U\times X)}^{{\mathbb{Z}}_{+}}\mid {\forall }_{k\in {\mathbb{Z}}_{+}}:\left(3\right)holds\}.$$

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 Figure 1.

2.2. 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^{\prime }$ and $U^{\prime }$ be a cover of X and U, respectively, by non-empty sets. A system $S^{\prime }=(X^{\prime },U^{\prime },G^{\prime })$ is called discrete approximation of S if

$${\Omega }_2\cap G({\Omega }_1,\mathrm{Y})\ne \varnothing \Rightarrow {\Omega }_2\in G^{\prime }({\Omega }_1,\mathrm{Y})$$

(4)

holds whenever ${\Omega }_1,{\Omega }_2\in X^{\prime }$, $\mathrm{Y}\in U^{\prime }$. Elements of $X^{\prime }$ and $U^{\prime }$ are called cells and input symbols, respectively.

Property (4) is crucial for the “approximation" of the plant behaviour: If for two cells ${\Omega }_1$, ${\Omega }_2\in X^{\prime }$ and an input symbol $\mathrm{Y}\in U^{\prime }$ it holds ${\Omega }_2\notin G^{\prime }({\Omega }_1,\mathrm{Y})$ then the contrapositive of (4) implies that the plant does not produce signals $(u,x)\in \mathcal{B}\left(S\right)$ such that

$$(u\left(k\right),x\left(k\right),x(k+1))\in \mathrm{Y}\times {\Omega }_1\times {\Omega }_2$$

(5)

for some $k\in {\mathbb{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^{\prime },U^{\prime },G^{\prime })$ of the plant $(X,U,G)$, the two quantizers (“analog-to-digital converters”)

$$Q_{\mathrm{state}}:X\to \mathcal{P}\left(X^{\prime }\right)\mathrm{and}{Q}_{\mathrm{input}}:U\to \mathcal{P}\left(U^{\prime }\right),$$

(6)

which are induced by the covers of X and U, respectively. (I.e., $\Omega \in Q_{\mathrm{state}}\left(a\right)$ if and only if $a\in \Omega$, and analogously for $Q_{\mathrm{input}}$). Then, for every point in time $k\in {\mathbb{Z}}_{+}$ the values

$$u\left(k\right),x\left(k\right)\mathrm{and}x(k+1)$$

are measured and converted to

$$\mathrm{Y}\left(k\right),\Omega \left(k\right)\mathrm{and}\Omega (k+1),$$

respectively, by means of the quantizers. More concretely, $\mathrm{Y}\left(k\right)\in Q_{\mathrm{input}}\left(u\left(k\right)\right)$ and $\Omega \left(s\right)\in Q_{\mathrm{state}}\left(x\left(s\right)\right)$ for $s\in \{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

$$\Omega (k+1)\in G^{\prime }(\Omega \left(k\right),\mathrm{Y}\left(k\right))$$

(7)

is performed. A fault is signaled if the latter set membership relation is violated. The algorithm is illustrated in Figure 2.

2.3. Computing Discrete Approximations

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

2.3.1. 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^{\prime },U^{\prime },G^{\prime })$ 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^{\prime }=U$ is set and a finite cover $X^{\prime }$ of X is chosen. The cover typically contains two distinguished subsets ${\overline{X}}^{\prime }$ and $X^{\prime }\backslash {\overline{X}}^{\prime }$, 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 $\Omega \in {\overline{X}}^{\prime }$ the set $G(\Omega ,\mathrm{Y})$ is overapproximated. Specifically, a set $\widehat{G}(\Omega ,\mathrm{Y})$ that satisfies $G(\Omega ,\mathrm{Y})\subseteq \widehat{G}(\Omega ,\mathrm{Y})$ is calculated. Here, a proper superset is typically needed; the reachable set $G(\Omega ,\mathrm{Y})$ 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^{\prime }$ is defined, by establishing (4) through intersection tests

$${\Omega }_2\cap \widehat{G}({\Omega }_1,\mathrm{Y})$$

(8)

for every ${\Omega }_1\in {\overline{X}}^{\prime }$, ${\Omega }_2\in X^{\prime }$ and $\mathrm{Y}\in U^{\prime }$. To be specific, if (8) is non-empty then ${\Omega }_2\in G^{\prime }({\Omega }_1,\mathrm{Y})$. The trivial choice $G^{\prime }(\Omega ,\mathrm{Y})=X^{\prime }$ is made for $\Omega \in X^{\prime }\backslash {\overline{X}}^{\prime }$ (“overflow” cell).

2.3.2. 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 $\widehat{G}({\Omega }_1,\mathrm{Y})$ 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.

2.4. 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 by

$$\dot{x}\left(t\right)\in f(x\left(t\right),u\left(t\right))+W$$

(9)

where f is known and possesses the previous continuity properties and where $W\subseteq {\mathbb{R}}^n$ is unknown (“gray box” model), but is encoded in the simulation model and $0\in W$. (Works like References [17,32,33] model uncertainties or disturbances by W in (9)). In this case, a discrete approximation $(X^{\prime },U^{\prime },G^{\prime })$ for the nominal dynamics, that is, for the plant $(X,U,G)$ obtained from sampling (9) for $W=\left\{0\right\}$, 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 $(\Omega ,\mathrm{Y})\in X^{\prime }\times U^{\prime }$ two sets of test points $x_i\in \Omega$, $u_j\in \mathrm{Y}$ are chosen, where $(i,j)\in I\times 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=(\mathrm{Y},{\Omega }_1,{\Omega }_2)$ satisfying $(u_j,x_i,y_{i,j})\in T$ is determined. If ${\Omega }_2\notin G^{\prime }({\Omega }_1,\mathrm{Y})$ then $G^{\prime }$ is redefined so that ${\Omega }_2\in G^{\prime }({\Omega }_1,\mathrm{Y})$. 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 $(\Omega ,\mathrm{Y})\in X^{\prime }\times U^{\prime }$ two sets of test points $x_i\in \Omega$, $u_j\in \mathrm{Y}$ are chosen, where $(i,j)\in I\times 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

$$G^{\prime }(\Omega ,\mathrm{Y})=\{Q_{\mathrm{state}}\left(y_{i,j}\right)\mid (i,j)\in I\times J\}.$$

(10)

See Figure 3c. In the case that the model additionally depends on parameters an additional set of test points $d_k\in P$, $k\in 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_{\mathrm{state}}\left(y_{i,j,k}\right)$ for all $(i,j,k)\in I\times J\times K$ are used in the right hand side of (10) accordingly.

Finally, a remark on input inaccuracies is given. Formally, the control input $u\left(t\right)$ 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^{\prime }$ be strictly overlapping by an amount accounting for the expected input inaccuracies. The previously described computational techniques remain unchanged in this case.

3. 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\ge 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.

3.1. 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_{\mathrm{aero}}$, the commanded rotational speed of the propeller $N_{\mathrm{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_{\mathrm{Mot}}$ and the battery state of charge $\mathrm{SoC}$.

The propeller, which produces a thrust $T\in \mathbb{R}$, is modelled by means of lookup tables of dimensionless thrust and power coefficients $C_T(V_{\mathrm{aero}},N)$ and $C_P(V_{\mathrm{aero}},N)$. With these, provided thrust and required power are calculated depending on the aerodynamic velocity $V_{\mathrm{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 ${\eta }_{\mathrm{Mot}}+\mathrm{Cont}(N,P)$ [20]:

$$P=U_{\mathrm{Bat}}{I}_{\mathrm{Bat}}{\eta }_{\mathrm{Mot}}+\mathrm{Inv}.$$

The battery voltage $U_{\mathrm{Bat}}$ is assumed to be constant during acceleration phase. Increasing of the motor temperature $T_{\mathrm{Mot}}$ is modelled by

$${\dot{T}}_{\mathrm{Mot}}=\frac{d{T}_{\mathrm{Mot}}}{dt}=\left(P\left(1-{\eta }_{\mathrm{Mot}}+\mathrm{Cont}\right)-{\dot{Q}}_{\mathrm{cool}}\right)\frac{1}{k_{T_{\mathrm{Mot}}}}$$

(11)

with the thermal conductivity of the motor windings. The cooling flux $\dot{Q}={\alpha }_{\mathrm{cool}}·S_{\mathrm{cool}}\left(T_{\mathrm{Mot}}-T_{\mathrm{Air}}\right)$ is determined with heat transfer coefficient ${\alpha }_{\mathrm{cool}}$, the cooling surface $S_{\mathrm{cool}}$, the thermal conductivity $k_{T_{\mathrm{Mot}}}$ and the air temperature $T_{\mathrm{Air}}$, where all terms except the motor temperature $T_{\mathrm{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 $\mathrm{RPM}$ with respect to limitations $I_{\mathrm{Mot}}\left(t\right)\le \overline{I}$, $P_{}\left(t\right)\le \overline{P}$ of the motor. Here, $I_{\mathrm{Mot}}$ and $P_{}$ denotes the motor current and the motor power, respectively, and some bounding constants $\overline{I}\in \mathbb{R}$ and $\overline{P}\in \mathbb{R}$. The resulting propeller $\mathrm{RPM}$ $N\left(t\right)$ equals

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{N}_{\mathrm{cmd}}\hfill & \mathrm{if}{I}_{\mathrm{Mot}}(N_{\mathrm{cmd}},t)<\overline{I}\wedge P_{}(N_{\mathrm{cmd}},t)<\overline{P}\hfill \\ N\left(I_{\mathrm{Mot},\max }\right)\hfill & \mathrm{if}{I}_{\mathrm{Mot}}(N_{\mathrm{cmd}},t)\ge \overline{I}\wedge P_{}(N_{\mathrm{cmd}},t)<\overline{P}\hfill \\ N\left(P_{\max }\right)\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

Finally, the decrease of the state of charge of the battery $\mathrm{SoC}$ is calculated with the nominal capacity $C_{\mathrm{nom}}$ of the battery:

$$\dot{SoC}=\frac{d\mathrm{SoC}}{dt}=\frac{I_{\mathrm{Bat}}}{C_{\mathrm{nom}}}.$$

(12)

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:{\mathbb{R}}^3\times {\mathbb{R}}^2\times P\to {\mathbb{R}}^3$. The parameter vector $p={\left[\mu V_{\mathrm{Wind}}\right]}^{\top }$ in (2) includes the surface friction coefficient $\mu$ and the headwind velocity $V_{\mathrm{Wind}}$. Regarding the later fault scenarios, the surface friction coefficient is realistically bounded by $\mu :=\left[0.005,0.04\right]$ [24,25] and the headwind amounts to ${\mathrm{V}}_{\mathrm{Wind}}:=\left[0,2\right]\left[\mathrm{m}/\mathrm{s}\right]$. 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

$$x\left(t\right)={\left[\begin{array}{ccc}V\left(t\right)& T_{\mathrm{Mot}}\left(t\right)& \mathrm{SoC}\left(t\right)\end{array}\right]}^{\top }$$

(13)

where the components represent the aircraft kinematic velocity $V=V_{\mathrm{aero}}+V_{\mathrm{Wind}}$, the motor temperature $T_{\mathrm{Mot}}$ and the state of charge of the battery ($\mathrm{SoC}$), respectively. The control vector is given by

$$u\left(t\right)={\left[\begin{array}{cc}{C}_L\left(t\right)& N_{\mathrm{cmd}}\left(t\right)\end{array}\right]}^{\top }$$

(14)

where $C_L$ and $N_{\mathrm{cmd}}$, respectively, is the commanded lift coefficient and commanded propeller $\mathrm{RPM}$ (rotational speed of the propeller).

Aerodynamic lift and drag forces L and D are computed using standard formulas, for example,

$$D=C_D\frac{\rho }{2}{V}_{\mathrm{aero}}^2{S}_{\mathrm{ref}}$$

(15)

with drag coefficient $C_D$ depending on $C_L$ by a polar, constant air density $\rho$ and reference surface of the wing $S_{\mathrm{ref}}$. The kinematic acceleration $\dot{V}$ is given by

$$\dot{V}=\frac{T-D-\mu (mg-L)}{m}$$

(16)

with gravitational acceleration g, the mass m of the aircraft and the friction coefficient $\mu$. Therefore, it is assumed that (1) thrust T and drag D are co-linear with $\dot{V}$, (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_{\mathrm{Mot}}$ and $\mathrm{SoC}$ taken from Reference [20].

3.2. 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\in \mathbb{R}$ is added to the state derivative ${\dot{T}}_{\mathrm{Mot}}$ from time $t_2$. This fault will cause an unusual faster increase of the motor temperature $T_{\mathrm{Mot}}$ with (11). Fault 2 is an increasing friction due to a problem with the tyres, so a factor $e_1\in \mathbb{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_{\mathrm{Bat}}$ increases. This is realized by a factor $e_3\in \mathbb{R}$, which is multiplied to $I_{\mathrm{Bat}}$ from time $t_3$. This fault will lead to an unusual steeper decrease of the battery state of charge $\mathrm{SoC}$ (12). See

Table 1. Numerical values and occurrence times of the faults.

Fault i	Increasing of …	${\mathit{e}}_{\mathit{i}}$	${\mathit{t}}_{\mathit{i}}\left[\mathbf{s}\right]$
1	motor temperature ${\dot{T}}_{\mathrm{Mot}}$	2.5	6
2	drag D	2.5	7
3	battery current $I_{\mathrm{Bat}}$	2	2

for the values of the aforementioned constants.

3.3. 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.

3.3.1. Computing the Discrete Approximation

The discrete approximation $S^{\prime }=(X^{\prime },U^{\prime },G^{\prime })$ 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^{\prime }$ follows. To begin with, all states, controls and parameter ranges are given in

Table 2. Ranges of the states, controls and the parameter assumed for the construction of the discrete approximation.

		Unit	Interval	Resolution
States
Velocity	V	$\mathrm{m}/\mathrm{s}$	$\left[0,25\right]$	0.5
Motor temp.	$T_{\mathrm{Mot}}$	°C	$\left[35,45\right]$	0.5
State of charge	$\mathrm{SoC}$	−	$\left[0.94,0.96\right]$	0.004
Controls
Lift coefficient	$C_L$	−	$\left[0.9,1.1\right]$	0.067
Propeller RPM	$N_{\mathrm{cmd}}$	RPM	$\left[2400,2600\right]$	20
Parameter
Friction coeff.	$\mu$	−	$[0.005,0.04]$	0.035
Headwind	${\mathrm{V}}_{\mathrm{Wind}}$	$\mathrm{m}/\mathrm{s}$	$[0,2]$	2

. The cover $X^{\prime }$ of ${\mathbb{R}}^3$ is formed by translated copies of $\left[0,0.5\right]\times \left[0,0.5\right]\times \left[0,0.004\right]$ (covering the operating range $\left[0,25\right]\times \left[35,45\right]\times \left[0.94,0.96\right]$) and some unbounded hyper-rectangles to cover ${\mathbb{R}}^3$. The cover $U^{\prime }$ of U is given by translated copies of $\left[0,0.067\right]\times \left[0,20\right]$. The shape of the cells and input symbols can be deduced from Table 2. The cover $P^{\prime }$ of the parameter range is given by the single hyper-rectangle $\left[0.005,0.04\right]\times \left[0,2\right]$.

To obtain the transition function $G^{\prime }$ of the discrete approximation $S^{\prime }$ the sampling time for evaluating the black box model is chosen as $\tau =1\mathrm{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.7GHz-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.

3.3.2. Simulation of the Takeoff Run

For a simulation of the takeoff run the initial state

$$x\left(0\right)={\left[\begin{array}{ccc}0& 35& 0.95\end{array}\right]}^{\top }$$

is assumed with the constant control $u\left(t\right)=u_0$ with

$$u_0={\left[\begin{array}{cc}0.95& 2500\end{array}\right]}^{\top }$$

is applied. The friction parameter $\mu$ 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\mathrm{m}/\mathrm{s}$ and a period time of about $15\mathrm{s}$.

3.3.3. Implementation of the Fault Detection Algorithm

The algorithm in Figure 2 is implemented as follows. From time $0\mathrm{s}$ the state $x\left(t\right)$ and input $u\left(t\right)$ are measured and converted to $Q_{\mathrm{state}}\left(x\left(t\right)\right)$ and $Q_{\mathrm{input}}\left(u\left(t\right)\right)$, 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\tau +k_2·0.05$$

with $k_1\in \mathbb{N}$ and $k_2\in \{0,\dots ,19\}$. From time $t=\tau =1\mathrm{s}$, fault detection begins with evaluating (7) according to Figure 2. (Before the chosen sampling time $\tau$ fault detection is not feasible by concept since the state at time $\tau$ is located in the future.) The actual states $x\left(t\right)$ are converted similarly by (6) and the corresponding images at time $t-\tau$ are taken from memory and used in (7) (cf. Figure 2). Cells and input symbols older than $t-\tau$ are discarded. The program sequence for automated takeoff monitoring is depicted in Figure 6.

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\mathrm{m}\mathrm{s}$. 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_{\mathrm{state}}\left(x\left(t\right)\right)$ is saved in the image of the transition function for the pair $(Q_{\mathrm{state}}\left(x(t-\tau )\right),Q_{\mathrm{input}}\left(u(t-\tau )\right))$ or not.

3.4. Experimental Evaluation

Below, the aforementioned three fault scenarios are simulated and discussed. The Figure 7, Figure 8 and Figure 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).

Here, the fault can be detected very early and clear at about half a second after the fault appears. The fault detection is ‘disabled’ at about 12 s (red line) since the ground speed V exceeds the upper limit of the operating range (cf. Table 2). In other words, the conversion of the continuous state signal returns overflow cells. Furthermore, it should be noted that the state progress without fault (gray dotted line) does not run within the thresholds included by the FSA. The state progression and the fault detection logic are simulated with the active fault for Figure 7.

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\mathrm{s}$. At $8.1\mathrm{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 $\mu$ and $V_{\mathrm{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\mathrm{s}$ appears first. From $12.5\mathrm{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_{\mathrm{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.

State progression and the fault detection of Fault 3 is shown in 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\mathrm{s}$.

Here, no fault can be detected in contrast to the cases of faults 1 and 2. Despite the fault, the $\mathrm{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 $\mathrm{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.

4. 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 $\mathrm{RPM}$ speed or the battery current. Moreover, the choice of the sampling time (denoted by $\tau$ 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.