2.1 Nonlinear model of aircraft engines

Aircraft engines are characterised by their high complexities and nonlinearities, and can be described by

(1a) \begin{align}\dot{\boldsymbol{x}} & = \boldsymbol{f}(\boldsymbol{x},\boldsymbol{u})\end{align}

(1b) \begin{align}\boldsymbol{y} & = \boldsymbol{g}(\boldsymbol{x},\boldsymbol{u})\end{align}

where ${\boldsymbol{x}=[x_1, \cdots,x_n]^{\textrm{T}} \in {\mathcal R}^{ n}}$ represents state vector, ${\boldsymbol{u} =[u_1, \cdots,u_m]^{\textrm{T}}}\in {\mathcal R}^{ m}$ represents the input vector, $ {\boldsymbol{y} =[y_1, \cdots,y_r]^{\textrm{T}}\in {\mathcal R}^{ r}}$ represents the output vector. $n, m$ , and $r$ are the dimensions of state, input and output, respectively. $\boldsymbol{f}({\cdot})$ and $\boldsymbol{g}({\cdot})$ are nonlinear vector functions, respectively.

For a twin-shaft aircraft engine, state vector $x$ can include low-pressure spool speed $n_l$ and high-pressure spool speed $n_h$ ; control vector $u$ may contain main fuel flow $W_{fb}$ , nozzle throat area $A_8$ ; output vector $y$ can be composed of low-pressure spool speed $n_l$ , total temperature of low-pressure turbine $T_5$ .

A typical general structure of twin-shaft turbofan aircraft engine is shown in Fig. 1 [Reference Richter1]. Fan, low-pressure compressor (LPC) and low-pressure turbine (LPT) are connected by a single shaft and thus rotate synchronously. High-pressure compressor (HPC) and high-pressure turbine (HPT) are located adjacent to the combustor and are connected by a separate shaft. Numbers are used to identify engine stations, only part of numbers are shown in Fig. 1. For example, numbers 20–24 are used to designate points between fan inlet and LPC outlet, and numbers 41–48 is used for HPT. Fuel flow (WF) is delivered by a pump, variable bleed valve (VBV) and variable stator vanes (VSV), which are the actuators utilised in this type of aircraft engine. For some aircraft engines, nozzle throat area is adjustable as well.

Figure 1. General structure of aircraft engine [Reference Richter1].

Hence, for a twin-shaft turbofan aircraft engine, low-pressure spool speed $n_l$ and high-pressure spool speed $n_h$ can be chosen as state vector ${\boldsymbol{x}}$ ; fuel flow $W_{fb}$ and nozzle throat area $A_8$ may be contained in control vector ${\boldsymbol{u}}$ ; some measurable variables, such as high-pressure spool speed $n_h$ , total temperature of low-pressure turbine $T_5$ , can be used to compose output vector ${\boldsymbol{y}}$ , while some calculated variables, such as thrust, compressor surge, can also be combined into output vector ${\boldsymbol{y}}$ [Reference Xiao, Du, Hu and Jiang20, Reference Xiao, Sattarov, Liu and Lin21, Reference Yu, Liu, Bao and Xu29].

Obviously, the dynamic model of aircraft engine in (1) is too complex, which should be avoided in designing control laws to minimise computation load and meet real-time requirements more efficiently. Therefore, linearised aircraft engine dynamic model is often utilised in controller design, particularly for prolonged operation at specific points like the rated working point or cruise point. Generally, the linearised dynamic model of (1) is given as

(2a) \begin{align}\dot{\boldsymbol{x}} & = A\boldsymbol{x} + B\boldsymbol{u}\end{align}

(2b) \begin{align}\boldsymbol{y} & = C\boldsymbol{x} + D\boldsymbol{u}\end{align}

In general, such a linearised model (2) is valid at steady-state operating points, primarily. To enhance the feasibility of model (2) across a wider range, one approach is to incorporate nonlinear terms [Reference Wang, Ye and Hu30], e.g.,

(3a) \begin{align}\dot{\boldsymbol{x}} & = A\boldsymbol{x} + B\boldsymbol{u} + \boldsymbol{G}(\boldsymbol{x})\end{align}

(3b) \begin{align}\boldsymbol{y} & = C\boldsymbol{x} + D\boldsymbol{u}\end{align}

where ${\boldsymbol{G}}({\cdot})$ can be a very general nonlinear term.

Through this extension, aircraft engine model becomes nonlinear again. It is hoped that nonlinear control design technique can be more convenient to carry out, compared with designing controller based on (1), while the control performance can be improved even over a large flight envelope, compared with designing controller based on (2).

Then, let ${\boldsymbol{f}}({\boldsymbol{x}}) \triangleq A{\boldsymbol{x}}+G({\boldsymbol{x}})$ , ${\boldsymbol{h}}({\boldsymbol{x}},{\boldsymbol{u}})\triangleq C{\boldsymbol{x}}+D{\boldsymbol{u}}$ , then system (3) can be written as

(4a) \begin{align}\dot{\boldsymbol{x}} & = \boldsymbol{f}(\boldsymbol{x}) + B\boldsymbol{u}\end{align}

(4b) \begin{align}{\boldsymbol{y}} &= \boldsymbol{h}(\boldsymbol{x},\boldsymbol{u})\end{align}

2.2 Uncertain nonlinear cascaded model of aircraft engines with actuator dynamics and faults

Taking into account manufacturing tolerances, aging of aircraft engine, flight destabilisations, unmodeled nonlinear dynamics, and other factors, uncertainties must be accounted for in the dynamic model of aircraft engine. Thus, when considering uncertainties, corresponding uncertain model of system (4) is

(5a) \begin{align}\dot{\boldsymbol{x}}& = \boldsymbol{f}(\boldsymbol{x}) + B\boldsymbol{u} + \boldsymbol\eta \end{align}

(5b) \begin{align}{\boldsymbol{y}}& = \boldsymbol{h}(\boldsymbol{x},\boldsymbol{u})\end{align}

