Robust fault tolerant control allocation for a modern over‐actuated commercial aircraft

Correspondence Christopher Edwards, College of Engineering, Mathematics and Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK. Email: c.edwards@exeter.ac.uk Abstract This paper presents a novel form of control allocation, designed within a sliding mode framework, for the fault tolerant control of over-actuated systems. The control allocation is designed in such a way as to allow a subset of the actuators to remain inactive under nominal fault-free conditions. In the event that the active set of actuators becomes unable to provide the desired performance, an adaption process takes place which allows the inactive actuators to compensate. A computationally light gradient descent algorithm is proposed to govern the adaption which guarantees that, if possible, actuator saturation is avoided and system performance is maintained, even in the event of severe actuator faults/failures. Rigorous conditions are derived, in terms of the faults/failures, uncertainties in fault reconstruction information and the adaptive process, which ensure sliding occurs in a finite time and that the resulting motion is stable. To demonstrate the effectiveness of the control scheme, a high-fidelity blended wing body aircraft model is also proposed in this paper; this particular configuration of aircraft is nominally unstable, with poor control authority and a large amount of redundancy, making it a suitable candidate for testing reconfigurable fault tolerant control laws in the presence of input constraints.

priority groups; then in the event that one group begins to saturate, the next group in the "chain" is used. This is a popular methodology in the aerospace community. Recently, an adaptive law was proposed in [8] which allows a suitable CA law to be calculated on-line, encompassing actuator non-linearities. Many CA methods utilise on-line optimisation to ensure the control distribution satisfies some criteria. In [9], a parameter within the CA scheme is chosen through an unconstrained optimisation process to minimise any saturation in the system. Contrastingly, quadratic programming is explored in [10] and [11], where the difference in the achieved and commanded "virtual" control is minimised, using the saturation limits as constraints. In the event that redistribution of the "virtual" control alone is not sufficient to prevent saturation, [12] proposes an optimisation scheme which minimises any degradation to the system's performance through optimising a Lyapunov function.
BWB aircraft are a natural candidate for demonstrating the efficacy of CA schemes because of their unique actuator configuration. The following works have explicitly considered the use of CA for BWBs. In [13], an analysis of the coupling between the BWB's control surfaces is presented. A discussion on the nonlinear interactions between the control surfaces, and how CA can be used to accommodate them, can be found in [11,14]. To help alleviate the effects of coupling, the work in [15] proposes weighting the lateral and longitudinal CA modules differently, allowing different actuators to be prioritised for specific roles. In the event of actuator faults and failures, these roles are relaxed avoiding any need to reconfigure the controller. The work in [16] proposes a method of designing the CA to minimise the trim-drag of the aircraft.
The main focus of this paper however is to consider the effects of actuator faults and failures. One of the main requirements for fault tolerant control is robustness to both uncertainties and disturbances. Multiple levels of Sliding Mode Controllers (SMCs) are designed for the control of a tailless aircraft in [17] and [18]: one SMC provides a robust "virtual" control law which is then distributed through a second SMC to the actuator suite, ensuring that the "virtual" control is robustly achieved. Works which have considered the use of robust control in conjunction with CA include [19] and [20]. In [19] and [20], H ∞ control and back-stepping are used, respectively, whereas SMC is used in [3] and [21] for FTC of aerospace systems. Higherorder SMCs (HOSMCs) have also been considered for FTC. These feature continuous control laws and therefore avoid the phenomenon of "chattering," which can be problematic for conventional SMC. In [22], an HOSMC is compared with a conventional SMC for the control of a fighter aircraft subject to different uncertainties and disturbances -including actuator faults -demonstrating the HOSMC's ability to alleviate the "chattering" experienced by the conventional SMC scheme. An adaptive HOSMC is proposed in [23] which does not require prior knowledge of the bounds on any disturbances/uncertainties to ensure that sliding is induced in finite time. In [24], Non-Linear Dynamic Inversion is used in conjunction with an HOSMC to further reduce the effects of uncertainty in the system. Some robust control methods have been designed to prevent saturation, without the need for redistribution. An Integral SMC scheme is presented in [25] which blends high-and low-performing controllers; an optimisation scheme chooses the best possible performing controller that keeps actuators within their limits. In a similar vein, [26] proposes a nominal reference model which is adapted on-line to change the system's performance. The closed-loop performance and the control distribution are designed through a constrained optimisation as part of a Model Predictive Control (MPC) law in [27]. However, due to the complexity of the optimisation, MPC is computationally expensive, making it difficult to accomplish in real time. In [5] and [28], variations of Model Reference Adaptive Control (MRAC) are proposed, which, although do not prevent saturation, provide some robustness against saturation. These schemes are shown to be effective in controlling a flying wing type aircraft (similar to the BWB). In [5] a "minimal control" synthesis is presented which does not rely on the identification of model parameters, whereas in [28], the reference model is of a reduced order, based on the tracking error system.
In this paper, a novel fault tolerant SMC and an adaptive CA scheme are developed. As compared to the work in [3] (where the CA is either fixed or dependent on the actuators' health), the proposed scheme in this paper is more elaborate and includes an additional weighting in the CA which provides further design freedom to prioritise the use of specified actuators. Furthermore, although the scheme in [3] can handle actuator faults/failures, the scheme does not consider actuator saturation: this paper proposes a scheme which allows the control effort to be redistributed in such a way as to ensure, if possible, that saturation is avoided. Another of the paper's main contribution is the modification of the CA structure in [9], which allows the coupling in a BWB's actuator suite to be more effectively handled. This is achieved through allowing a subset of the actuator suite to produce the entire control effort in nominal conditions. If this subset begins to saturate, then an adaptive process allows the entire actuator suite to be used. The SMC is designed to be robust to uncertainties in the fault reconstruction (which is used in the CA).
The structure of the paper is as follows: Section 2 presents an overview of the problem addressed in this paper. In Sections 3 and 4, the novel CA scheme is proposed and then appropriate SMC laws are designed and analysed. Section 5 provides the details of a high-fidelity 6 Degree of Freedom (DoF) non-linear BWB model which will form the basis of the control design, which is discussed in Section 6, and the simulation results, which are shown in Section 7. Section 8 provides some concluding remarks.
The notation used in this paper is quite standard: ℝ represents the set of all real numbers and I represents an appropriately sized identity matrix. Unless specified otherwise, ‖ ⋅ ‖ represents the Euclidean 2 norm. The vector x = col (x 1 , … , x m ) represents the concatenation of the components into a column vector and X = diag(x 1 , .., x m ) represents a matrix whose diagonal elements are the x i s.

PROBLEM STATEMENT
Consider an n th order over-actuated Linear Time Invariant (LTI) systemẋ where the system's state, input and the controlled outputs are, respectively, denoted by x(t ) ∈ ℝ n , u(t ) ∈ ℝ m and y(t ) ∈ ℝ l . Here, it is assumed that all the states are available. In this paper, it is also assumed that actuator redundancy is available, that is the system is over-actuated (m > l ), meaning that there are strictly more actuators than outputs to be controlled. The corresponding system matrices are: is a time-varying matrix, which represents any actuator faults and failures present in the system. Each scalar diagonal component w f i (t ) signifies the health of the i th actuator and as such is bounded so that 0 ≤ w f i (t ) ≤ 1. A healthy actuator is indicated by w f i = 1, whilst w f i = 0 represents a completely failed actuator-anywhere in between these values indicates a fault. In reality, the health of the actuators are not precisely known. For practical purposes, they can be approximated online by a Fault Detection and Isolation (FDI) system (see, e.g. the sliding mode approach proposed in [29] or the Kalman filter approach in [30]); however, this introduces uncertainty into the feedback system. Therefore, consider the following expansion of W f (t ) in the state equations of (1) so thaṫ whereŴ f (t ) ∈ ℝ m×m represents the approximation of W f (t ) (as calculated by an FDI scheme) and Δ(t ) ∈ ℝ m×m is a diagonal matrix denoting the uncertainty in the fault reconstruction.
To capture the effects of actuator saturation, it is assumed that the individual components of the control signal in (2), u(t ), are bounded satisfying where u c (t ) is the signal generated by the controller and u(t ) is the achieved control signal (i.e. that which is seen by the plant). The maximum and minimum limits of the i th actuator are denoted by u i max and u i min , respectively; typically in aerospace systems, these are known and represent the minimum and maximum deflections of flaps or the minimum and maximum thrust produced by an engine.

CONTROL ALLOCATION
Consider the system presented in (2). For the purpose of control design, assume the states can be reordered and partitioned to achieve the following form: where B 2 ∈ ℝ l ×m contains the majority of the control effort and B 1 ∈ ℝ (n−l )×m is considered "small" in comparison [3]. For the sake of simplicity, it is assumed that ‖B 2 ‖ = 1: such a result is always possible by scaling the last l states of the system. Commonly, in the CA literature, a reduced order "virtual" control signal v(t ) ∈ ℝ l is related to the physical control signal u c (t ) through a pseudo-inverse structure. In [9,31], the following form is proposed: where W p1 ∈ ℝ m×m is a fixed diagonal matrix which constitutes design freedom, thus allowing the use of certain actuators to be prioritised. Through direct evaluation, it can be verified that this choice of u c (t ) satisfies and therefore demonstrates that, despite actuator faults and failures, v(t ) is mapped to the bottom partition of (2) and directly represents the effect of the control on the last l channels ofẋ(t ): as long as there is no uncertainty in the system. Although it is possible to use W p1 in (5) to try to accommodate saturation limits, it is not straight-forward as the relationship between W p1 and u c (t ) is non-linear. This approach would also remove any facility to prioritise actuators subject to other criteria, such as (in an aerospace context) reducing structural loading or increasing efficiency. As an alternative, consider the addition of an adaptive component to (5): specifically where Θ(t ) = diag( 1 (t ), … , m (t )) is an adaptive parameter, B † 2 is the pseudo inverse of B 2Ŵ f defined in (5) and B ‡ 2 ∈ ℝ m×l is chosen as any appropriate right pseudo inverse of B 2Ŵ f (this will be discussed later). It can be verified that, independently of Θ(t ) and B ‡ 2 (t ), (7) satisfies (6) and that u c (t ) varies linearly with respect to Θ(t ). Furthermore, when Θ(t ) = 0, the expression (7) reduces to (5), thus retaining the original CA form under nominal conditions. A suitable choice of Θ(t ) is discussed in Section 7.2. Remark 1. The choice of structure in (7) is inspired by the adaptive formulation in [6]. This has been modified in [9] to accommodate actuator faults and failures. In [9], the pseudo-inverse B ‡ 2 (t ) is chosen such that B ‡ 2 (t ) = B † 2 (t ) and a full-rank restriction is placed on the design of W p1 . In this paper, B ‡

(t ) is chosen as
where W p2 ∈ ℝ m×m is a second priority weighting matrix. By Remark 2. If the restriction on W p1 in [9] is relaxed, (5) can be designed to operate within a subset of the actuator suite since the i th actuator can be removed from the distribution of (5) by setting W p1 i = 0. This may be desirable in cases where actuators should be reserved for their primary functions which do not fit in the controllers remit; for example, in aerospace systems the main purpose of engines are for longitudinal control but they could also be used for lateral control through differential thrust.
In the event of saturation, the adaptive parameter Θ(t ) redistributes the control effort through the adaptive component of (7) and this can be tuned through choice of W p2 in (8), to utilise a different set of actuators. In this paper, the entire actuator suite is considered for adaption by making W p2 full rank. This allows the CA structure in (7) to use a reduced set of actuators in nominal conditions (when Θ(t ) = 0) whilst allowing the entire actuator suite to be used in the event of saturation (by setting Θ(t ) ≠ 0).

Consider a set of allowable actuator faults and failures
Remark 3. Within the set , the norms of the pseudoinverse structures B † 2 and B ‡ 2 , defined in (5) and (8), respectively, are bounded. This is due to the boundedness properties of pseudo-inverses described in [32]. Therefore, the following bounds are always guaranteed to exist. The smallest possible values of 0 and 1 can be calculated efficiently through the use of Genetic Algorithms to search the entire set  (which is compact). The gains 0 and 1 will be exploited in the stability analysis later in the paper.

OPTIMISATION ALGORITHM
This section focusses on the choice of Θ(t ) in (7): specifically an adaptive algorithm will be proposed suitable for running in real time.

Cost function
For any given actuator suite, there exists a set of "virtual" control signals, v(t ) ∈ , which is achievable subject to the actuator suite's health and saturation limits. Assumption: The virtual control v(t ) always remains in the set ; in other words, there must always exist a value of Θ(t ) = Θ * which, when substituted into (7), ensures the saturation error is equal to zero. The problem considered here is one of perfectly allocating v(t ), that is, finding the unknown parameter Θ * . Substituting from (7), into (10), for a generic choice of Θ(t ), and for where ) .
To develop an algorithm to find Θ * on-line, define the n th estimate of Θ * , at a given time t , as Θ n , and define the error function (10) when evaluated at Θ n as e n = Λ(Θ n − Θ * ).
In (13), Λ is considered fixed and represents the value of Λ(t ) in (12) at the current instant of time. Assuming that the saturation limits in (3) are known, e n in (13) can be evaluated on-line (without knowing Θ * ) and therefore can be used to iteratively find Θ * . Furthermore, note that since rank(Λ) ≤ m − l, e = 0 ⇏ = * . To develop a iterative scheme to ensure e n → 0, consider the quadratic cost function It is clear that f (Θ n ) ≥ 0 and for any v(t ) ∈ , the minimum of (14) is at f (⋅) = 0. Note that f (Θ n ) = 0 ⇏ Θ n = Θ * but it does ensure that u c (t ) = u(t ) and hence v(t ) is perfectly allocated.

Gradient decent method
To find the minima of the function f (⋅), the following Gradient Descent algorithm, based on the work in [9], is adopted. Consider the following update law for selecting Θ n+1 in terms of the previous estimate Θ n : where the scalar a n represents the step length whilst P n ∈ ℝ m represents the search direction [33]. The step length is chosen to minimise f (Θ n+1 ) = f (Θ n + a n P n ) at each iteration. In such algorithms, the search direction is commonly chosen as the steepest negative gradient: here this is given by where Θ n i represents the i th diagonal element of Θ n . The value of ∇ f (Θ n ) can be explicitly calculated as from the definition of f ( n ) in (14). Evaluating f ( n+1 ) in (14), and utilising (15), yields It is clear that minimising f ( n+1 ) can be achieved by minimising the term f n . Using the definitions of n from (15) and P n from (16) and (17), f n can be expressed as Thinking of the right-hand side of (19) as a quadratic scalar equation, in terms of a n , it is clear that f n has a minimum at Theorem 1. The scalar a n , as defined in (20) always remains bounded and so by updating Θ n+1 according to where P n are chosen as described in (17), will always result in a decrease of the cost-function (14).
Proof. The proof is similar to Proposition 1 in [9]. The algorithm is guaranteed to converge to a local minimum where f (Θ n ) = 0, but does not guarantee Θ n → Θ * . □ Remark 4. In Section 7.2, the computational load of the algorithm is discussed; it is demonstrated that the algorithm is sufficiently "light" to run in real time on common, existing employed micro-processors.

SLIDING MODE CONTROL DESIGN
In this section, a SMC procedure will be used to design a robust control law for v(t ) to control (2) through the CA structure (7).

Sliding surface stability analysis
Assuming there is no saturation in the system (u(t ) = u c (t )), substituting (7) into (2) yields the following faulty "virtual" system:ẋ Define a coordinate transformation x ↦ T r x =x, where then in the new coordinates the system in (22) becomeṡx Under "ideal" operating conditions, that is, where the actuators are all healthy and perfectly monitored ( , and there is no adaption (Θ(t ) = 0); it can be verified through direct evaluation that which means (24) reduces to the regular form [34]: This system will be used for the basis of the control design. Assuming that the system pair (Â,B v ) is controllable, (Â 11 ,Â 12 ) is also controllable and a reduced order switching function can be designed [34]. Consider the hyperplane defined by where the switching function s(t ) :∈ ℝ n →∈ ℝ l is given by Here, the matrix M ∈ ℝ l ×(n−l ) represents the design freedom: a minimal requirement being that M must be chosen so that To analyse the dynamics of (24), when confined to the hyperplane , define a second coordinate trans- In the new coordinate system, (24) can be rewritten as and where In particular, from (31),Ã 11 =Â 11 −Â 12 M and is therefore Hurwitz stable. When confined to the hyperplane , the system is said to be in an ideal sliding motion [34]. This is characterised bẏs(t ) = s(t ) = 0. Substituting these values into the bottom partition of (31) and solving for v(t ) yields the control effort necessary to maintain sliding. This gives the so-called equivalent control, v eq (t ), which in this case can be calculated as Substituting (34) for v(t ) in the top partition of (31) yieldṡx which describes the equations of motion governing sliding. The reduced order dynamic (35) can be thought of as the negative feedback interconnection of the non-linear feedback gaiñ and the LTI planṫx which has the associated transfer functioñ For the subsequent analysis define Since M in (29) is chosen to stabiliseÃ 11 , it follows thatG (s) is stable and therefore 3 is finite.

Theorem 2. The sliding motion presented in (35) is asymptotically stable for any fault
where 0 and 1 are defined in (9), 2 is defined in (39), 3 is defined in (40) and Θ max is a bound on the adaptive matrix such that Θ max > ‖Θ(t )‖, where Θ(t ) is defined in (7).
Proof. For the equivalent control in (34) to be unique, the gain Ψ(t )-as defined in (33) must be invertible. This is guaranteed if the following inequality holds true: It can be shown that the CA structure Φ(t ), from (7), satisfies Similarly, using (43), (39) and the fact that ‖B 2 ‖ = 1, it can be verified that Rearranging (44), it can be shown that is a sufficient condition for (42) to be satisfied. Because in gen- [35], the following bound on ‖Ψ(t ) −1 ‖ can be calculated using (45) and the definition of Ψ(t ) from (33): From the definition of the non-linear feedback gain Ω(t ) in (36) and the bound on ‖Ψ(t ) −1 ‖ from (46) Finally, from the small gain theorem [36], if the inequality is satisfied, then the closed loop system described by (37) and (36) is stable. Using the bound on Φ(t ) from (47) and (40), it follows that and therefore, the condition that is sufficient to guarantee stability. Rearranging (50) and substituting into (43) yields the inequality in (41). □

Sliding mode control laws
In this paper, the SMC law for the "virtual" control has been selected to have the general form where the linear component is given by and the non-linear component is given by

Theorem 3. The reachability condition is satisfied if the function (x, t ) in (53) is chosen to satisfy
where is a positive design scalar and Φ max is defined in (43). This ensures that a sliding motion is guaranteed to occur in a finite time and to continue to be sustained for all subsequent time.
Proof. The equation that governs the dynamics of the switching function can be obtained by substituting the 'virtual' control law (51) for v(t ) in the bottom partition of (31). This yieldṡ Using the component-wise definitions of v(t ) in (52) and (53), Equation (55) can be simplified tȯ From the definition of Ψ(t ) in (33) and using the fact that ‖B 2 ‖ = 1, it can be seen that Using the definitions of 2 from (39) andΔ max from (41), Equation (57) can be further simplified to finally through substituting (54) into (58), it can be verified that this choice of (t, x) satisfies the reachability condition [34] s Ṫs ≤ − ‖s‖.
This guarantees s is forced to zero in a finite time and remains zero for all subsequent time i.e there exists a time t = t s at which s(t s ) = 0, and for all time t > t s , s(t ) = 0. □

NON-LINEAR BLENDED WING BODY AIRCRAFT MODEL
This section outlines the development of a non-linear 6 Degree of Freedom (6DoF) BWB model. In [37], several BWB configurations were presented and analysed, which formed the basis of the work in [38]: however, both of these endeavours consider the case of static aerodynamic coefficients. In reality, these quantities are dependent on the aircraft's state. To ensure that the BWB's dynamics are properly captured, the aircraft geometry presented in [37,38] has been reproduced and analysed in TORNADO-a Vortex Lattice Method (VLM) [39]. The results from the analysis have been used to generate lookup tables for the aerodynamic coefficients.

Aircraft geometry
In Table 1, details are provided of seven distinct aerofoil sections. The root aerofoil is represented by section 1 and the winglets are represented by section R. Interpolating between these sections, along the span, and mirroring at the root provides a complete description of the aircrafts' geometry. For the centre-wing (wing section 1), a reflex camber aerofoil (NACA-25116) is chosen to provide a positive (nose-up) pitching moment, aiding trimming of the aircraft. The outer-wing (wing sections 4-6) features supercritical aerofoils (SC2-0610 and SC2-0612) to provide better high-speed performance. The inner and outer wings are blended together through a symmetrical aerofoil (NACA-0016) in the centre (wing sections 2 and 3). The twist distribution is designed to reduce the angle of attack at the wing-tips, reducing the chances of "tip-stall," and at the centre-wing, to better shape the span-wise lift distribution. The aircraft's positive dihedral helps to provide some stability in the roll axis. Additional geometric and inertial properties  Table 2. Note that the Centre of Gravity (CoG), Centre of Pressure (CoP) and the aircraft's mass are considered constant in this paper (although in reality these will be dependent on the aircraft's payload and the flying conditions).

Equations of motion
Two different reference frames are commonly used to describe the motion of an aircraft. The earth reference frame (X e , Y e , Z e ) is used to describe the position of an aircraft's CoG with respect to a fixed point on Earth, whilst the body reference frame (X b , Y b , Z b ) describes the orientation and motion of the aircraft with respect to its CoG. Within the body reference frame, the following standard equations of motion [40,41] can be derived: where u b , v b and w b represent the aircraft's velocities (ms −1 ) along the X b , Y b and Z b axes, respectively. The Euler angles , and represent the aircraft's rotational position (rad ) around the X b , Y b and Z b axes, respectively 1 -these are widely known as the roll, pitch and yaw angles. The rate of change of these angles (rads −1 ) are denoted by p, q and r. The aircraft's mass is denoted bym, the moments of inertia are denoted by J ii , where ii is a pairing of axis (in the body reference frame), and g represents the acceleration due to gravity (ms −2 ). The body forces (F X , F Y , F Z ) (N ) and moments (M l , M m , M n ) (Nm −1 ) are related to the aerodynamic coefficients C i (x, u) through where F eng and M eng represent the forces and moments produced by the aircraft's engines. The matrix DCM Z e eb is a partial directional cosine matrix which transforms a force in the Z e direction (in this case the aircraft's weight) into the body reference frame: this is given by (63) 1 In this paper, rotational angles will follow the "right hand rule" and therefore positive will be defined as a clockwise rotation around an axis when viewed from the origin.  Table 1 labeled) The geometric parametersb andc represent the aircraft's wingspan and the mean aerodynamic chord, respectively, while the dynamic pressure (Pa) is denoted byq and can be calculated throughq wherērepresents air density (Kgm −3 ) and represents the true airspeed of the aircraft (ms −1 ).

