Fault-Tolerant Control Scheme for the Sensor Fault in the Acceleration Process of Variable Cycle Engine

Abstract

This paper presents a fault-tolerant control scheme for the sensor fault in the acceleration process of the variable cycle engine. Firstly, an adaptive equilibrium manifold model with multiple inputs and multiple outputs is established. Combined with the Kalman filter bank, sensor fault diagnosis is carried out to realize the diagnosis and signal reconstruction of the engine in the case of a single sensor and double sensor faults. On this basis, isolation and group isolation are used to diagnose sensor faults and reconstruct signal in speed closed-loop control. Then, the control plan of the acceleration process is optimized based on the target shooting method, aiming to simulate the variation of various variables in the engine acceleration process more accurately, so as to verify the feasibility of the sensor fault-tolerant control scheme. Finally, a hardware-in-loop simulation platform is built based on the idea of distributed control, and the fault-tolerant control scheme of the sensor proposed previously is verified based on this platform. The results show that the proposed scheme can accurately diagnose the sensor faults and reconstruct the signal within 0.2 s, and the actual speed can rise from 67.87% to 99.9% in 4 s, ensuring the safe and rapid completion of the acceleration process.

1. Introduction

As the “brain” of an aeroengine, the role of the control system is to make the engine work stably and reliably under any changing conditions and give full play to its performance benefits [1]. Visible, stability and reliability are the fundamental guarantee of all other higher control requirements. Modern aero engines adopt a full authority digital electronic control (FADEC) system responsible for complex tasks including control, monitoring and management. Its structure consists of numerous sensors, actuators and electronic control components. Due to the long-term work in high temperature, high pressure and intense vibration, these components are prone to failure during system operation. Once a fault occurs, the whole control system may collapse, resulting in severe consequences. According to the data, the sensor fault of aeroengine is higher than 80% of the total system fault [2].

Fault-tolerant control (FTC) means that when the system fails, the remaining control tasks can be completed safely according to the specified performance index, or the performance index can be reduced but still within the acceptable range. At present, the most fault-tolerant design of engine control systems is based on hardware redundancy. The development of non-hardware redundancy is mainly focused on sensor fault diagnosis and signal reconstruction. In the 1970s, Wallhagen, R. E. and D. J. Arpasi proposed a passive FTC technology with analytic redundancy to reconstruct the fault system for engines with sensor faults [3]. Since the 1980s, NASA has designed an analytical redundancy scheme to improve engine control reliability, which can detect and distinguish between hard and soft faults of the control system and reconstruct the system for sensor faults [4]. Around the 90s, Merrill et al. used a set of Kalman estimators in fault detection and isolation. By monitoring the residuals of each estimator, specific faults that have occurred can be detected and isolated [5,6]. Since the end of the last century, intelligent algorithms have been applied in sensor fault diagnosis and FTC systems and successfully verified in turbojet, turbofan and other engines [7,8]. Xianghua Huang used a self-association neural network to construct the analysis redundancy of sensors to ensure the reliability of aircraft engine closed-loop control [9]. Yang et al. used wavelet transform to diagnose sensor faults, and used BP neural network to predict the output data of fault sensors [10]. Zhao et al. obtained a set of modal functions by empirical mode decomposition of the output signal of the sensor, and took its sample entropy as the fault feature of the sensor, and classified the fault state of the sensor by sparse representation-based classifier [11]. Chen and Liu used radial basis function (RBF) neural network and general regression neural network (GRNN) respectively to predict timing data. If the difference between the predicted value and the actual output value of the sensor exceeds the set threshold, the sensor is considered to be faulty, and the sensor signal is reconstructed with the predicted value [12,13]. Sheng et al. proposed an intelligent FTC system based on an Online Sequential Extreme Learning Machine (OS-ELM) for sensor faults. The system can realize online fault diagnosis and signal reconstruction without establishing a system model [14]. These latest research results are mostly combined with artificial intelligence algorithm for sensor fault diagnosis, most of them stay at the theoretical level. The complexity of the algorithm reduces the speed of diagnosis and makes it difficult to apply it in engineering practice. Besides, most of the research on sensor FTC is based on steady-state control. In the process of engine transition, because the speed and other parameters are sensitive parameters closely related to the engine performance, to improve the stability and reliability of the transition, closed-loop control is often carried out. Therefore, the sensor FTC problem in the transition process also exists, but the relevant research is still relatively rare.

Because of this, this paper presents an acceleration control scheme based on fault tolerance of the sensor for variable cycle engine (VCE), as shown in Figure 1. First, in Section 2, an improved multi-input multi-output adaptive equilibrium manifold model (AEMM) with dual scheduling variables is established. In Section 3, the AEMM is combined with the Kalman filter bank to design the single sensor fault and the double sensor faults diagnosis scheme. Based on this, the actual engine parameter estimation under the sensor fault is obtained. In Section 4, A multi-target method is used to optimize the engine acceleration process control plan to verify the FTC scheme. In Section 5, the above procedure is verified by hardware-in-the-loop simulation, and the FTC of the accelerated process is proved from the perspective of engineering realizability. Section 6 concludes with a summary of the paper.

2. On-Board Adaptive Model

2.1. Equilibrium Manifold Model

Equilibrium manifold model (EMM) is a model developed from the linearization model [15].

Consider the following smooth nonlinear system:

$$\left\{\begin{array}{l}\dot{x}=f(x,u)\\ y=g\left(x,u\right)\end{array}\right.,$$

(1)

where $f(x,u)$ and $g(x,u)$ are smooth functions. Do Taylor expansion at any equilibrium point and abandon the higher-order terms, and a linear model of the original system can be obtained:

$$\left\{\begin{array}{l}\dot{x}=A\left(x-x_e\right)+B\left(u-u_e\right)\\ y-y_e=C\left(x-x_e\right)+D\left(u-u_e\right)\end{array}\right..$$

(2)

The set of all equilibrium points in a system constitutes an equilibrium family ${\Gamma }_e$ and ${\Gamma }_e$ can be considered an equilibrium manifold.

$${\Gamma }_e=\left\{\left(x_e,u_e,y_e\right)|f\left(x_e,u_e\right)=0,y_e=g\left(x_e,u_e\right)\right\}.$$

(3)

These equilibria can be scheduled with a variable of the same dimension as the control quantity $u$, which is defined as the scheduling variable $\alpha$. $\alpha$ is composed of relevant variables of the current working point and can be expressed as $\alpha =p\left(x,y,u\right)$. Then we can get:

$$\left\{\begin{array}{l}{x}_e=x_e\left(\alpha \right)\\ u_e=u_e\left(\alpha \right)\\ y_e=y_e\left(\alpha \right)\end{array}\right..$$

(4)

Similarly, the state matrix and control matrix of the linearized model of the system at each equilibrium point can also be expressed as the matrix scheduled with $\alpha$, so as to obtain the dynamic model of the system near the equilibrium manifold:

$$\left\{\begin{array}{l}\dot{x}=A\left(\alpha \right)\left(x-x_e\left(\alpha \right)\right)+B\left(\alpha \right)\left(u-u_e\left(\alpha \right)\right)\\ y-y_e\left(\alpha \right)=C\left(\alpha \right)\left(x-x_e\left(\alpha \right)\right)+D\left(\alpha \right)\left(u-u_e\left(\alpha \right)\right)\end{array}\right..$$

(5)

To make the linearized model family (5) a nonlinear model, it is necessary to establish the relationship between scheduling variables $\alpha$ and the current working point. Considering the algebraic relationship between $y$ and $x$, $u$, only the following design is needed:

$$\alpha =p\left(x,u\right).$$

(6)

By combining Equations (5) and (6), a nonlinear model scheduled by $x$ and $u$ can be obtained, that is, an equilibrium manifold model. It can be seen that the EMM is derived from the linearized model family (5) and the scheduling relation (6). It has the characteristic of variable parameter of linear parameter varying (LPV) model whose coefficient matrix all changes with the working state, and the introduction of Equation (4) makes the model always change near the equilibrium manifold, which overcomes the deficiency of steady-state error in LPV model.

In terms of the structure of the model, the EMM can be regarded as a Taylor expansion model at the equilibrium point of real-time change. The Equation (6) essentially establishes a mapping of the current operating point to the Taylor expansion point on the equilibrium manifold. Figure 2 shows the three mapping modes of equal input $u=u_e\left(\alpha \right)$, equal partial state $x_i=x_{ie}\left(\alpha \right)$, and the orthogonal expansion $\sqrt{{\left(x_e\left(\alpha \right)-x\right)}^2+{\left(u_e\left(\alpha \right)-u\right)}^2}$. The Equation (6) can be obtained from these mappings. The three designs have different EMM and different errors from Equation (1). Determine the value of the scheduling variable utilizing equal input shown in Figure 2 and select fuel flow $W_f$ and nozzle area $A_8$ as the scheduling variables. Then, there is $W_{fe}\left(\alpha \right)=W_f$, $A_{8e}\left(\alpha \right)=A_8$ and $\alpha =\left(W_f,A_8\right)$. Represent $A\left(\alpha \right)$ as $\left[\begin{array}{cc}{a}_{11}(\alpha )& a_{12}(\alpha )\\ a_{21}(\alpha )& a_{22}(\alpha )\end{array}\right]$ and the EMM of the above process can be expressed as:

$$\left\{\begin{array}{l}{\dot{n}}_h=a_{11}\left(\alpha \right)\left(n_h-n_{he}\left(\alpha \right)\right)+a_{12}\left(\alpha \right)\left(n_l-n_{le}\left(\alpha \right)\right)\\ {\dot{n}}_l=a_{21}\left(\alpha \right)\left(n_h-n_{he}\left(\alpha \right)\right)+a_{22}\left(\alpha \right)\left(n_h-n_{le}\left(\alpha \right)\right)\\ \alpha =u=\left(W_{fe},A_{8e}\right)\end{array}\right..$$

(7)

It can be seen that, due to the selection of scheduling variables and scheduling relations, the control matrix $B\left(\alpha \right)$ will not appear in the equilibrium manifold expansion model, and the EMM is reduced to a group of differential equations without initial input. It should be noted that the complete dynamic matrix can be obtained by taking partial derivatives of the corresponding variables [16]. In other words, the EMM and the nonlinear model have the same state-space matrix at the same equilibrium point, that is, their linearization models are the same. Using this property, the B matrix at any equilibrium point can be obtained.

Due to the strong nonlinear characteristics of the engine, just as different input combinations can obtain different linearized models, a nonlinear object also has multiple equilibrium manifold dynamic structures, which means different combinations of inputs give different $A\left(\alpha \right)$. That is to say, if a set of inputs is applied at a certain equilibrium point and identification data is obtained, then the identified model can only be used in the neighborhood of the set of inputs (i.e., the point with little state change). In addition, since the EMM loses the B matrix and only retains a matrix, different input combinations are easy to identify and obtain a matrix with great differences. Compared with the traditional MIMO linearization model, it is easier to generate dynamic characteristics with remarkable differences. To solve this problem, one way is to build different EMMs for different input combinations. We can define the following variables:

$$\theta =\frac{\Delta W_f/{\displaystyle \ell }{W}_f}{\Delta W_f/\ell W_f+\Delta A_8/\ell A_8},$$

(8)

where ${\displaystyle \ell }{W}_f$ and $\ell A_8$ represent the length of the subinterval, $\Delta W_f$ and $\Delta A_8$ represent the change value of the input quantity. Therefore, this variable can be regarded as the proportion of the change of fuel flow in the shift of total input quantity and can be called the modification of control quantity.

Then the EMM of the engine with a wide range of variable cycles can be obtained through identification and the parameterized representation of the model coefficient matrix under each $\theta$ [17]. It should be noted that parameter identification of coefficient matrix needs to be based on nonlinear component level model (CLM) of the engine. The establishment method of CLM of VCE has been very mature. For example, detailed steps are given in reference [18], which will not be elaborated in detail in this paper.

The following is to test the operation effect of the model by changing the control quantity randomly in an extensive range. The control quantity shown in Figure 3 is applied to the CLM and the EMM. The operation results of the two models are shown in Figure 4. It can be seen that the model has high steady-state accuracy and good dynamic accuracy.

2.2. On-Board Adaptive Equilibrium Manifold Model

Referring to the conventional linear state variable model with degradation, the following EMM with degradation is established:

$$\left\{\begin{array}{l}\dot{x}=A\left(\alpha \right)\left(x-x_e\left(\alpha \right)\right)+L\left(\alpha \right)h\\ y-y_e\left(\alpha \right)=C\left(\alpha \right)\left(x-x_e\left(\alpha \right)\right)+M\left(\alpha \right)h\end{array}\right.,$$

(9)

where h is the engine health parameter, representing the degradation of the flow/efficiency of each engine component. $\alpha$ is the scheduling variable. Since the health parameters are independent of the scheduling variables, they are not scheduled by them, while the coefficient matrices L and M are related to the engine operating points and are scheduled by the scheduling variables.

Considering that performance degradation is a slow development process during the service period of the engine, it does not require high dynamic characteristics and focuses on its steady-state value. Therefore, based on the EMM without degradation, the matrix expressions are retained, and the health parameters are changed in the component level model to cause changes in state quantities and output quantities. The degradation coefficient matrices L(α) and M(α) are obtained by the changing steady-state values identification to form the final EMM with degradation.

An EMM with degradation is not enough to be called an adaptive model. It needs to have the ability to update health parameters (degradation) to match with the actual engine running state. For this reason, a Kalman filter is introduced to estimate the starter degradation from the measurable output deviation containing the measured noise. Firstly, model (9) is written into a standard form. For this reason, the health parameters are extended to state quantities.

$$\left\{\begin{array}{l}\left[\begin{array}{c}\dot{x}\\ \dot{h}\end{array}\right]=A_{aug}\left[\begin{array}{c}x-x_e\left(\alpha \right)\\ h\end{array}\right]+\omega \\ y-y_e\left(\alpha \right)=C_{aug}\left[\begin{array}{c}x-x_e\left(\alpha \right)\\ h\end{array}\right]+\upsilon \end{array}\right.,$$

(10)

where $A_{aug}=\left[\begin{array}{cc}A\left(\alpha \right)& L\left(\alpha \right)\\ 0& 0\end{array}\right]$ and $C_{aug}=\left[\begin{array}{cc}C\left(\alpha \right)& M\left(\alpha \right)\end{array}\right]$. $\omega$ and $\upsilon$ are the system noise matrix and measurement noise matrix respectively, and the covariance matrices are $Q$ and $R$. Then, the Kalman filter uses the probability distribution of the measured value and the optimal estimate to obtain the optimal state estimate:

$$\left\{\begin{array}{l}\left[\begin{array}{c}\dot{\widehat{x}}\\ \dot{\widehat{h}}\end{array}\right]=A_{aug}\left[\begin{array}{c}\widehat{x}-x_e\left(\alpha \right)\\ \widehat{h}\end{array}\right]+K\left(y-\widehat{y}\right)\\ \widehat{y}-y_e\left(\alpha \right)=C_{aug}\left[\begin{array}{c}\widehat{x}-x_e\left(\alpha \right)\\ \widehat{h}\end{array}\right]\end{array}\right.,$$

(11)

where $\widehat{x}$, $\widehat{y}$, $\widehat{h}$ are estimated values of state quantity, output quantity and degradation quantity. $K$ is called the Kalman filter gain, and $K=P{C}_{aug}{R}^{-1}$. $P$ is the solution of the following Riccati equation:

$$A_{aug}P+P{A}_{aug}^T-P{C}_{aug}^T{R}^{-1}{C}_{\mathrm{aug}}P+Q=0.$$

(12)

Thus, the AEMM is formed. The ability of health parameter estimation of this model is verified at the design point below. Set high-pressure turbine component efficiency in degrading by 2%, and the steady-state deviation of the output of the CLM and the output without degradation is respectively $\Delta n_l=-1.98\%$, $\Delta n_h=-1.93\%$, $\Delta P_{t3}=-52.84$ $\mathrm{Kpa}$ and $\Delta T_{t3}=-10.0$ K. The steady-state output is added with Gaussian noise and fed into the Kalman filter to obtain the optimal estimation of health parameters and measured values, as shown in Figure 5 and Figure 6.

It can be seen that the Kalman filter can accurately estimate the efficiency degradation of the high-pressure turbine components at 2% within 1.5 s, and the estimation accuracy of other health parameters is also within 0.1%. At the same time, the estimation of the measured values also has high accuracy.

3. Sensor Fault Diagnosis

In this section, the AEMM is combined with the Kalman filter bank to design a fault diagnosis scheme for single sensor and double sensors. In the case of the single sensor fault, filter banks are designed. Each filter isolates one of the sensors, and the output estimation is carried out by using the remaining sensor information. In the case of the two sensor faults, they are grouped and the single-sensor diagnosis mode is adopted.

3.1. Single Sensor Fault

Assuming that m measurable parameters are used to estimate engine performance, m Kalman filters are designed accordingly. Each filter isolates the measured value of one of the sensors, and the remaining $(m-1)$ measured values are used for state/output estimation. The structure is shown in Figure 7. Take the ith Kalman filter as an example. Its input contains the remaining $(m-1)$ measurement values except for the ith sensor. After the failure of the ith sensor, the estimated values of the other filters are deviated from the actual measured values due to the access of the fault sensor signal. When the deviation of all the measured values is added up in a certain way, a relatively significant change can be obtained, and the faulty sensor can be identified. Moreover, since the ith filter is not connected to the fault sensor, the error between its estimated value and the actual measured value is small so that it can be used for fault signal reconstruction. The Kalman filter for the ith sensor is:

$$\left\{\begin{array}{l}\left[\begin{array}{c}\dot{\widehat{x}}\\ \dot{\widehat{h}}\end{array}\right]=A_{aug}\left[\begin{array}{c}\widehat{x}-x_e\left(\alpha \right)\\ \widehat{h}\end{array}\right]+K^i\left(y^i-{\widehat{y}}^i\right)\\ {\widehat{y}}^i-y_e^i\left(\alpha \right)=C_{aug}^i\left[\begin{array}{c}\widehat{x}-x_e\left(\alpha \right)\\ \widehat{h}\end{array}\right]\end{array}\right.$$

(13)

where $y^i$ is the subset of measurement parameters excluding the ith sensor, matrix $C_{aug}^i$ is the submatrix of matrix $C_{aug}$ after removing the ith row, and K is the gain of the filter.

Define the residual between the filter input $y^i$ and the optimal output estimate ${\widehat{y}}^i$:

$$e^i=y^i-{\widehat{y}}^i$$

(14)

when the sensor is fault-free and the filtering tends to be stable, the residual should obey the multi-dimensional standard normal distribution. However, when a fault occurs, the residual error does not possess this property, and its expectation is often not zero. Define the following variables:

$${\mathit{WSSR}}^i=W_r^i{\left(e^i\right)}^T{\left({\sigma }_i^2\right)}^{-1}\left(e^i\right)$$

(15)

where WSSR is the weighted sum of squared residual. When there is no fault, WSSR obeys the ${\chi }^2$-distribution of $(m-1)$ degrees of freedom. $W_r^i$ is the weight, and it can adjust the output of the filter. Its value will affect the diagnostic effect, resulting in misdiagnosis and missed diagnosis.

Select sensor measurement parameters as $y=(n_l,n_h,P_{t3},T_{t3},P_{t21},T_{t21},P_{t6},T_{t6})$. Compared with AEMM established previously, the four additional sensors are the total pressure and temperature at the fan outlet and the total pressure and temperature at the mixing chamber inlet. State quantity and health parameters remain unchanged. The noise matrices Q and R are reselected, and 8 Kalman filters are designed. The input of each filter is another 7 measured values. Set the efficiency of high-pressure turbine components to degrade by 2%, and obtain the steady-state deviations of 8 sensors under normal conditions: $\Delta n_l=-1.98\%$, $\Delta n_h=-1.93\%$, $\Delta P_{t3}=-52.84\mathrm{Kpa}$, $\Delta T_{t3}=-10.0\mathrm{K}$, $\Delta P_{t21}=-0.28\mathrm{Kpa}$, $\Delta T_{t21}=-2.35\mathrm{K}$, $\Delta P_{t6}=-2.96\mathrm{Kpa}$, $\Delta T_{t6}=21.18\mathrm{K}$. At 2 s, a certain bias signal is injected into the outlet temperature sensor of the high-pressure compressor ($T_3$). As can be seen from Figure 8, except for the fourth filter, the WSSR value of all the other filters underwent drastic abrupt changes and then gradually approached a stable value, while the WSSR value of the fourth filter remained at a low level. Therefore, it can be determined that the fourth sensor, the temperature sensor $T_{t3}$, has failed. Then, a drift fault of 1 K/s is injected into the sensor. As can be seen from Figure 9, the WSSR of the other filters also drifts, while only the fourth filter remains unchanged. Therefore, the designed Kalman filter bank can diagnose all kinds of single sensor faults effectively.

Subsequently, this filter is used to re-estimate the degradation, and the result is shown in Figure 10, which is in line with the actual situation. Finally, the actual degradation is brought into the on-board AEMM to reconstruct the fault sensor signal, as shown in Figure 11, where the short red line is the fault sensor signal, and the solid black line is the corresponding sensor signal reconstructed from the on-board model.

3.2. Dual Sensor Faults

Although the probability of simultaneous failure of two or more sensors is relatively low, due to the high-reliability requirements of engine control and fault diagnosis system, on the basis of being able to diagnose a single fault accurately, this section further studies the fault diagnosis and signal reconstruction of two sensors.

When the dual-sensor fault occurs, the Kalman filter design scheme that isolates one of the sensors mentioned above will not work because all the filters use the fault information. The output residual will make the WSSR of all the filters change. In this case, a conservative but effective scheme is to isolate the two sensors, but the number of Kalman filters required will increase as the number of combinations. So, when there are more sensors, the idea of grouping sensors can be used to design. For example, there are 2 N sensors, which are divided into groups A and B. Each group has N sensors. Kalman filter banks are independently designed for these two groups to carry out the single fault diagnosis. If both sensor faults occur in one group, the WSSR values of the filters in this group will deviate from the normal values, while the other group will keep low values. In this case, the signals in this group can be completely abandoned, and the optimal state estimation and signal reconstruction can be carried out with another group of filters. If sensor faults occur in both groups A and B, only one filter will remain low in both groups. The two fault sensors can be diagnosed and isolated at this time, and the remaining sensors can be used for optimal state estimation and signal reconstruction. Figure 12, Figure 13, Figure 14 and Figure 15 are the results of adopting this scheme. Group A contains the signals as $y_A=(n_l,n_h,P_{t3},T_{t3},P_{t21},T_{t21})$ while $y_B=(P_{t3},T_{t3},P_{t21},T_{t21},P_{t6},T_{t6})$. Firstly, the fault sensors are set as the high spool speed sensor ($N_h$) and the total pressure sensor at the inlet of the mixing chamber ($P_{t6}$). As shown in Figure 12 and Figure 13, group A and group B each have a filter with a low WSSR value so that the fault sensor can be easily diagnosed. Secondly, the fault sensors are set as the total pressure sensor ($P_{t6}$) and temperature sensor at the inlet of the mixing chamber ($T_{t6}$). As shown in Figure 14 and Figure 15, it can be seen that the WSSR of all the filters in Group A remains low, while all the filters in Group B are offset. Therefore, it is judged that both sensor faults occur in Group B, and all sensors in Group A are normal. The methods for optimal state estimation and signal reconstruction using normal sensors are similar to those in the previous section and will not be repeated here.

4. Optimization Design of Acceleration Control Plan

4.1. Optimization Scheme of Control Plan Based on Target Shooting Method

In the 1980s, H.G.Bock proposed the idea of a direct multiple target shooting method, whose essence is to transform the optimal control problem into a parameter optimization problem [16,17]. Consider the following optimal control problem:

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

(16)

The variation range of control variables and state variables is constrained as:

$$\left\{\begin{array}{l}{l}_1\le u\left(t\right)\le h_1\\ l_2\le x\left(t\right)\le h_2\end{array}\right..$$

(17)

For initial state and terminal state, there are the following boundary conditions:

$$\left\{\begin{array}{l}x(0)=x_0\\ x(T)=x_T\end{array}\right..$$

(18)

The optimal performance index is:

$$J={\displaystyle {\int }_0^{T}L\left(x\left(t\right),u\left(t\right),t\right)dt}.$$

(19)

Terminal time T can be fixed or free. For the case where T is free, a new time variable $\tau \in \left[0,1\right]$ can be introduced. Let T be a variable parameter in a range and define $t=T\tau$, $\frac{dt}{d\tau }=T$ and $\frac{dT}{d\tau }=0$. At the same time, the augmented state variable $x_{l+1}$ and state equation are introduced:

$${\dot{x}}_{l+1}=T\cdot L\left(x\left(\tau \right),u\left(\tau \right),T\tau \right),x_{l+1}\left(0\right)=0.$$

(20)

Thus, the above optimal trajectory problem is equivalent to the following Mayer problem [18] with fixed terminal time.

$$\min f_1\left(x_{l+1}\left(1\right),T\right).$$

(21)

The Equation of state of the system is as follows:

$$\dot{x}=Tf\left(x\left(\tau \right),u\left(\tau \right),T\tau \right)=F\left(x\left(\tau \right),u\left(\tau \right),T,\tau \right).$$

(22)

The starting and ending constraints become:

$$\left\{\begin{array}{l}x(0)=x^0\\ x(1)=x^f\end{array}\right..$$

(23)

Control variables, state variables, and process time constraints are shown in the following Equation:

$$\left\{\begin{array}{l}{L}_1\le u\left(\tau \right)\le H_1\\ L_2\le x\left(\tau \right)\le H_2.\\ L_3\le T\le H_3\end{array}\right.$$

(24)

Next, the Mayer problem can be solved by using a direct multiple target shooting method. The basic idea of the target shooting method is to transform the infinite dimensional problem into a finite dimension problem. The optimization process is divided into several intervals, and the design values of control variables and state variables are given in the interval nodes. In this way, the constraints on the whole process state and control variables can be transformed into constraints on these design values, and the performance indicators can also be transformed into functions of the design values. The following matching conditions should be considered: Set the change rule of the control quantity in each interval. It can be in the form of a constant value or polynomial. Then, the state variables at the end of each interval can be obtained from the guess value of the state variables at the beginning of the interval and the system state equation as shown in Equation (22). This value needs to be equal to the estimated value of the state here, which will be met by continuously revising the design value through the optimization algorithm. The specific steps of the target shooting method are as follows:

(1)

Divide the fixed time interval [0, 1] into m equal parts to get the node: ${\tau }_i=\frac{i}{m}$, where $i=0,1,\dots ,m$.

(2)

Parameterize the control quantity. Introduce a set of vectors $P_i\in R^q$ where $i=0,1,\dots ,m$ and define:

$$u\left(\tau \right)=u_i\left(P_0,P_i,\dots ,P_m\right),\tau \in \left[{\tau }_i,{\tau }_{i+1}\right],$$

(25)

where $u_i\left(\tau \right)$ can be a constant value or some interpolation function.

(3)

Set up the initial value problem. A set of vectors $S_i\in R^n$ is selected as an estimate of the state variables $x_i\left({\tau }_i\right)$ at the node ${\tau }_i$. Then we have m initial value problems (IVP):

$${\dot{x}}_i\left(S_i,u_i,T,\tau \right)=F\left(x_i\left(S_i,u_i,T,\tau \right),u_i\left(\tau \right),T\right),\tau \in \left[{\tau }_i,{\tau }_{i+1}\right],$$

(26)

where the initial value is:

$$x_i\left(S_i,u_i,T,{\tau }_i\right)=S_i,i=0,1,\dots ,m-1.$$

(27)

(4)

Constituting Nonlinear Programming (NLP). The optimization objective is:

$$\min f_1\left(S_m^{*},T\right),$$

(28)

where $S_m^{*}$ is the estimate of $x_{l+1}\left(1\right)$. The endpoints of each interval meet the matching conditions:

$$x\left(S_i,u_i,T,{\tau }_i\right)-S_{i+1}=0,i=0,1,\dots ,m-1.$$

(29)

The starting and terminal states meet the boundary conditions:

$$\left\{\begin{array}{c}{S}_0=x_0\\ S_m=x_T\end{array}\right..$$

(30)

Constraints on state quantity, control quantity and terminal time in the process are transformed into constraints on parameters at the end of the interval:

$$\left\{\begin{array}{l}{l}_1\le P_i\le h_1\\ l_2\le S_i\le h_2, i=0,1,\dots ,m\\ l_3\le T\le h_3\end{array}\right..$$

(31)

Finally, the Mayer optimization problem is transformed into a general NLP problem as shown in the following Equation.

$$\left\{\begin{array}{l}\min F_1\left(y\right)\\ F_2\left(y\right)=0,\\ F_3\left(y\right)\le 0\end{array}\right.$$

(32)

where $F_1$, $F_2$ and $F_3$ respectively represent the objective, equality constraint and inequality constraint function, and the optimization variable $y$ is:

$$y={\left(S_0^T,S_1^T,\dots ,S_m^T,P_0^T,P_1^T,\dots ,P_m^T,T\right)}^T.$$

(33)

For the solution of the NLP problem, any NLP algorithm can be used, such as sequential quadratic programming, active set, interior point method, etc. In this paper, the Quasi-Newton method is used. The optimization variables are iteratively updated by the following formula:

$$y^{k+1}=y^k+t_k\Delta y^k,t_k\in \left[t_{\min },t_{\max }\right],$$

(34)

where the increment $\Delta y^k$ is determined by the Kuhn-Tucker condition of the quadratic approximation problem below:

$$\left\{\begin{array}{l}\min \frac{1}{2}\Delta y^T{B}^k\Delta y+\nabla F_1^k\Delta y\\ F_2^k+\nabla F_2^k\Delta y=0\\ F_3^k+\nabla F_3^k\Delta y\ge 0\end{array}\right.,$$

(35)

where $F_i^k$ and $\nabla F_i^k(i=1,2,3)$ are the values and gradients of the function at $y^k$, and $B^k$ is the approximation of the Hessian matrix (${\nabla }_y^{2}L$) of the following Lagrange equation:

$$L\left(y,\lambda \right)=F_1\left(y\right)-{\lambda }_1^T{F}_2\left(y\right)-{\lambda }_2^T{F}_3\left(y\right).$$

(36)

The above parameters can be updated by the following formula:

$$\left\{\begin{array}{l}{B}^{k+1}=B^k+U\left(B^k,\Delta y^k,{\gamma }^k\right)\\ {\gamma }^k=\nabla L\left(y^{k+1},{\lambda }_1^k,{\lambda }_2^k\right)-\nabla L\left(y^k,{\lambda }_1^k,{\lambda }_2^k\right)\end{array}\right.,$$

(37)

where $U\left(B^k,\Delta y^k,{\gamma }^k\right)$ represents some correction method. The more commonly used is BFGS correction.

4.2. Accelerated Process Control Plan Design of VCE

Acceleration is the process of an engine from one thrust level to another. An important index to measure the performance of the acceleration process is acceleration time, which is mainly determined by the residual power of the engine. The main control variables of the acceleration process are the fuel flow and the nozzle area, among which the fuel quantity is the biggest factor affecting the engine acceleration performance. Increasing the fuel flow can increase the gas temperature in front of the turbine ($T_{t4}$), thus increasing the power capacity of the turbine. However, it should be noted that $T_{t4}$ is limited by the strength of turbine blades and oil-rich flameout of the combustion chamber, and the sharp increase in fuel quantity will also cause engine speed to exceed limit and pose a certain threat to the mechanical load of the rotor. On the other hand, the acceleration process will lead to the decline of the compressor stability margin, so the compressor stability margin is often increased by opening the nozzle, deflating air, adjusting the guide vane of the compressor, etc.

Select the starting point of the accelerated process state $x_0$ and the control quantity $u_0$:

$$\left\{\begin{array}{l}{x}_0:n_l=67.87\%,n_h=84.37\%\\ u_0:W_f=0.33 \mathrm{kg}/\mathrm{s},A_8=0.088 {\mathrm{m}}^2\end{array}\right..$$

(38)

Set the terminal state quantity $x_T$ and terminal control quantity $u_T$:

$$\left\{\begin{array}{l}{x}_T:n_l=99.90\%,n_h=98.50\%\\ u_T:W_f=0.80 \mathrm{kg}/\mathrm{s},A_8=0.096 {\mathrm{m}}^2\end{array}\right..$$

(39)

The constraints in the optimization process are as follows:

$$\left\{\begin{array}{l}{n}_h<100\%\\ n_l<100\%\\ T_{t4}<1850 \mathrm{K}\\ {\mathrm{SM}}_{\mathrm{CDFS}}>5\%\\ {\mathrm{SM}}_{\mathrm{HPC}}>5\%\\ W_f<0.84 \mathrm{kg}/\mathrm{s}\\ A_8<0.1 {\mathrm{m}}^2\\ {\dot{W}}_f<0.25 \mathrm{kg}/{\mathrm{s}}^2\\ {\dot{A}}_8<0.01 {\mathrm{m}}^2/\mathrm{s}\end{array}\right.,$$

(40)

where the last two terms represent the limit on the rate of change of the control quantity. The optimization objective is the shortest acceleration time:

$$\min J={\displaystyle {\int }_0^{T}dt}.$$

(41)

The AEMM is used as the optimization model of the accelerated process. Then the problem is modified to apply it to the target shooting method. Divide acceleration process into m intervals and the optimization variables $y={\left(S_0^T,S_1^T,\dots ,S_m^T,P_0^T,P_1^T,\dots ,P_m^T,T\right)}^T$ can be obtained. The starting and ending constraints can then be translated into the following equality constraints:

$$F_1{(y)}_1=\left\{\begin{array}{c}{S}_0-x_0\\ S_m-x_T\end{array}=0.\right.$$

(42)

In addition, the state quantity of the starting point and the end point of each subinterval should be guaranteed to match.

$$F_1{(y)}_2=x\left(S_i,u_i,T\right)-S_{i+1}=0,i=0,1,\dots ,m-1.$$

(43)

These two expressions together form equality constraints $F_1(y)$. Note that the fourth-order Runge-Kutta method is used here to calculate the state quantities at the end of each interval. In the process, the restrictions on the high and low spool speed, the gas temperature in front of the turbine, the surge margin, the change rate of the control quantity, etc. are also easily expressed in the form of $F_2(y)<0$. For example, the limit on the rate of change of fuel flow can be written as:

$$F_2{(y)}_1=\frac{(P_{i+1,W_f}-P_{i,W_f})\cdot T}{m}-0.3<0, i=0,1,\dots ,m-1,$$

(44)

where $P_{i,W_f}$ represents the parameterized results of fuel flow at each node. The optimization goal is transformed into:

$$\min F_3(y)=T.$$

(45)

At this point, the problem has been transformed into an NLP problem as shown in Equation (32).

Solve the above problems, and the results are shown in Figure 16, Figure 17 and Figure 18, where CDFS represents the core drive fan stage of VCE.

It can be seen that the low spool speed changes from 67.87% to 99.9%, and the high spool speed changes from 84.37% to 98.5%. The acceleration time is about 2.4 s. The acceleration process is fast and stable without overshoot, and no overshoot occurs. At the beginning of acceleration, the dominant limiting condition is the rate of change of the control quantity, while at the later stage of acceleration, the dominant limiting condition is $A_8$ and the upper limit of $T_4$. On the whole, the proposed accelerated process control plan meets the expected objectives.

4.3. Sensor Fault Simulation

Use the variation rule of low spool speed and nozzle throat area optimized by target shooting method in 4.2 as control instructions for the acceleration process, and the acceleration control scheme of low spool speed closed-loop and nozzle throat area open-loop is adopted. The Kalman filter bank is designed by selecting sensors of low spool speed, high spool speed, high-pressure compressor outlet total pressure and total temperature, which constitute the sensor fault diagnosis module. Next, the sensor fault simulation in the acceleration process is carried out. The result is shown in Figure 19, Figure 20, Figure 21 and Figure 22.

At 3 s, the acceleration control begins. At about 3.8 s, 4% bias fault occurres in the low spool speed sensor, resulting in a large difference from the command speed. Therefore, the controller rapidly increases the amount of fuel to reduce the error between the two, as shown in Figure 22a. It can be seen from the changes of high spool speed and turbine front temperature that this adjustment leads to great changes in the actual engine state. From Figure 19, at 3.8 s the WSSR value of $n_h$, $P_{t3}$, $T_{t3}$ filters change greatly, and only the $n_l$ filter keeps a low value, so it is quickly determined that the low spool speed has been faulty. It is worth noting that the WSSR values of each filter also changes after 3 s, which is mainly caused by the difference between the on-board model and the real engine in the acceleration process. Compared with the change caused by sensor faults, the value is much smaller, so it does not affect the diagnosis. At about 4 s, the controller uses the analytical redundancy of low spool speed provided by the on-board model, and it can be seen that the engine states are restored. At about 5.4 s, the low spool speed output by the on-board model successfully keeps up with the control command, and the engine completes acceleration. In general, the scheme successfully achieves acceleration control under sensor faults. However, the sudden change of the control quantity caused by the sensor failure leads to the great change of the engine state, which is mainly because the actuation characteristics and small closed-loop control characteristics of the actuator are not considered in the digital simulation. Therefore, the hardware-in-the-loop test is used to further verify the results.

5. Hardware-in-the-Loop Simulation and Verification

The hardware-in-the-loop (HIL) test platform is shown in Figure 23 and Figure 24, which consists of an engine simulator, a distributed controller and a monitoring computer. The engine simulator includes NI myRIO and interface simulator. NI myRIO is responsible for the real-time operation of the engine model and simplified actuator model, and the interface simulator is responsible for simulating the signals of various sensors of the engine and collecting the current signals of actuator valves. The distributed controller comprises several intelligent nodes, including high and low spool speed acquisition node, pressure acquisition node, fuel flow control node, nozzle area control node, core control node and backup. Each intelligent node is an embedded system based on the ZYNQ development board and communicates with each other via the TTP/C bus. In addition to the tasks specified by TTP/C communication, each intelligent node also needs to complete the duties prescribed by the control system, which is conducted by two CPUs in the logical processing unit and the program processing unit in ZYNQ. The nodes are programmed and developed by Vivado’s integrated development environment and SDK software development suite. The monitoring computer is responsible for sending control instructions and monitoring the engine state, and is produced by LabVIEW. In terms of time schedule, the simulation step size of the engine simulator is 20 ms. In the distributed control system, the step size of the data acquisition node and the core control node is also 20 ms, while the step size of the actuator control node is 5 ms. Data transmission of each node on the TTP/C bus has a specified time sequence and is scheduled by a pre-planned Message Description List (MEDL).

The low spool speed, the high spool speed, the total pressure and the total temperature at the outlet of the high-pressure compressor are also selected for fault diagnosis and design. FTC results are shown in Figure 25, Figure 26, Figure 27 and Figure 28. At 3 s, the acceleration control starts. At this time, the speed of the engine and the speed of the on-board model is both good for tracking instructions. At about 3.8 s, a 5% bias occurs to the low spool speed sensor, as shown by the red dash in Figure 26. The controller rapidly increases the fuel flow control instruction, as shown in the solid black line in Figure 28. However, different from the result in 4.3, since the change rate of the actuator has reached the limit at this time, the actual fuel volume does not increase dramatically but maintains the maximum rate of increase. As shown in Figure 25, the low spool speed bias is quickly diagnosed by the Kalman filter, and the diagnosis effect is noticeable. At about 4 s, the reconstructed signal of the airborne model is switched for FTC and the fuel flow control instruction falls back. At 6.8 s, the reconstructed low spool speed signal keeps up with the control instruction and tends to a steady state, and the engine accelerates successfully. The fault diagnosis and signal reconstruction can be done within 0.2 s, and the whole FTC process can be carried out in less than 3 s. Parameters before and after the acceleration process are shown in

Table 1. Parameters before and after the acceleration process.

Variable	Before the Acceleration Process	After the Acceleration Process
Low spool speed (%)	67.87	99.9
High spool speed (%)	84.37	98.5
Stall margin of CDFS (%)	7.62	10.13
Inlet temperature of HPT (K)	1435.25	1823.69
Fuel flow (kg/s)	0.33	0.80
High spool speed (%)	84.37	98.5

. In general, the proposed scheme has good fault diagnosis ability and fault tolerance ability. The acceleration control under sensor fault is successfully realized, and the acceleration process does not cause a big change of engine state. Therefore, the validity of the FTC scheme proposed in this paper and its compatibility with the real-time embedded platform has been proved.

6. Conclusions

In this paper, an FTC scheme for the sensor fault in the transition process of VCE based on AEMM is proposed. The scheme adopts closed-loop control of speed and open-loop control of other control variables for transition control. To diagnose the sensor faults and study the control of acceleration process, a modeling study is first carried out. The fuel flow and the nozzle area are selected as the two scheduling variables to establish the EMM of VCE with multiple inputs and multiple outputs. Then the AEMM is obtained by introducing health parameters and Kalman filter. After that, the AEMM is combined with the Kalman filter bank to design a fault diagnosis scheme for the single sensor fault and double sensor faults. The simulation results of bias and drift fault cases show that the proposed scheme can accurately diagnose the fault sensors except when two fault sensors are grouped in the same group. But in either case, the scheme can accurately estimate the state and reconstruct the signal. Then, the target shooting method is used to optimize the engine acceleration control plan, and based on this, the fault simulation of the low spool speed sensor is carried out. Finally, a HIL simulation platform is built to verify the FTC method considering the actuation characteristics and other limitations of the actuator in the actual process. The results show that the fault sensor can be diagnosed and signal reconstruction can be carried out within 0.2 s by using the proposed scheme.

The method proposed in this paper can complete the FTC of sensor fault in a very short time, and is not too complicated, and has strong engineering feasibility and application value. At the same time, this method is also suitable for general turbofan engine fault-tolerant control. But there are still many areas that need to be further studied. For example, when two faulty sensors appear in the same group, only accurate signal reconstruction can be carried out, but the location of the faulty sensor cannot be diagnosed. And in addition to the acceleration process of the VCE, the FTC of the mode switching process should be paid more attention to. In this paper, only the fault-tolerant design of the sensor is carried out, and there may also be actuator fault, component fault and other forms. Therefore, the FTC problem of VCE needs more and more in-depth research.