where ${\boldsymbol \eta} $ is the lumped uncertainty vector, including unmodelled nonlinear dynamics, parameter uncertainties and so on. For example, a description of ${\boldsymbol \eta} $ can be ${\boldsymbol \eta} = \Delta{\boldsymbol{f}}({\boldsymbol{x}})+\Delta B{\boldsymbol{u}}+{\boldsymbol \varrho}(t)$ , where $\Delta{\boldsymbol{f}}({\boldsymbol{x}})$ is the uncertainty of ${\boldsymbol{f}}({\boldsymbol{x}})$ , $\Delta B{\boldsymbol{u}}$ is the uncertainty of $B{\boldsymbol{u}}$ , ${\boldsymbol \varrho}(t)$ is external disruption with respect to time.

In order to analyse the dynamic performance of actuators, first-order inertial lags are used for actuators in the design process [Reference Xiao and Ye31], that is

(6) \begin{equation}\dot{\boldsymbol{u}} = - {\boldsymbol\alpha ^{ - 1}}\boldsymbol{u} + {\boldsymbol\alpha ^{ - 1}}{\boldsymbol{u}_0}\end{equation}

where ${\boldsymbol{u}_0}$ is the control signal before first-order inertial lags, and ${\boldsymbol \alpha}={\textbf {diag}}\{\alpha_1, \alpha_2, ...\alpha_m\}$ is diagonal matrix. $\alpha_i \gt 0, (i=1,\cdots, m)$ are time constants.

Actuators’ faults should be considered, based on (5) and (6), gives

(7a) \begin{align}\dot{\boldsymbol{x}} & = \boldsymbol{f}(\boldsymbol{x}) + B\boldsymbol{u} + {{\boldsymbol\eta }}\end{align}

(7b) \begin{align}\dot{\boldsymbol{u}} & = - {\boldsymbol\alpha ^{ - 1}}\boldsymbol{u} + {\boldsymbol\alpha ^{ - 1}}({{\boldsymbol{u}_0} + {\boldsymbol{u}_{fault}}})\end{align}

(7c) \begin{align}\boldsymbol{y} & = \boldsymbol{h}(\boldsymbol{x},\boldsymbol{u})\end{align}

where ${\boldsymbol{u}_{fault}}$ is actuator fault vector.

Let $\bar{\boldsymbol{x}} = \left[ {\bar{\boldsymbol{x}}_1^{\textrm{T}},\bar{\boldsymbol{x}}_2^{\textrm{T}}{]^{\textrm{T}}} \buildrel \Delta \over = } \right[{{\boldsymbol{x}}^{\textrm{T}}},{{\boldsymbol{u}}^{\textrm{T}}}{]^{\textrm{T}}}$ , ${{\boldsymbol{f}}_1}({{{\bar{\boldsymbol{x}}}_1}}) \buildrel \Delta \over = {\boldsymbol{f}}({\boldsymbol{x}})$ , ${{\boldsymbol{g}}_1}({{{\bar{\boldsymbol{x}}}_1}}) \buildrel \Delta \over = B$ , ${{\boldsymbol{f}}_2}(\bar{\boldsymbol{x}}) \buildrel \Delta \over = - {\boldsymbol\alpha ^{ - 1}}{\boldsymbol{u}}$ , ${{\boldsymbol{g}}_2}(\bar{{\boldsymbol{x}}}) \buildrel \Delta \over = {\boldsymbol\alpha ^{ - 1}}$ , then (7) can be re-written to

(8a) \begin{align}\dot{\bar{\boldsymbol{x}}}_1 & {}={\boldsymbol{f}}_1({\bar{\boldsymbol{x}}}_1)+{\boldsymbol{g}}_1({\bar{\boldsymbol{x}}}_1){\bar{\boldsymbol{x}}}_2+{\boldsymbol \eta} \end{align}

(8b) \begin{align}\dot{\bar{\boldsymbol{x}}}_2 & {}={\boldsymbol{f}}_2({\bar{\boldsymbol{x}}})+{\boldsymbol{g}}_2({\bar{\boldsymbol{x}}})({{\boldsymbol{u}}}_0+{\boldsymbol{u}}_{fault}) \end{align}

(8c) \begin{align}{\boldsymbol{y}} & {}={\boldsymbol{h}}({\bar{\boldsymbol{x}}}) \end{align}

Undoubtedly, (8) shows that aircraft engine model with actuator dynamics is a species of nonlinear cascade systems [Reference Krstic, Kokotovic and Kanellakopoulos32, Reference Wang, Liu, Liu, Li and Wang33], arising mismatched uncertainties and actuator faults [Reference Jiang, Staroswiecki and Cocquempot34, Reference Zhang, Jiang and Cocquempot35]. Its structure is shown in Fig.2.

State ${\bar{\boldsymbol{x}}_2}$ in subsystem (8a) is supposed virtual control [Reference Tang, Dong and Gao36], which can be devised based on the nominal form of system (8).

The nominal system of (8) can be written as

(9a) \begin{align}\dot{\bar{\boldsymbol{x}}}_1 & {}={\boldsymbol{f}}_1({\bar{\boldsymbol{x}}}_1)+{\boldsymbol{g}}_1({\bar{\boldsymbol{x}}}_1){\bar{\boldsymbol{x}}}_2 \end{align}

(9b) \begin{align}\dot{\bar{\boldsymbol{x}}}_2 & {}={\boldsymbol{f}}_2({\bar{\boldsymbol{x}}})+{\boldsymbol{g}}_2({\bar{\boldsymbol{x}}}){{\boldsymbol{u}}}_0\end{align}

(9c) \begin{align}{\boldsymbol{y}} {} & ={\boldsymbol{h}}({\bar{\boldsymbol{x}})}\end{align}

2.3 Control objective

The primary purpose of this paper is to design a novel sliding surface and the corresponding fault-tolerant control method. It will make the following quadratic performance indicator (10) to be minimised when the closed-loop system is approaching the designed sliding surface.

(10) \begin{equation}J=\int^\infty_t{\mathbb{Q}({\bar{\boldsymbol{x}}}_1)+\big( {\bar{\boldsymbol{x}}}_1^{\textrm{T}}Q{\bar{\boldsymbol{x}}}_1+{\bar{\boldsymbol{x}}}_2^{\textrm{T}}R{\bar{\boldsymbol{x}}}_2\big)} d\tau,\end{equation}