Actuator suite
The BWB is equipped with 10 elevons spanning the trailing edge of the wing. These are aggregated into six control surfaces ( f 1 , … , f 6 ) (rad ). The aircraft also features two winglet rudders (r l , r r )(rad ) and three engines-one centered and the other two symmetrically offset (T l , T c , T r ) (N ). The actuator suite is shown in Figure 1. The sign convention used here describes a positive deflection of the elevons as upward (negative in the Z b axis), whilst for the rudders a positive deflection is defined as leftward (negative in the Y b axis). The BWB's three engines are positioned relative to the aircraft's CoG. Each engine is 0.5 m higher than the CoG and inclined upward by 2.5 • . The outboard engines are displaced in the Y b axis by 5 m and canted inward by 2 • . The forces and moments produced by the engines on the aircraft, as used in Equations (61) and (62), can be approximated as: and where T l , T c and T r represent the thrust (N) produced by the engines as labeled in Figure 1.

Aerodynamic coefficients
The aerodynamic coefficients C i , used in Equations (61) and (62), can be described through the equations where is the Angle of Attack (AoA) (rad ), is the sideslip angle (rad ) and V t denotes the true airspeed of the aircraft (ms −1 ). The coefficients in (68) are actually functions of the system state such that C i j = C i j (x); this notation has been simplified in (68) to improve readability. The coefficients C i denote the changes of C i caused by the aerodynamics control surfaces. The total effect of the control surfaces on the aircraft's aerodynamic derivatives (68) can be expressed as where C i m represents the derivatives of the coefficient C i with respect to the m th control surfaces' perturbation from the neutral point, denoted by m. The BWB aircraft, as described in §6.1, was rendered in the VLM program TORNADO [39] to study its aerodynamic properties. From experimentation, it was found that the most significant changes in the aerodynamic coefficients were due to changes in V t and (since m and CoG are assumed constant); therefore, it was decided to base the lookup tables on these parameters. Using a grid of 24 different values of V t (ranging from 50ms −1 to 280ms −1 ) and 61 different values of (ranging from −5 • to 10 • ), a total of 1464 different operating conditions were analysed using the VLM. Results from the analysis were then transformed into lookup tables to approximate the functions C i j (x) as used in (68). The data from some of these lookup tables are shown in Figure 2, in which it can be seen the aerodynamic coefficients vary with the aircrafts state (as opposed to remaining constant). Out of the 48 total lookup tables, those shown in Figure 2 were considered the best cross-section of the entire set, both in terms of their individual shapes and the type of coefficients they represent.