where $\mathbb{Q}({\cdot})$ is positive semi-definite matrix, which is associated with the lumped uncertainty, $Q,R$ are weight matrices with appropriate dimensions. Obviously, the performance index (10) reflects the requirement on uncertainty, regulation and control law, simultaneously.

Figure 2. Cascaded structure of aircraft engine control system with actuator dynamics.

Figure 3. Engine line and operating limits on compressor map [Reference Richter1].

3.1 ADP-based sliding surface design

For the system (8), a ADP-based sliding mode function is designed as

(11) \begin{equation}{\boldsymbol{s}} = {\bar{\boldsymbol{x}}_2} + \frac{1}{2}{R^{ - 1}}{\boldsymbol{g}}_1^T({{{\bar{\boldsymbol{x}}}_1}}){(\nabla {\boldsymbol\sigma _c}( {\bar{{\boldsymbol{x}}}}_1))^{\textrm{T}}}{\hat{\boldsymbol{W}}_c}\end{equation}

where ${\boldsymbol{s}} = {[{s_1}, \cdots, {s_m}]^{\textrm{T}}} \in {\mathcal{R}^m}$ is sliding mode variable for system (8), ${\boldsymbol\sigma _c}({\cdot})$ is an activation function vector of critic neural network (critic NN) in ADP, $\nabla {\boldsymbol\sigma _c}({\cdot})$ is gradient with respect to ${\bar{\boldsymbol{x}}_1}$ denoted as $\nabla {\boldsymbol\sigma _c} \buildrel \Delta \over = \frac{{\partial {\boldsymbol\sigma _c}}}{{\partial {{{\bar{{\boldsymbol{x}}}}}_1}}}$ , and ${\hat{\boldsymbol{W}}_c}$ is an estimation of ideal weight vector in ADP. The corresponding sliding surface is $\mathcal{S} = \{ \bar{\boldsymbol{x}}|{\boldsymbol{s}}( {\bar{\boldsymbol{x}}} ) = \boldsymbol 0\} $ .

3.2 Sliding mode fault-tolerant controller design

Let $\widehat {\bar{\boldsymbol{x}}}_2^{{*}} \buildrel \Delta \over = - \frac{1}{2}{R^{ - 1}}{\boldsymbol{g}}_1^T({{{\bar{\boldsymbol{x}}}_1}}){(\nabla {\boldsymbol\sigma _c}({{{\bar{\boldsymbol{x}}}_1}}))^T}{\hat{\boldsymbol{W}}_c}$ , then (11) turns to

(12) \begin{equation}s = {\bar{\boldsymbol{x}}_2} - \hat {\bar{\boldsymbol{x}}}_2^{{*}}\end{equation}

hence, the derivative of sliding mode variable, i.e. $\dot{\boldsymbol{s}}$ , can be written as (13) according to (11) and (8b).

(13) \begin{align}\dot{\boldsymbol{s}}& {}=\dot{\bar {\boldsymbol{x}}}_2-\dot{\hat{ {\bar {\boldsymbol{x}}}}}_2^*\nonumber\\& {}={\boldsymbol{f}}_2( {\bar {\boldsymbol{x}}})+{\boldsymbol{g}}_2( {\bar {\boldsymbol{x}}})({\boldsymbol{u}}_0+{\boldsymbol{u}}_{fault})-\dot{\hat{\bar{\boldsymbol{x}}}}_2^*\end{align}

The novel sliding mode fault-tolerant controller ${\boldsymbol{u}}_0$ with ADP-based sliding surface $\mathcal{S}$ is constructed as

(14) \begin{equation}{\boldsymbol{u}}_0={\boldsymbol{g}}_2^{-1}( {\bar {\boldsymbol{x}}})\bigg(\dot{\hat{\boldsymbol{x}}}_2^*-{\boldsymbol{f}}_2({\boldsymbol{x}})-{\boldsymbol{K}}_s \big(\textbf{sgn}({\boldsymbol{s}})\big)\bigg)-\hat{\boldsymbol{u}}_{fault},\end{equation}

with adaptive fault estimation law

(15) \begin{equation}\dot{\hat{\boldsymbol{u}}}_{fault}={\boldsymbol{K}}_f {\boldsymbol{g}}_2^{\textrm{T}}({\bar{\boldsymbol{x}}}){{\boldsymbol{s}}},\end{equation}

where ${{\boldsymbol{K}}_s},{{\boldsymbol{K}}_f}$ are positive control gain matrix (14) and learning rate matrix in adaptive law (15), respectively. ${\mathbf{sgn}}(s) = {[{\mathbf{sgn}}({{s_1}}), \cdots, {\mathbf{sgn}}( {{s_m}} )]^{\textrm{T}}}$ , ${\mathbf{sgn}}({{s_i}}),( {i = 1, \cdots, m} )$ are sign functions. ${\hat{\boldsymbol{u}}_{fault}}$ is estimation of actuator fault vector.

3.3 Robustness analysis

3.3.1 Stability of sliding surface

When analysing stability of sliding surface $\mathcal{S}$ , nominal system (9) is considered. According to (9a), ${\bar{\boldsymbol{x}}_2}$ can be considered as virtual control input of subsystem (10). Considering the performance index (21), optimal cost function (OCF) can be defined as

(16) \begin{equation}J^*({\bar{\boldsymbol{x}}}_1)=\min_{{\bar{\boldsymbol{x}}}_2\in \textsf{A}(\Omega)}\int^\infty_t{\mathbb{Q}({\bar{\boldsymbol{x}}}_1)+\big( {\bar{\boldsymbol{x}}}_1^{\textrm{T}}Q{\bar{\boldsymbol{x}}}_1+{\bar{\boldsymbol{x}}}_2^{\textrm{T}}R{\bar{\boldsymbol{x}}}_2\big)} d\tau,\end{equation}

where ${\textsf{A}}(\Omega)$ is admissible control set on the compact set ${{\Omega }}$ , and the Hamilton-Jacobi-Bellman (HJB) equation of this system is given as

(17) \begin{equation}\begin{array}{*{20}{l}}{0 = \mathop {{\textrm{min}}}\limits_{{{\bar{\boldsymbol{x}}}_2} \in \textsf{A}(\Omega)} H\left( {{{\bar{\boldsymbol{x}}}_1},{{\bar{\boldsymbol{x}}}_2},\nabla {J^{{*}}}({{{\bar{\boldsymbol{x}}}_1}})} \right)}\end{array}\end{equation}

where $H({\cdot})$ is the Hamiltonian function of this optimal control problem

(18) \begin{align}H({{\bar{\boldsymbol{x}}_1},{\bar{\boldsymbol{x}}_2},\nabla J({{\bar{\boldsymbol{x}}_1}})}) = &\ \mathbb{Q}({{\bar{\boldsymbol{x}}_1}}) + ({\bar{\boldsymbol{x}}_1^{\textrm{T}}Q{\bar{\boldsymbol{x}}_1} + \bar{\boldsymbol{x}}_2^{\textrm{T}}R{\bar{\boldsymbol{x}}_2}})\nonumber\\& + {(\nabla J({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}( {{{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) + {{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}}){\bar{\boldsymbol{x}}_2}} ).\end{align}

with $J(0) = 0$ .

Therefore, by differentiating (16) with respect to $\bar{\boldsymbol{x}}_2^{{*}}$ , the optimal control policy is

(19) \begin{equation}\bar{\boldsymbol{x}}_2^{{*}} = - \frac{1}{2}{R^{ - 1}}{\boldsymbol{g}}_1^{\textrm{T}}({{\bar{\boldsymbol{x}}_1}})\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}})\end{equation}

Considering control law (19) is unbounded, a generalised nonquadratic utility function[Reference Abu-Khalaf and Lewis39] can be used to resolve problem of input constraints.

Considering (17), by substituting the optimal control policy (19) into Hamiltonian function (18), HJB Equation (17) can be given as

(20) \begin{align}0 = &\ \mathbb{Q}({{\bar{\boldsymbol{x}}_1}}) + ( {\bar{\boldsymbol{x}}_1^{\textrm{T}}Q{\bar{\boldsymbol{x}}_1} + \bar{\boldsymbol{x}}{{_2^{{*}}}^{\textrm{T}}}R\bar{\boldsymbol{x}}_2^{{*}}} ) + {(\nabla J({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}({{{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) + {{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}})\bar{\boldsymbol{x}}_2^{{*}}})\nonumber\\ = &\ \mathbb{Q}({{\bar{\boldsymbol{x}}_1}}) + \bar{\boldsymbol{x}}_1^{\textrm{T}}Q{\bar{\boldsymbol{x}}_1} + {(\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}{{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}})\nonumber\\ & - \frac{1}{4}{(\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}{{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}}){R^{ - 1}}{\boldsymbol{g}}_1^T({{\bar{\boldsymbol{x}}_1}})\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}})\end{align}

which means that $H( {{\bar{\boldsymbol{x}}_1},\bar{\boldsymbol{x}}_2^{{*}},\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}})} ) = 0$ .

The choice of $\mathbb{Q}({\bar{\boldsymbol{x}}}_1)$ is important in robust stabilisation scheme. Here, $\mathbb{Q}({\bar{\boldsymbol{x}}}_1)$ is specified as

(21) \begin{equation}\begin{array}{*{20}{l}}{\mathbb{Q}({{\bar{\boldsymbol{x}}_1}}) = \frac{1}{8}{{(\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}}))}^{\textrm{T}}}\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}}) + 2{{\bar{\boldsymbol \eta} }^2}}\end{array}\end{equation}

and $\mathbb{Q}({\bar{\boldsymbol{x}}}_1) \geq 0$ meets. Based on the help of this form, the following lemma showing stability of the sliding-mode dynamics of the systems (9) is derived as follows.

Lemma 1. [Reference Du, Jiang, Ma and Cheng25] Consider nominal subsystem (9a) and quadratic performance index (10) with the specified form of $\mathbb{Q}({{\bar{\boldsymbol{x}}_1}})$ (21), it is supposed that the HJB Equation (20) is a solution of ${J^{{*}}}({{\bar{\boldsymbol{x}}_1}})$ . $\bar{\boldsymbol{x}}_2^{{*}}$ is given by (19), and the sliding surface defined by $\boldsymbol{s} = {\bar{\boldsymbol{x}}_2} - \bar{\boldsymbol{x}}_2^{{*}} = 0$ exists. On this sliding surface, sliding-mode dynamics of this mismatched uncertain nonlinear system described by (8) is asymptotic stability.

Proof of Lemma 1. On the basis of the Theorem 1 in [Reference Wang, Liu, Mu and Zhang40], it is evident that the optimal control policy (19) can guarantee asymptotic stability of the uncertain nonlinear subsystem (8a). On the sliding surface, ${\bar{\boldsymbol{x}}_2} = \bar{\boldsymbol{x}}_2^*({{\bar{\boldsymbol{x}}_1}})$ is hold and explicitly state ${\bar{\boldsymbol{x}}_2}$ is asymptotically convergent to the equilibrium point as well.

Because the analytical solution of the HJB Equation (20) is rather difficult to calculate, ADP with NN approximation is used in the following to find the siding surface by obtaining a numerical solution for (20).

Assumption 3. On the compact set $\Omega $ , The ideal NN weights ${\boldsymbol{W}_c}$ are bounded, the NN activation functions and their gradient are bounded, the NN approximation error and its gradient are bounded, i.e., $\parallel {{\boldsymbol{W}}_c}\parallel \le {\lambda _W}$ , $\parallel {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}})\parallel \le {\lambda _\sigma }$ , $\parallel \nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}})\parallel \le {\lambda _{\nabla \sigma }}$ , $\parallel {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}})\parallel \le {\lambda _\varepsilon }$ , and $\parallel \nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}})\parallel \le {\lambda _{\nabla \varepsilon }}$ , where ${\lambda _W},{\lambda _\sigma },{\lambda _{\nabla \sigma }},{\lambda _\varepsilon },{\lambda _{\nabla \varepsilon }}$ are existing unknown boundaries.