Linear model
For the purposes of control design, the non-linear model described in the previous sections was linearised using the MATLAB functions: linmod and operspec. The linearisation point was chosen at a typical cruise condition with an airspeed of 200ms −1 at an altitude of 3000m. The linearised model showed very little cross-coupling between the lateral and longitudinal dynamics: this allows the problems of lateral and longitudinal control to be addressed independently. In this paper, only lateral control will be considered. The lateral dynamics can be described in the form of (1) where the A and B matrices are shown in (70). The (controlled) output matrix is given by The lateral system's states x(t ), controlled outputs y(t ), and inputs u(t ) are given by where the elements of u(t ) represent perturbations from the trim conditions (not the neutral point) of the corresponding controls in Figure 1. The trim value of u(t ), generated through

CONTROL DESIGN
In this section, the specifics of a controller for the lateral control of a non-linear BWB are discussed. The controller will be based upon applying the theory discussed in Sections 3 and 5 to the linear BWB model from Section 6.5.

Control allocation design
The CA will be based on the distribution structure in (7) where the individual components are defined in (5) and (8). For the nominal component of the CA the associated priority matrix is selected as meaning that only the rudders and the two central elevons (r l , f 3 , f 4 and r r ) are used under nominal conditions, leaving the rest of the controls to focus on longitudinal control. The priority weighting associated with the adaptive component is selected as which prioritises the use of actuators with the greatest range of motion and allows the adaptive component to access the entire actuator suite in the event of the nominal set saturating.

On-line optimisation
The optimisation scheme proposed in Section 4.2 is sufficiently computationally light to be executed in real time. In the following simulations, the optimisation scheme is run every 0.01 s for a total of 100 iterations. The continuous time controller is thus essentially implemented in discrete time at a sample rate of 0.01 s. A "forgetting factor" is implemented on Θ(t ) which allows the adaptive parameter to return to zero after a period of saturation has ended: thus returning to the nominal CA law.
To demonstrate the efficiency of this algorithm, the number of floating point operations required to run it was counted using [42]. Consider the case where m = 10 and l = 2 (i.e. the BWB example studied in this paper). Since at each individual controller update time-step Λ is fixed, it therefore only needs to be calculated once per time-step. This requires a total of approximately 13,000 floating point operations. Calculating the value of Θ n+1 , through (15), requires 1400 operations per iteration. Assuming that the optimisation process is run every 0.01 s for 100 iterations, a total of 15 × 10 6 operations are required each second-this is representative of the simulations presented in this paper. An Intel(R) Core(TM) i7 -9700K processor running at 3.6 GHz is capable of 44GFLOPs (10 9 floating point operations per second), and therefore in this particular example, the optimisation uses less than 0.04% of the processor's capacity. The algorithm is also suitable to run on smaller processors. The Raspberry Pi is capable of 0.132GFLOPs [43] and the ARM-core Cortex-A9 (a common mobile phone processor) is capable of 0.372GFLOPs [44]; therefore, in this example, the optimisation algorithm uses 11.9% and 4.2% of the respective processors' capacity. Note that although this optimisation algorithm is implementable in real time, it still requires a large amount of computation compared to other calculations undertaken to update the control signal: for example, calculating u(t ) for a given value of Θ(t ) and x(t ), through (51) and (7), requires only 334 floating point operations per time-step (3.34 × 10 −6 GFLOPs).