Assumption 4. For the nominal subsystem (9a) with the quadratic performance index (10) and substituting the optimal control policy (19) into the subsystem (9a), let ${J_s}({{\bar{\boldsymbol{x}}_1}})$ be a continuously differentiable function and simultaneously a Lyapunov function candidate, it is assumed that ${\dot J_s}({{\bar{\boldsymbol{x}}_1}})$ can be written as

(22) \begin{equation}{\dot J_s}({{\bar{\boldsymbol{x}}_1}}) = {(\nabla {J_s}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}\left( {{{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) + {{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}})\bar{\boldsymbol{x}}_2^{{*}}} \right) \lt 0.\end{equation}

and the positive definite matrix ${{\Gamma }} \in {\mathcal{R}^n}$ ensures that (23) is true, where ${\lambda _{min}}({{\Gamma }})$ represents the minimal eigenvalue of the matrix ${{\Gamma }}$ .

(23) \begin{equation}{(\nabla {J_s}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}\left( {{{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) + {{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}})\bar{\boldsymbol{x}}_2^{{*}}} \right) = - {(\nabla {J_s}({{\bar{\boldsymbol{x}}_1}}))^T}\Gamma \nabla {J_s}({{\bar{\boldsymbol{x}}_1}}) \le - {\lambda _{min}}(\Gamma)\parallel \nabla {J_s}({\bar{\boldsymbol{x}}}_1){\parallel ^2}\end{equation}

Given the sake of universal approximation property of NN, approximate ${J^{{*}}}({{\bar{\boldsymbol{x}}_1}})$ by NN with only one hidden layer can be given as

(24) \begin{equation}{J^{{*}}}({{\bar{\boldsymbol{x}}_1}}) = {\boldsymbol{W}}_c^T{\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}) + {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}})\end{equation}

where ${{\boldsymbol{W}}_c} \in {\mathcal{R}^p}$ is the ideal weight, ${\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}) \in {\mathcal{R}^p}$ is the activation function, $p$ is the number of neurons, and ${\varepsilon _c}({\bar{\boldsymbol{x}}}_1)$ is the unknown approximation error of this NN.

Thus, partial derivative of the OCF (24) is written as

(25) \begin{align}\nabla {J^{{*}}}({{\bar{\boldsymbol{x}}_1}})& = \frac{{\partial {J^{{*}}}({\bar{\boldsymbol{x}}}_1)}}{{\partial {{{\bar{\boldsymbol{x}}}}_1}}}\nonumber\\& = {(\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}))^T}{{\boldsymbol{W}}_c} + \nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}).\end{align}

and the optimal control policy is yielded as

(26) \begin{equation}\bar{\boldsymbol{x}}_2^{{*}} = - \frac{1}{2}{R^{ - 1}}{\boldsymbol{g}}_1^{\textrm{T}}( \bar{\boldsymbol{x}} )(( {\nabla {\boldsymbol\sigma _c}( \bar{\boldsymbol{x}} ){)^{\textrm{T}}}{{\boldsymbol{W}}_c} + \nabla {\varepsilon _c}( \bar{\boldsymbol{x}} )} ).\end{equation}

Define two non-negative matrices $\mathcal{A}({{\bar{\boldsymbol{x}}_1}}),\mathcal{B}({{\bar{\boldsymbol{x}}_1}})$ as

\begin{align*} & \mathcal{A}({{\bar{\boldsymbol{x}}_1}}) \buildrel \Delta \over = \nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}}){R^{ - 1}}{\boldsymbol{g}}_1^{\textrm{T}}({{\bar{\boldsymbol{x}}_1}})({(\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}\\&\mathcal{B}({{\bar{\boldsymbol{x}}_1}}) \buildrel \Delta \over = \nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}})({(\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}.\end{align*}

With the quadratic performance index (10) and control policy (26) denoted by the NN, the Hamiltonian function could be given by

(27) \begin{align}H( {{\bar{\boldsymbol{x}}_1},{{\boldsymbol{W}}_c}} ) = &\ \bar{\boldsymbol{x}}_1^TQ{\bar{\boldsymbol{x}}_1} + {\boldsymbol{W}}_c^{\textrm{T}}\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) - \frac{1}{4}{\boldsymbol{W}}_c^{\textrm{T}}\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c} + 2{\bar{\boldsymbol \eta} ^2}\nonumber\\& + \frac{1}{8}{\boldsymbol{W}}_c^{\textrm{T}}\mathcal{B}({\bar{\boldsymbol{x}}}_1){{\boldsymbol{W}}_c} + {e_{cH}} = 0\end{align}

where the term ${e_{cH}}$ is the residual error of the NN expression.

(28) \begin{align}{e_{cH}} = &\ {(\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}))^T}{{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) - \frac{1}{2}{(\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}))^T}{{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}}){R^{ - 1}}{\boldsymbol{g}}_1^T({{\bar{\boldsymbol{x}}_1}})({(\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}))^T}{{\boldsymbol{W}}_c}\nonumber\\& - \frac{1}{4}{(\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}{{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}}){R^{ - 1}}{\boldsymbol{g}}_1^{\textrm{T}}({{\bar{\boldsymbol{x}}_1}})\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}) + \frac{1}{4}{(\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}{(\nabla {\boldsymbol\sigma _c}({{x_1}}))^{\textrm{T}}}{{\boldsymbol{W}}_c}\nonumber\\& + \frac{1}{8}{(\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}\nabla {\varepsilon _c}({{\bar{\boldsymbol{x}}_1}})\end{align}

In the view that the ideal NN weights are unknown, the estimation of NN weight vector ${\hat {\boldsymbol{W}}_c}$ is utilised to approximate the OCF (24)

(29) \begin{equation}{\hat J^{{*}}}({{\bar{\boldsymbol{x}}_1}}) = \hat {\boldsymbol{W}}_c^{\textrm{T}}{\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}})\end{equation}

then an approximation of the partial derivative of the OCF (24) is