A practical control law with tracking
Here, reference tracking capabilities are introduced by augmenting system (2) with integral action states [34]. This guarantees perfect tracking of fixed reference signals at steady state.
The augmented system state is given byx(t ) = col (x r (t ), x(t )) whereẋ In (76), r (t ) is a (continuous) external reference signal [34] and y(t ) = Cx(t ) represents the controlled outputs, as defined in (1). The overall augmented system iṡx where the system matrices arē If the control signal u(t ) is chosen to stabilise (77), then perfect tracking of a fixed reference signal is achieved [34]. It can be shown by direct evaluation that if (A, B, C ) has no transmission zeros at the origin and (A, B) is controllable, the augmented system (Ā,B) is controllable. A small "feed-forward" modification to v l (t ) in (52) is required to include the reference signal: specifically To avoid the phenomenon of "chattering" [34], the following continuous approximation of (53) is used in the simulations: where is a small positive design scalar. In this paper, = 0.0025 and a fixed value of is chosen where = 0.2. These values were selected through intensive simulation testing. Firstly, is selected large enough to satisfy (54), and therefore ensure sliding occurs (but not so unnecessarily large as to cause excessive "chattering"). Finally, is chosen to be as small as possible to prevent "chattering" whilst maintaining a reasonable approximation of an "ideal" sliding motion.