(30) \begin{equation}\nabla {\hat J^{{*}}}({{\bar{\boldsymbol{x}}_1}}) = \frac{{\partial {{\hat J}^{{*}}}({{\bar{\boldsymbol{x}}_1}})}}{{\partial {x_1}}} = {(\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}))^{\textrm{T}}}{\hat {\boldsymbol{W}}_c}.\end{equation}

the approximate optimal control policy can be expressed as

(31) \begin{equation} \hat{\bar{\boldsymbol{x}}}_2^{{*}} = - \frac{1}{2}{R^{ - 1}}{\boldsymbol{g}}_1^{\textrm{T}}({{\bar{\boldsymbol{x}}_1}}){{(\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}))}^{\textrm{T}}}{{\hat{\boldsymbol{W}}}_c}.\end{equation}

and the approximate Hamiltonian function is given as

(32) \begin{align}\hat H( {{\bar{\boldsymbol{x}}_1},{{\boldsymbol{W}}_c}} ) = &\ \bar{\boldsymbol{x}}_1^{\textrm{T}}Q{\bar{\boldsymbol{x}}_1} + \hat {\boldsymbol{W}}_c^{\textrm{T}}\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}})\nonumber\\& - \frac{1}{4}\hat {\boldsymbol{W}}_c^T\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){\hat {\boldsymbol{W}}_c} + 2{\bar{\boldsymbol \eta} ^2}({{\bar{\boldsymbol{x}}_1}}) + \frac{1}{8}\hat {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){\hat {\boldsymbol{W}}_c}\end{align}

Define ${e_c} = \hat H( {{\bar{\boldsymbol{x}}_1},{{\boldsymbol{W}}_c}} ) - H( {{\bar{\boldsymbol{x}}_1},{{\boldsymbol{W}}_c}} )$ , and ${\tilde {\boldsymbol{W}}_c} = {{\boldsymbol{W}}_c} - {\hat {\boldsymbol{W}}_c}$ . According to (27) and (32), the ${e_c}$ with respect to ${\tilde {\boldsymbol{W}}_c}$ is expressed as

(33) \begin{align}{e_c} = &\ \hat H( {{\bar{\boldsymbol{x}}_1},{{\boldsymbol{W}}_c}} ) - H( {{\bar{\boldsymbol{x}}_1},{{\boldsymbol{W}}_c}} )\nonumber\\ = & \ - \tilde {\boldsymbol{W}}_c^T\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) - \frac{1}{4}\tilde {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){\tilde {\boldsymbol{W}}_c} + \frac{1}{2}\tilde {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c}\nonumber\\& + \frac{1}{8}\tilde {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){\tilde {\boldsymbol{W}}_c} - \frac{1}{4}\tilde {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{B}({\bar{\boldsymbol{x}}}_1){{\boldsymbol{W}}_c} - {e_{cH}}.\end{align}

Now, it is required to train the critic NN and to find an updated law of weights by minimising the simple cost standard ${E_c} = ({1/2})e_c^2$ .

An improved weight updated rule [Reference Wang, Liu, Mu and Zhang40] is chosen to avoid the difficulty of obtaining a valid initial control policy.

Then, based on Assumption 4, ${J_s}({\bar{\boldsymbol{x}}}_1)$ is given as ${J_s}({{\bar{\boldsymbol{x}}_1}}) = 0.5\bar{\boldsymbol{x}}_1^{\textrm{T}}{\bar{\boldsymbol{x}}_1}$ , and a weight updated law of the critic NN is designed as

(34) \begin{equation}\dot{\hat{{\boldsymbol{W}}}}_{c}=-\alpha_{c}\bigg(\frac{\partial E_c}{\partial\hat{{\boldsymbol{W}}}_{c}}\bigg)-\alpha_{s}\bigg(\frac{\partial[(\nabla J_{s}({\bar{\boldsymbol{x}}}_1))^{\textrm{T}}({\boldsymbol{f}}_1({\bar{\boldsymbol{x}}}_1)+{\boldsymbol{g}}_1({\bar{\boldsymbol{x}}}_1)\hat{{\bar{\boldsymbol{x}}}}_2^*)]}{\partial \hat{{\boldsymbol{W}}}_{c}}\bigg)\end{equation}

where ${\alpha _c} \gt 0$ and ${\alpha _s} \gt 0$ are the basic update rate of critic NN and the learning rate of the additional stabilisation term, respectively.

By substituting the optimal control strategy (31) into (34), the updated law can be expressed as

(35) \begin{equation}\dot{\hat{{\boldsymbol{W}}}}_{c}=-\alpha_{c}\bigg(\frac{\partial E_c}{\partial\hat{{\boldsymbol{W}}}_{c}}\bigg)+\frac{1}{2}\alpha_{s}\nabla{\boldsymbol \sigma}_{c}({\bar{\boldsymbol{x}}}_1){\boldsymbol{g}}_1({\bar{\boldsymbol{x}}}_1) R^{-1}{\boldsymbol{g}}_1^{\textrm{T}}({\bar{\boldsymbol{x}}}_1)\nabla J_{s}({\bar{\boldsymbol{x}}}_1)\end{equation}

Let the weight estimation error (WEE) ${\tilde {\boldsymbol{W}}_c}$ is ${\tilde {\boldsymbol{W}}_c} = {{\boldsymbol{W}}_c} - {\hat {\boldsymbol{W}}_c}$ , due to ${{\boldsymbol{W}}_c}$ will be constant at last, $\dot{\tilde{{\boldsymbol{W}}}}_{c}=-\dot{\hat{{\boldsymbol{W}}}}_{c}$ and the dynamics of the WEE can be written as

(36) \begin{equation}\dot{\tilde{{\boldsymbol{W}}}}_{c}= {}\alpha_{c}\bigg(\frac{\partial E_c}{\partial\hat{{\boldsymbol{W}}}_{c}}\bigg)-\frac{1}{2}\alpha_{s}\nabla{\boldsymbol \sigma}_{c}({\bar{\boldsymbol{x}}}_1){\boldsymbol{g}}_1({\bar{\boldsymbol{x}}}_1)R^{-1}{\boldsymbol{g}}_1^{\textrm{T}}({\bar{\boldsymbol{x}}}_1)\nabla J_{s}({\bar{\boldsymbol{x}}}_1),\end{equation}

where

\begin{align*}\frac{{\partial {E_c}}}{{\partial {{\boldsymbol{W}}_c}}} & = {}{{e_c}\frac{{\partial {e_c}}}{{\partial {{\boldsymbol{W}}_c}}}}\nonumber\\ & = {}{{e_c}( {\nabla {\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{f}}_1}({{\bar{\boldsymbol{x}}_1}}) - \frac{1}{2}\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c} + \frac{1}{4}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c}}).}\end{align*}