Sliding surface stability analysis
The sliding surface was designed using the LMI approach outlined in [45], which aims to solve an LQR-like problem whilst minimising the norm 3 . The weighting matrix used for this design is given as where the first two entries (those associated with the integral action states) are weighted heavily to provide good tracking performance.
Using genetic algorithms, suitable values of 0 and 1 can be found as 3.015 and 1.705, respectively. By direct evaluation, the values of 2 and 3 can be calculated as 0.0149 and 0.0176, respectively. Substituting these values into the equality (41) defines the region of guaranteed stability, as shown in Figure 3. The figure shows that the system is most robust to uncertainty when Θ(t ) = 0, at this point stability is guaranteed as long as ‖Δ(t )‖ < 0.29. As ‖Θ(t )‖ increases, the allowable value of ‖Δ(t )‖ decreases exponentially and reaches the value of zero when ‖Θ(t )‖ = 6.01.

SIMULATION RESULTS
The following simulation results demonstrate the effectiveness of the proposed control scheme. The simulations consist of the non-linear BWB (from §6) undertaking a "wing rocking" manoeuvre in three different scenarios: nominal (with perfect actuator health and Θ(t ) = 0), faulty without the adaptive CA element (Θ(t ) = 0) and in a faulty scenario with the adaptive CA. Three different cases will be considered in this section to demonstrate the efficacy of the proposed scheme in maintaining system performance when subjected to actuator faults/failures (especially when dealing with saturation): specifically this paper considers with the fixed, non-optimised CA from [3]-this equivalent to (7) when Θ(t ) = 0 throughout the simulation.
Here, Θ(t ) is not fixed and is chosen on-line using the gradient descent algorithm.

Nominal case
The results for the nominal case are achieved by setting W f (t ) = W f (t ) = I and the adaptive parameter Θ(t ) in (7) is fixed at zero. In Figures 4 and 5, nominal performance is shown. Note that with perfect actuator health, the nominal CA in Figure 5 does not violate any of the actuator's saturation limits (dashed black lines).  This includes complete failures of r l and f 4 which are part of the subset used by the nominal CA in (5). Setting w f i = 0 (from the start of the simulation) essentially simulates a total loss effectiveness whereby the i th actuator is non-responsive to any command signal and remains locked at the trim point, whilst setting 0 < w f i < 1 simulates a "loss of effectiveness" fault. To simulate the uncertainty introduced by an FDI scheme, the diagonal components of Δ(t ) are set as sine waves with: random magnitudes (between 0 and 0.1), random frequencies (between 3rads −1 and 0rads −1 ) and random phase angles. The diagonal components ofŴ f (t ) are bounded such that 0 ≤ŵ f i (t ) ≤ 1.
As in the nominal case Θ(t ) is fixed at zero, effectively reducing the controller to the previous SMC scheme with fixed CA [3,15]. Once faults and failures occur in the nominal actuator set, the performance deteriorates and the manoeuvre induces saturation. This can be seen in Figure 6 and is most clear in the roll angle ( ) response shown in Figure 6a: for example, in response to the first negative reference signal (at 30 s), the rise time is increased from 4.28 to 4.71 s and there is a 66% overshoot (compared to 0% in the nominal case). There is also some clear degradation in the yaw ( ) response shown in Figure 6b. There is a 45% increase in the overshoot between 32 − 34 s and 57 − 59 s; there is also a 2 • increase in overshoot between 38 − 44 s and 62 − 70 s (compared to 0.1 • in the nominal case, giving a 2000% increase). The reason for this can be seen in Figure 7 where the remaining healthy elevon ( f 3 shown in Figure 7b) saturates in an attempt to compensate for the failed actuator ( f 4 shown in Figure 7b). As f 3 enters saturation, the controller attempts to compensate again, by increasing the signal to f 3 , pushing the system deeper into saturation.