thus

(37) \begin{align}\dot{\tilde{{\boldsymbol{W}}}}_{c}= &\ {}\alpha_{c}\bigg(-\tilde{{\boldsymbol{W}}}_{c}^{\textrm{T}}\nabla{\boldsymbol \sigma}_{c}({\bar{\boldsymbol{x}}}_1){\boldsymbol{f}}_1({\bar{\boldsymbol{x}}}_1) -\frac{1}{4}\tilde{{\boldsymbol{W}}}_{c}^{\textrm{T}}\mathcal{A}({\bar{\boldsymbol{x}}}_1)\tilde{{\boldsymbol{W}}}_{c}\nonumber\\& + \frac{1}{2}\tilde{\boldsymbol{W}}_{\textrm{c}}^{\textrm{T}}\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c} + \frac{1}{8}\tilde {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){\tilde {\boldsymbol{W}}_c} - \frac{1}{4}\tilde {\boldsymbol{W}}_c^{\textrm{T}}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c} - {e_{cH}})\nonumber\\& \times \left( {\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{f}}_1} + \frac{1}{2}\mathcal{A}({{x_1}}){{\tilde{\boldsymbol{W}}}_c} - \frac{1}{2}\mathcal{A}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c} - \frac{1}{4}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){{\tilde{\boldsymbol{W}}}_c} + \frac{1}{4}\mathcal{B}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{W}}_c}} \right)\nonumber\\& - \frac{1}{2}{\alpha _s}\nabla {\boldsymbol\sigma _c}({{\bar{\boldsymbol{x}}_1}}){{\boldsymbol{g}}_1}({{\bar{\boldsymbol{x}}_1}}){R^{ - 1}}{\boldsymbol{g}}_1^{\textrm{T}}({{\bar{\boldsymbol{x}}_1}})\nabla {J_s}({{\bar{\boldsymbol{x}}_1}}).\end{align}

Therefore, the following Lemma 2 illustrates the system states’ property on the novel sliding surface.

Lemma 2. In view of the nominal system described by (9) under the provided ADP-based sliding surface (11), the states of system (8) on the sliding surface and the weight error dynamics are uniformly ultimately bounded (UUB) when the sliding mode exists.

Proof of Lemma 2. The definition of UUB can be found in [Reference Vamvoudakis and Lewis41], while according to the theorem 2 in [Reference Wang, Liu, Mu and Zhang40], one can see that state ${\bar{\boldsymbol{x}}_1}$ on the sliding surface specified by (11) and dynamics of WEE are both UUB. For the sake of brevity, the proof process will not be proved repeatedly. Besides, in term of the sliding surface $\mathcal{S} = \{ \bar{\boldsymbol{x}}|{\boldsymbol{s}}( \bar{\boldsymbol{x}} ) = 0\} $ , state ${{\bar{\mathbf{x}}}_2}$ of the system on this sliding surface is also UUB.

Remark 3 It is easy to find that Assumption 4 is a general assumption, which can see details in [Reference Wang, Liu, Zhang and Li42, Reference Du, Jiang and Ma43]. Besides, though the use of ADP faces the problem of the training taking too much time, it can be avoided by off-line training to make sure the real-time performance of the whole control system.

3.3.2 Reaching ability of sliding surface

Theorem 1. For the system described by (8) with mismatched uncertainty, actuator dynamics and faults, under the proposed control law (14), with the ADP-based sliding surface described by (11), and adaptive fault estimation law (15), the closed-loop system’s states are UUB and its fault estimation error is bounded.

Proof [Proof of Theorem 1] Select Lyapunov function as

(38) \begin{equation}V(t) = \frac{1}{2}{\boldsymbol{s}^{\textrm{T}}}\boldsymbol{s} + \frac{1}{2}\tilde{\boldsymbol{u}}_{{\textrm{fault}}}^{\textrm{T}}\boldsymbol{K}_f^{ - 1}{\tilde{\boldsymbol{u}}_{fault}}\end{equation}

where ${\tilde u_{fault}}$ is fault estimation error, which is defined as ${\tilde{\boldsymbol{u}}_{fault}} = {\boldsymbol{u}_{fault}} - {\hat{\boldsymbol{u}}_{fault}}$ .

One the basis of Assumption 2, along the sliding mode variable ${\mathbf{s}}$ , the derivative of the Lyapunov function (38) is

(39) \begin{align}\dot{V}(t)& {}={\boldsymbol{s}}^{\textrm{T}}\dot{{\boldsymbol{s}}}-\tilde{{\boldsymbol{u}}}_{fault}^{\textrm{T}}{\boldsymbol{K}}_f^{-1}\dot{\hat{{\boldsymbol{u}}}}_{fault}\nonumber\\& {}={\boldsymbol{s}}^{\textrm{T}}\big({\boldsymbol{f}}_2({\bar {\boldsymbol{x}}})+{\boldsymbol{g}}_2({\bar {\boldsymbol{x}}})({\boldsymbol{u}}_0+{\boldsymbol{u}}_{fault})-\dot{\hat{\bar{{\boldsymbol{x}}}}}_2^*({\bar {\boldsymbol{x}}}_1)\big)-\tilde{{\boldsymbol{u}}}_{fault}^{\textrm{T}}{\boldsymbol{K}}_f^{-1}\dot{\hat{{\boldsymbol{u}}}}_{fault}\end{align}

Based on the control law (14) and the actuator fault compensation law (15), (39) will be expressed as