Faulty case-optimised CA
In this scenario, the same fault/failure case is considered as with the faulty fixed CA scenario; however, now Θ(t ) is optimised on-line. The system response in Figure 8 demonstrates  the effectiveness of the optimised CA scheme, providing a system response that is identical to the nominal case-despite the actuator faults and failures. Figures 9-11 show the control signals generated by the optimised scheme. It can be seen that, whilst the majority of the control effort is still provided by the remaining actuators within the nominal set ( f 3 and r r ), actuators from outside of the nominal set are used to help prevent saturation. In Figures 9d and 10d, it can be seen that the respective actuators both reach their saturation limits at 30 and 55 s, at which point other actuators, which do not belong to the nominal set, begin to provide some of the control effort, namely: f 1 in Figure 9b and the two engines T l and T r in Figure 11.  In Figure 12, the product of the on-line optimisation is presented and shows that when the nominal CA begins to saturate, the adaptive component weights the saturating actuators negatively, effectively passing on the control effort to the other remaining actuators. Once the period of saturation ends, the effects of the "forgetting factor" can be seen as the weights begin to return to zero, retaining the nominal CA.

Discussion
In Figure 13a, a comparison of the tracking errors for the three different scenarios are shown. The RMS (root mean squared) value of the tracking error signals are shown in Table 3. These values show that, when subjected to faults/failures, the fixed CA has nearly 99% worse tracking performance in and 61% worse tracking performance in . In comparison, the tracking performance of the optimised CA, when subjected to the same set of faults/failures, is basically identical to the fault-free case. Note that the fixed CA controller (developed previously in [3]) has the same "virtual" control law as the proposed scheme, and therefore provides a fair comparison between the new scheme and one that already exists in the literature. The results from the three different scenarios (Figures 4-13) show the shortcomings of the fixed CA controller and the advantages of the proposed scheme. This can most clearly be seen in Table 3. Although the fixed CA has the ability to deal with faults/failures, the controller has no method of accommodating saturation limits. The results clearly demonstrate that the more severe the fault/failure situation, the more likely actuator saturation is to occur due to the extra control effort assigned to the remaining healthy actuators. In the fixed CA case, this results in performance loss. The proposed scheme not only has the ability to handle faults/failures (through the same mechanism as the fixed CA) but also features an efficient adaptive component which allows it to adapt to saturation. In the case of severe faults/failures, this scheme was able to maintain the nominal fault-free performance by preventing saturation.

CONCLUSION
This paper has presented a novel control allocation method which exclusively uses a subset of the actuator suite, until saturation occurs. A sliding mode controller is designed for use with control allocation, and rigorous conditions are derived to guarantee the stability of the closed-loop system. The scheme is shown to be capable of dividing the actuator suite (of an over-actuated system) such that a specified set of actuators is exclusively used; in the event that this set begins to saturate, the entire actuator suite is made available. An on-line optimisation is proposed to adapt the control allocation to ensure the total required control effort is met, whilst keeping actuators within their specified limits. The scheme is proven to be computationally "light" enough to run in real time. The effectiveness of the scheme is demonstrated using a high-fidelity blended wing body aircraft model. Simulation results show the scheme's ability to maintain the nominal fault-free performance despite being presented with severe, uncertain faults/failures and saturation limits. In comparison, when an existing fixed control allocation law is used, some actuators are prone to saturate which results in a substantial amount of performance loss.