(40) \begin{align}\dot{V}(t)& {}=-{K_s}\|{\boldsymbol{s}}\|+{\boldsymbol{s}}^{\textrm{T}}{\boldsymbol{g}}_2(x)\tilde{{\boldsymbol{u}}}_{fault}-\tilde{{\boldsymbol{u}}}_{fault}^{\textrm{T}}{\boldsymbol{K}}_f^{-1}\dot{\hat{{\boldsymbol{u}}}}_{fault} {}\nonumber\\ & = - {K_s}\parallel s\parallel\nonumber \\ & \le 0.\end{align}

Therefore, it is clear that all states of closed-loop system are bounded. In terms of the Lyapunov theorem and Lemma 2, the sliding modes exist and the overall system exists in UUB states.

Here complete the proof.

3.4 Intelligent parameters optimisation for ADP-based SMFTC

In this section, GWO will be employed to search the optimal control designable parameters in (14) and (15) under the multi-target performance indicators of aircraft engine are

(41a) \begin{align}J_{GWO_1} {} & =\int_0^t{{\boldsymbol{x}}^{\textrm{T}}{\boldsymbol{x}}}d\tau\end{align}

(41b) \begin{align}J_{GWO_2} & {}=\int_0^t{{\boldsymbol{u}}_0^{\textrm{T}}{\boldsymbol{u}}_0}d\tau\end{align}

where $J_{GWO_1}$ reflects steady-state control accuracy requirement, and $J_{GWO_2}$ expresses the interest on input energy consumption, in view that the state deviation from the equilibrium point is expected to be as small as possible, while a minimum input energy consumption is also desired at the same time.

Generally speaking, GWO is a meta-heuristic optimisation algorithm, which is encouraged by the primary phases of grey wolf hunting.

The main mathematical equations of describing the hunting process are:

(42) \begin{align} &\vec X({{t_{GWO}} + 1}) = {{\vec X}_p}({{t_{GWO}}}) - \vec A\vec D \nonumber\\ & {}{{\vec A}} = 2\vec a{{\vec r}_1} - \vec a,{{\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;}}\vec C = 2{{\vec r}_2},{{\;\;\;\;\;\;\;\;}}\vec D = \left| {\vec C{{\vec X}_p}({{t_{GWO}}}) - \vec X({{t_{GWO}}})} \right| \nonumber\\ & {}{{\vec D}_\alpha } = \left| {{{\vec C}_1}{{\vec X}_\alpha } - \vec X} \right|,{{\;\;\;\;\;\;\;\;}}{{\vec D}_\beta } = \left| {{{\vec C}_2}{{\vec X}_\beta } - \vec X} \right|,{{\;\;\;\;\;\;\;\;}}{{\vec D}_\delta } = \left| {{{\vec C}_3}{{\vec X}_\beta } - \vec X} \right|, \\ & {}{{\vec X}_1} = {{\vec X}_\alpha } - {{\vec A}_1}{{\vec D}_\alpha },{{\;\;\;\;\;\;\;\;}}{{\vec X}_2} = {{\vec X}_\beta } - {{\vec A}_2}{{\vec D}_\beta },{{\;\;\;\;\;\;\;\;}}{{\vec X}_3} = {{\vec X}_\delta } - {{\vec A}_3}{{\vec D}_\delta },\nonumber \\ & {}\vec X({{t_{GWO}} + 1}) = {\frac{{{\vec X}_1} + {{\vec X}_2} + {{\vec X}_3}}{3}} \nonumber \end{align}

where the current iteration is represented by ${t_{{{{\;}}_{GWO}}}}$ , the position vector of the prey and a grey wolf are given by ${\vec X_p}( {{t_{GWO}}} )$ and $\vec X({{t_{GWO}}})$ , respectively. The subscripts ${({\cdot})_\alpha }$ , ${({\cdot})_\beta }$ and ${({\cdot})_\delta }$ are in behalf of the ${\boldsymbol\alpha _{{{}_{GWO}}}}$ (best candidate solution), ${\boldsymbol\beta _{{{}_{GWO}}}}$ and ${\boldsymbol\delta _{{{}_{GWO}}}}$ grey wolves, respectively. $\vec A$ and $\vec C$ are coefficient vectors, with $\vec a$ is decreased from 2 to 0 in the course of iterations, linearly. ${\vec r_1}$ and ${\vec r_2}$ are stochastic vectors in the range of $[{0,1}]$ . $\vec D$ , ${\vec X_1}$ , ${\vec X_2}$ and ${\vec X_3}$ are instrumental variables.

The optimisation imitates the hunting process of grey wolves, the detailed principle of GWO can be seen in Refs [Reference Xiao, Sattarov, Liu and Lin21, Reference Mirjalili, Mohammad Mirjalili and Lewis44].

According to the SMC fundamental theory [Reference Gao38], ${\boldsymbol{K}_s}$ is crucial for the control performance of aircraft engine. Meanwhile, the ${\boldsymbol{K}_f}$ is important for the fault estimation, as shown in the control law (14) and adaptive law (15). Thus, ${\boldsymbol{K}_s}$ and ${\boldsymbol{K}_f}$ are key designable parameter matrices in ADP-based SMFTC controller for obtaining the minimum performance indexes (41).

Therefore, when GWO is applied to reach the parameter optimisation of ADP-based SMFTC for aircraft engine, the $\vec X({{t_{GWO}}})$ described by (42) for each ${t_{GWO}}$ can be $\vec X = \left[ {{K_{{s_1}}}, \cdots, {K_{{s_n}}},{K_{{{\boldsymbol{f}}_1}}}, \cdots, {K_{{f_n}}}} \right]$ , where ${K_{{s_i}}},{K_{{f_i}}},( {i = 1, \cdots, n})$ are diagonal elements of ${\boldsymbol{K}_s},{\boldsymbol{K}_f}$ , respectively.

After obtaining the optimal designable parameters by GWO intelligent optimisation algorithm, the proposed ADP-based SMFTC method, is called as ADP-based intelligent SMFTC, and it is referred to as “ADP-based ISMFTC” method for short in the following.

The block diagram of the presented ADP-based ISMFTC method for aircraft engine with actuator dynamics and faults, is given in Fig. 4.

Figure 4. The block diagram of ADP-based ISMFTC method for aircraft engine with actuator dynamics and faults.