Control-Oriented Modeling for Nonlinear MIMO Turbofan Engine Based on Equilibrium Manifold Expansion Model

Abstract

This paper investigates the control-oriented modeling for turbofan engines. The nonlinear equilibrium manifold expansion (EME) model of the multiple input multiple output (MIMO) turbofan engine is established, which can simulate the variation of high-pressure rotor speed, low-pressure rotor speed and pressure ratio of compressor with fuel flow and throat area of the nozzle. Firstly, the definitions and properties of the equilibrium manifold method are presented. Secondly, the steady-state and dynamic two-step identification method of the MIMO EME model is given, and the effects of scheduling variables and input noise on model accuracy are discussed. By selecting specific path, a small amount of dynamic data is used to identify a complete EME model. Thirdly, modeling and simulation at dynamic off-design conditions show that the EME model has model accuracy close to the nonlinear component-level (NCL) model, but the model structure is simpler and the calculation is faster than that. Finally, the linearization results are obtained based on the properties of the EME model, and the stability of the model is proved through the analysis of the eigenvalues, which all have negative real parts. The EME model constructed in this paper can meet the requirements of real-time simulation and control system design.

1. Introduction

In recent years, as an ideal power plant for aircraft, turbofan engine has received extensive attention [1,2,3,4,5]. The continuous improvement of aircraft performance and the expansion of flight envelope have put forward higher requirements for the control system of turbofan engines. Data show that 80% of the work for the aero-engine control system design is devoted to modeling and understanding the characteristics of engine [6]. Therefore, control-oriented modeling plays a significant role in the design. However, the turbofan engine is a complex aerodynamic thermodynamic system with multi-variable, nonlinear, and time-varying characteristics, so the control-oriented modeling of the turbofan engine faces great challenges.

The modeling needs to resolve the contradiction between model accuracy and real-time. Early work on aero-engine modeling is to build nonlinear component-level (NCL) models [7,8,9]. The input-output relationships are established based on the characteristics of each component of the turbofan engine. Then, the NCL model is established according to the relationship and various constraints. At present, the NCL model has been applied in real-time simulation and used to complete the real-time digital verification of the control system [10,11]. However, the NCL model cannot meet the requirements of the closed-loop control system design for the turbofan engine due to its complexity [12]. Since most modern control methods are based on linear theory, the state variable model is often used in linear system theory. It is a linear model established in the state space, which can be obtained by derivative method [13] and fitting method [14,15]. Unfortunately, although the state variable model solves the application of linear control method in turbofan engines, it only has good fitting effect near the steady-state point used in the model identification. This shows that the state variable model is not effective for a wide range of operating conditions. At the same time, it is difficult to guarantee the good performance of turbofan engines in full envelope or in a wide range of operating conditions. In order to deal with these problems, the linear parameters-varying (LPV) model has been widely discussed in academia [16]. The modeling methods of LPV and quasi-LPV are discussed in Refs. [16,17], and the main difference between them is that quasi-LPV considers the change of state variables. There are three main modeling methods for the LPV model: the Jacobian linearization approach [18], the state transformations techniques [19], and the function substitution method [20]. Huang et al. [21] studied the identification of the LPV model with two scheduling variables. Rotondo et al. [22] introduced the modeling of the quasi-LPV for MIMO system. Lu and Huang [23] proposed a new lifecycle real-time model to describe the dynamic behavior of turbofan engines and studied engine/model mismatch compensation and performance degradation adaptive problems. The LPV model was converted to a switched convex polytopic form with hysteresis switching logic, and a switched LPV model of aero-engine rotor speed was obtained in Ref. [24]. With the development of LPV technology, many studies have appeared on the control method of turbofan engines in a wide range of operating conditions. However, the main problem of LPV method is that equilibrium points of engines change greatly with the scheduling variable in the actual work process, which brings great challenges to the accuracy and stability of the LPV model.

In order to obtain a control-oriented model that guarantees the accuracy and stability in a wide range of working conditions, the equilibrium manifold expansion (EME) model is proposed [25,26,27,28,29]. The ultimate goal of the equilibrium manifold method is to obtain the linearized family of nonlinear plants expanded along the equilibrium manifold. Thus, the EME model contains steady-state and dynamic characteristics of the nonlinear system near the equilibrium manifold, which are beneficial for the modeling of turbofan engines with operating line. Different from the LPV model, the EME model considers the mapping between the EM and scheduling variables, thereby ensuring the high accuracy of the model. Moreover, the equilibrium manifold method has been well applied in the modeling of areo-engines. The comparison of simulation and engineering test data has a small error. Yu et al. [25] studied the design of the shock position control system based on the equilibrium manifold method. Sui and Yu [26] proposed the EME model, analyze the error of the model, and solve the ill-posed problem of the identification matrix by using the constraint conditions of the equilibrium manifold. Yu and Zhao [27,28] developed the aero-engine control based on the EME model, and establish the affine EME model to complete the feedback linearization control of the engine. Liu et al. [29,30] built a control-oriented surge margin EME model and discussed the adaptive problem based on the EME model. In order to further improve the accuracy of the EME model, a switched EME model was proposed, and the switching control of the aircraft engine was carried out after the stability was proven [31]. Chen and Zhao [32] applied the switching control to the acceleration and safety protection of turbofan engines based on the switched EME model. With the improvement of the performance requirements for turbofan engine, the engine must be developed in a geometrically adjustable direction to ensure both performance and safety. Unfortunately, the aforementioned EME models are identified from the SISO systems with fuel flow as an input variable, which cannot meet the control requirements of advanced turbofan engines in a wide range of working conditions. In particular, geometrically adjustable variables, such as throat area of the nozzle, expand the operating range of turbofan engines. Up to now, few studies about how to extend the EME model from SISO system to MIMO system have been reported. Therefore, constructing the MIMO EME mode of turbofan engine is a new perspective, which is the focus of this paper.

Compared with existing studies, the main contributions of this work are as follows. Firstly, the basic definition of EM and EME model is introduced in Section 2, and the mathematical properties of the EME model are given. Subsequently, the modeling of the EME model for MIMO turbofan engine is given in Section 3. The steady-state and dynamic data needed for the identification of the MIMO EME model are obtained through the NCL model simulation, and all the parameters are normalized. After constructing the EME model, the EME model of the MIMO turbofan engine is obtained by using the steady-state and dynamic two-step identification method. In the process of dynamic parameter identification, it is creatively proposed to use less identification data and discuss the influence of different paths on dynamic parameter identification. In addition, the EME model shows that it can ensure the recognition accuracy under the interference of input noise. The EME models of the MIMO turbofan engine with different scheduling variables are discussed. Besides, the validation and stability analysis of the EME model is stated in Section 4. This paper verifies the accuracy of the EME model at off-design points. The validity and reliability of EME model is demonstrated over the whole range of operating conditions. In addition, the stability analysis of the EME model is carried out based on the linearization results of the model. Finally, Section 5 presents conclusion.

2. Description of EM and EME Model

2.1. Definitions of the Model

Consider the nonlinear system in the following form:

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

(1)

where $x\in {\mathbf{R}}^n$ is the state variable, $u\in {\mathbf{R}}^m$ is the input variable, $y\in {\mathbf{R}}^r$ is the output variable, and $f\left(·\right)$ and $g\left(·\right)$ are corresponding smooth nonlinear functions.

Definition 1 

(Ref. [26]). The equilibrium manifold (EM) of the system (1) is a set of equilibrium points, namely:

$$\begin{array}{c}\hfill \left\{\left(x_e,u_e,y_e\right)\mid f\left(x_e,u_e\right)=0,y_e=g\left(x_e,u_e\right)\right\},\end{array}$$

(2)

where the subscript “e” represents “steady-state”.

The EM which contains a series of $n+m+r$ algebraic equations is parameterized by the scheduling variable $alpha$. Generally, the dimensional of scheduling variable $alpha$ is equal to the dimensional of input variable u, and the scheduling variable $alpha$ is selected from the steady-state value of input or state variables. Therefore, all the equations of the EM can be expressed as the functions of $alpha$, namely:

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

(3)

Definition 2 

(Ref. [27]). The equilibrium manifold expansion (EME) model of system (1) is a family of first-order Taylor expansion at different equilibrium points $\left(x_e\left(\alpha \right),u_e\left(\alpha \right)\right)$, and there appears a mapping $\alpha =h(x,u)$ between the scheduling variable $alpha$ and current operating point $\left(x,u\right)$. The EME model is expressed as:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\dot{x}={\left.\frac{\partial f}{\partial x}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}\left(x-x_e\left(\alpha \right)\right)+{\left.\frac{\partial f}{\partial u}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}\left(u-u_e\left(\alpha \right)\right)\hfill \\ y={\left.\frac{\partial g}{\partial x}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}\left(x-x_e\left(\alpha \right)\right)+{\left.\frac{\partial g}{\partial u}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}\left(u-u_e\left(\alpha \right)\right)+y_e\left(\alpha \right)\hfill \end{array},\right.\end{array}$$

(4)

where the content in Equation (4) can be simplified as:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\left.\frac{\partial f}{\partial x}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}=A\left(\alpha \right),{\left.\frac{\partial f}{\partial u}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}=B\left(\alpha \right)\hfill \\ {\left.\frac{\partial g}{\partial x}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}=C\left(\alpha \right),{\left.\frac{\partial g}{\partial u}\right|}_{x=x_e\left(\alpha \right),u=u_e\left(\alpha \right)}=D\left(\alpha \right)\hfill \end{array}\right].\end{array}$$

(5)

The research shows that the mapping is not arbitrary [26]. It must satisfy the condition that any equilibrium point maps to itself. That is, $\alpha =h\left(x_e\left(\alpha \right),u_e\left(\alpha \right)\right)$ must be satisfied in the mapping. Although the EME model of system (1) is obtained by Taylor expansion, it is still a nonlinear system when the scheduling variable is eliminated.

2.2. Properties of the Model

There are a lot of studies on equilibrium manifold modeling in SISO aircraft gas turbine engine. Compared with other models, the EME model has advantages in accuracy and real-time performance due to its property.

Property 1. 

The linearization results and the EME model (4) of nonlinear system (1) are the same at any equilibrium point [27].

Explanation:The EME model can be regarded as an approximate nonlinear result of system (1), which is valid only if the linearization results of the two are consistent. At the same steady state, the coefficient matrix $A\left(\alpha \right)$, $B\left(\alpha \right)$, $C\left(\alpha \right)$, and $D\left(\alpha \right)$ should be consistent.

Property 2. 

Equation (6) is a necessary and sufficient condition for EM (2) to become a family of linearized models [27].

$$\begin{array}{c}\hfill \left\{\begin{array}{c}A\left(\alpha \right)\frac{\partial x_e\left(\alpha \right)}{\partial \alpha }+B\left(\alpha \right)\frac{\partial u_e\left(\alpha \right)}{\partial \alpha }=0\hfill \\ C\left(\alpha \right)\frac{\partial x_e\left(\alpha \right)}{\partial \alpha }+D\left(\alpha \right)\frac{\partial u_e\left(\alpha \right)}{\partial \alpha }=\frac{\partial y_e\left(\alpha \right)}{\partial \alpha }\hfill \end{array}.\right.\end{array}$$

(6)

Explanation:Equation (6) is the constraint of the steady-state and dynamic structure of the EME model. The unique identification method and structure of the EME model may result in missing parameters of Equation (5). Therefore, the completeness of linearization results of the EME model can be guaranteed by the constraint of Equation (6).

3. Modeling of the EME Model for the MIMO Turbofan Engines

3.1. The Source of the Identification Data

In order to identify the EME model of the turbofan engine, it is necessary to carry out the simulation in a wide range of operating conditions. The NCL model used in this paper was established by Ref. [33]. A high degree of confidence component level model of a twinspool turbofan engine has been achieved. The turbofan engine mainly consists of the following components: inlet, fan, HPC, combustor, High Pressure Turbine (HPT), Low Pressure Turbine (LPT), bypass, and nozzle. Each component is modeled by aerothermodynamics calculations and solving a set of balance equations. The engine design operation data and characteristic maps of rotating components are used to construct the turbofan engine nonlinear model [34]. This NCL model of turbofan engine has been used in other studies [35,36,37,38]. The NCL model only provides the raw data needed for the EME model identification. Through the study of EME model identification methods, the high-precision mathematical models that can fit different types of turbofan engines will be obtained.

The EM of an MIMO turbofan engine develops from a simple curve to a complex surface form. In the simulation, the engine NCL model works at H = 0 km, Ma = 0. In this paper, input variables of MIMO turbofan engine are fuel flow $\left(W_f\right)$ and throat area of the nozzle $\left(A_n\right)$. At each $W_f$, the $A_n$ increases step by step during simulation. The specific change of input variables is shown in Figure 1. Then, the high-pressure rotor speed $\left(n_{HC}\right)$, low-pressure rotor speed $\left(n_{LC}\right)$, and engine pressure ratio ($EPR$) are considered to be the output variables of the turbofan engine. Figure 2, Figure 3 and Figure 4 show the output results of the NCL model in simulation. The inputs in Figure 1 is normalized, and the actual range of $W_f$ is 1.5906–2.4060 kg/s, while the actual range of $A_n$ is 0.2577–0.3144 ${\mathrm{m}}^2$. This process involves the acceleration of the engine. Meanwhile, the steady-state data which is required for the identification of the EME model is extracted from the dynamic process. In addition, all the data in this paper are normalized results.

Instead of simple data processing, it is significant to analyze the steady-state and dynamic performances of nonlinear plants and select the appropriate structure for the identification of the EME model. Firstly, when $W_f$ is small, $n_{HC}$ is hardly affected by $A_n$ in Figure 2. It indicates that the turbofan engine exhibits strong nonlinearity over a wide range of operating conditions. As the value of $W_f$ reaches around 0.3333, $n_{HC}$ rises with the increase of $A_n$ (greater than 0.6667). Subsequently, Figure 3 shows that both $W_f$ and $A_n$ have effects on $n_{LC}$ in 10 groups of dynamic operation. On the whole, larger $W_f$ and $A_n$ correspond to higher $n_{LC}$. As $W_f$ increases, there appears a prominent steady-state linear relationship between $A_n$ and $n_{LC}$. Next, Figure 4 indicates that $W_f$ effectively regulates the overall work capacity of the turbofan engine. The performance of the engine improves with the increase of $W_f$, and $EPR$ increases at the corresponding $A_n$. However, the engine is close to dangerous conditions in high performance, such as surge. Therefore, $A_n$ plays an important role in solving engine safety problems in high performance operation. Under the same $W_f$ condition, $EPR$ can be effectively reduced by increasing $A_n$, which provides control direction for avoiding compressor surge.

3.2. Identification of the MIMO EME Model

The turbofan engine is described as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\dot{X}=F(X,U)\hfill \\ Y=G(X,U)\hfill \end{array},\right.\end{array}$$

(7)

where $X={\left[n_{HC},n_{LC}\right]}^T$ is the vector of the state variable, $Y={\left[{EPR,n}_{HC},n_{LC}\right]}^T$ is the vector of the output variable, and $U={\left[W_f,A_n\right]}^T$ is the vector of the input variable. In this paper, the typical state variables $n_{HC}$ and $n_{LC}$ of turbofan engines are directly used as output variables.

The family of steady-state operating points of this engine is expressed as:

$$\begin{array}{c}\hfill \left\{\left(n_{HCe},n_{LCe},W_{fe},A_{ne},EPR\right)\mid \begin{array}{c}0=\dot{X}=F\left(n_{HCe},n_{LCe},W_{fe},A_{ne}\right)\\ Y_e=G\left(n_{HCe},n_{LCe},W_{fe},A_{ne}\right)\end{array}\right\},\end{array}$$

(8)

where $Y_e={\left[{EPR,n}_{HCe},n_{LCe}\right]}^T$.

The EM shown in the Equation (3) is parameterized by $\alpha =h\left(X,U\right)$. Thus, the EME model of the turbofan engine shown in Equation (4) can be obtained as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\dot{X}=A\left(\alpha \right)\mathsf{\Delta }X+B\left(\alpha \right)\mathsf{\Delta }U\hfill \\ Y=C\left(\alpha \right)\mathsf{\Delta }X+D\left(\alpha \right)\mathsf{\Delta }U+Y_e\hfill \\ \alpha =h(X,U)\hfill \end{array},\right.\end{array}$$

(9)

where $\mathsf{\Delta }X=X-X_e\left(\alpha \right)$, $\mathsf{\Delta }U=U-U_e\left(\alpha \right)$, $A\left(\alpha \right)=\left[\begin{array}{cc}{a}_{11}\left(\alpha \right)\hfill & a_{12}\left(\alpha \right)\hfill \\ a_{21}\left(\alpha \right)\hfill & a_{22}\left(\alpha \right)\hfill \end{array}\right]$, $B\left(\alpha \right)=\left[\begin{array}{cc}{b}_{11}\left(\alpha \right)\hfill & b_{12}\left(\alpha \right)\hfill \\ b_{21}\left(\alpha \right)\hfill & b_{22}\left(\alpha \right)\hfill \end{array}\right]$, $C\left(\alpha \right)=\left[\begin{array}{cc}{c}_{11}\left(\alpha \right)\hfill & c_{12}\left(\alpha \right)\hfill \\ c_{21}\left(\alpha \right)\hfill & c_{22}\left(\alpha \right)\hfill \\ c_{31}\left(\alpha \right)\hfill & c_{32}\left(\alpha \right)\hfill \end{array}\right]$, $D\left(\alpha \right)=\left[\begin{array}{cc}{d}_{11}\left(\alpha \right)\hfill & d_{12}\left(\alpha \right)\hfill \\ d_{21}\left(\alpha \right)\hfill & d_{22}\left(\alpha \right)\hfill \\ d_{31}\left(\alpha \right)\hfill & d_{32}\left(\alpha \right)\hfill \end{array}\right]$.

The number of scheduling variables $\alpha$ in the EME model is equal to the number of input variables. In this section, $n_{HC}$ and $A_n$ are selected as $\alpha$ to build the EME model of turbofan engine. So, the structure of $\alpha$ is shown as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\alpha }_1=n_{HC}=n_{HCe}\left(\alpha \right)\hfill \\ {\alpha }_2=A_n=A_n\left(\alpha \right)\hfill \\ \alpha =\left[{\alpha }_1,{\alpha }_2\right]\hfill \end{array}.\right.\end{array}$$

(10)

According to Equation (10), $\mathsf{\Delta }{n}_{HC}=n_{HC}-n_{HCe}\left(\alpha \right)=0$, and $\mathsf{\Delta }{A}_n=A_n-A_{ne}\left(\alpha \right)=0$. Meanwhile, $n_{HC}$ and $n_{LC}$ are both state variables and output variables. Thus, the EME model of turbofan engine shown in Equation (9) can be replaced by:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{n}}_{HC}=a_{12}\left(\alpha \right)\mathsf{\Delta }{n}_{LC}+b_{11}\left(\alpha \right)\mathsf{\Delta }{W}_f\hfill \\ {\dot{n}}_{LC}=a_{22}\left(\alpha \right)\mathsf{\Delta }{n}_{LC}+b_{21}\left(\alpha \right)\mathsf{\Delta }{W}_f\hfill \\ EPR=EP{R}_e\left(\alpha \right)+c_{12}\left(\alpha \right)\mathsf{\Delta }{n}_{LC}+d_{11}\left(\alpha \right)\mathsf{\Delta }{W}_f\hfill \\ n_{HC}=n_{HCe}\left(\alpha \right)\hfill \\ n_{LC}=n_{LCe}\left(\alpha \right)+\mathsf{\Delta }{n}_{LC}\hfill \\ \alpha =\left[{\alpha }_1,{\alpha }_2\right]=\left[n_{HC},A_n\right]\hfill \end{array}.\right.\end{array}$$

(11)

Here, the MIMO EME model of the turbofan engine is established according to Equation (11). The modeling procedure, which is called the steady-state and dynamic two-step method, is divided into two steps:

• Based on the steady-state results obtained by the simulation of the turbofan engine NCL model, the steady-state EM results of the engine shown in Equation (8) are identified.
• In the NCL model simulation process, the input variable signal is set to the staircase signal. The EME model coefficients in Equation (11) are identified by simulation results of the NCL model and the EM model.

Firstly, it is significant to obtain steady-state EM of turbofan engine for modeling the EME model that shown in Equation (11). Normally, the EM is identified by polynomial. Considering the accuracy of the identification, the engine EM is identified with the m-order polynomial as shown below.

$$\begin{array}{c}\hfill {\xi }_e=\sum _{i=0}^m\sum _{j=0}^{m-i}{k}_{nij}{\alpha }_1^i{\alpha }_2^j,\end{array}$$

(12)

where ${\xi }_e$ is the EM of the turbofan engine, which is distinguished by the subscript n.

Table 1. Subscript corresponding to different EM.

${\mathit{\xi }}_{\mathit{e}}$	n
$n_{HCe}$	h
$n_{LCe}$	l
$W_{fe}$	f
$A_{ne}$	s
$EPR$	p

shows the subscript n corresponding to different EM ${\xi }_e$.

When $n_{HC}$ and $A_n$ are selected as scheduling variables $\alpha$, the EM of the turbofan engine are $n_{LCe}$, $EP{R}_e$, and $W_{fe}$, which have the structure of Equation (12). The one-to-one correspondence between the input variables $\left(A_n,W_f\right)$ and the performance parameters $\left(n_{HC},n_{LC},EPR\right)$ of the turbofan engine is shown as the black circular mark in Figure 2, Figure 3 and Figure 4. These data are obtained by simulating the turbofan engine NCL model carried out in Section 3.1. Therefore, $n_{HC}$ and $A_n$ are taken as independent variables to extract all steady-state results from simulations, and the corresponding dependent variables $n_{LCe}$, $EP{R}_e$, and $W_{fe}$ are obtained. By using the least square method and selecting the 4-order polynomial $\left(m=4\right)$, the identification results of EM $\left(\alpha =[n_{HC},A_n]\right)$ for the turbofan engine are shown in

Table 2. The coefficients of the EM $\left(\alpha =[n_{HC},A_n]\right)$.

${\mathit{n}}_{\mathit{LCe}}$		${\mathit{W}}_{\mathit{fe}}$		${\mathit{EPR}}_{\mathit{e}}$
Coef 1	Result 1	Coef 2	Result 2	Coef 3	Result 3
$k_{l00}$	−0.0031	$k_{f00}$	−0.0278	$k_{p00}$	0.1875
$k_{l01}$	0.0952	$k_{f01}$	0.1859	$k_{p01}$	0.0355
$k_{l02}$	0.0029	$k_{f02}$	−0.4295	$k_{p02}$	−0.5316
$k_{l03}$	0.1822	$k_{f03}$	0.4005	$k_{p03}$	0.4544
$k_{l04}$	−0.0662	$k_{f04}$	−0.1310	$k_{p04}$	−0.1465
$k_{l10}$	0.1586	$k_{f10}$	1.0829	$k_{p10}$	0.8436
$k_{l11}$	1.1210	$k_{f11}$	0.4523	$k_{p11}$	0.2797
$k_{l12}$	−0.2957	$k_{f12}$	−0.2919	$k_{p12}$	−0.0796
$k_{l13}$	−0.0700	$k_{f13}$	−0.1283	$k_{p13}$	−0.2328
$k_{l20}$	2.0277	$k_{f20}$	1.6503	$k_{p20}$	1.6111
$k_{l21}$	−2.0877	$k_{f21}$	−2.8265	$k_{p21}$	−3.1073
$k_{l22}$	0.2379	$k_{f22}$	0.8793	$k_{p22}$	0.9245
$k_{l30}$	−1.9790	$k_{f30}$	−1.5050	$k_{p30}$	−1.6747
$k_{l31}$	1.1385	$k_{f31}$	1.4396	$k_{p31}$	1.7264
$k_{l40}$	0.5497	$k_{f40}$	0.2737	$k_{p40}$	0.2624

. Meanwhile, the comparison steady-state results and error analysis are shown in Figure 5, Figure 6 and Figure 7.

In Figure 5a, Figure 6a, and Figure 7a, the identification results of EM $\left(\alpha =[n_{HC},A_n]\right)$ are presented. Among them, the black points are the steady-state results of the NCL model, and the red asterisks are the fitting result. Meanwhile, the absolute errors of the EM with the NCL model’s steady-state results are given in sub-figure (b). Firstly, Figure 5 indicates that the absolute error range of $n_{LCe}$ is about −0.015–0.017. In the region composed of $n_{HC}$ (0–0.4) and $A_n$ (0–0.5), the fluctuation of the $n_{LCe}$ is caused by the strong nonlinearity of the engine, and the error changes from −0.0084 ($n_{HC}$ = 0.1034, $A_n$ = 0) to 0.0168 ($n_{HC}$ = 0.1885, $A_n$ = 0.1111) and then to −0.0132 ($n_{HC}$ = 0.3761, $A_n$ = 0.1111). The same phenomenon also occurs when $n_{HC}$ and $A_n$ become larger, but the fluctuation amplitude is smaller. Secondly, the absolute error range of $EP{R}_e$ is about −0.022–0.028, and the maximum value is 0.028 ($n_{HC}$ = 1, $A_n$ = 1) in Figure 6. There appears the same strong nonlinear region, in which the error of $EP{R}_e$ changes by a wide margin from −0.0219 ($n_{HC}$ = 0.009, $A_n$ = 0) to 0.0196 ($n_{HC}$ = 0.1885, $A_n$ = 0.1111) and then to −0.0207 ($n_{HC}$ = 0.3761, $A_n$ = 0.1111). Thirdly, Figure 7 shows the absolute error of $W_{fe}$. Its range is about −0.02–0.024, and the maximum one is 0.0246 ($n_{HC}$ = 1, $A_n$ = 1). Besides, the absolute error of $W_{fe}$ changes from −0.0176 ($n_{HC}$ = 0.009, $A_n$=0) to 0.0180 ($n_{HC}$ = 0.1885, $A_n$ = 0.1111) and then to −0.0188 ($n_{HC}$ = 0.3814, $A_n$ = 0.2222) in the same strong nonlinear region. Overall, the identification result of the $n_{LCe}$ is better than the $EP{R}_e$, and $W_{fe}$, which is less than 0.017.

The fourth-order polynomial used in this section meets requirement of the accuracy. Usually, when a low order polynomial is used to fit the EM of the system with strong nonlinearity, the steady-state accuracy may not be satisfied. It is considered that increasing order of polynomial is a method to improve the fitting accuracy of the EM. However, the EM with a high-order polynomial also brings problems, such as overfitting and complex expression forms, when the result of the EM is close to the data from the NCL model. It is not conducive to the application of the EME model. In addition, the future research can also adopt switching method and neural network to solve the accuracy problem of the EM.

Secondly, the dynamic and static two-step identification method is applied to fit the EME model. After obtaining the steady-state EM, the coefficient matrix $A\left(\alpha \right)$, $B\left(\alpha \right)$, $C\left(\alpha \right)$, and $D\left(\alpha \right)$ in Equation (11) is identified by the dynamic simulation results of the turbofan engine NCL model. In the simulation of the NCL model, the input variables $W_f$ and $A_n$ that contain two paths considering the sufficient excitation conditions for the identification are shown in Figure 8. As described above, the identification results of the EM require a large number of steady-state data to ensure the accuracy. However, this is not the case for the dynamic parameters of the EME model. If a large number of dynamic processes are adopted or the variation interval of input parameters is excessively reduced in the identification, the results of EME model may be unstable. Therefore, some simulation paths are provided in Figure 8a to obtain the data needed for the dynamic coefficients of the EME model. Path 1 and Path 2 are the diagonal lines of the feasible region of the engine input, and they contain a certain degree of dynamic characteristics. However, neither Path 1 nor Path 2 alone can identify stable dynamic parameters of the EME model. This situation is alleviated after combining Path 1 and Path 2 (Path 1 & 2). At the same time, through the well-designed dynamic parameter identification Path New, better identification results are obtained.

Table 3. Accuracy results of EME model under different dynamic paths.

Path	$\mathit{\alpha }$	Output	$\mathit{SSE}$	$\mathit{RMSE}$	${\mathit{R}}^2$
1&2	$\left(n_{HC},A_n\right)$	$n_{HC}$	0.9508	0.0103	0.9979
		$n_{LC}$	0.1539	0.0041	0.9996
		${\pi }_C$	0.4772	0.0073	0.9990
New	$\left(n_{HC},A_n\right)$	$n_{HC}$	0.2701	0.0035	0.9994
		$n_{LC}$	0.1623	0.0042	0.9996
		${\pi }_C$	0.1615	0.0042	0.9997
1&2 with noise	$\left(n_{HC},A_n\right)$	$n_{HC}$	0.9204	0.0101	0.9980
		$n_{LC}$	0.1912	0.0046	0.9995
		${\pi }_C$	0.5516	0.0078	0.9988
1&2	$\left(W_f,A_n\right)$	$n_{HC}$	0.3830	0.0065	0.9992
		$n_{LC}$	0.2505	0.0053	0.9994
		${\pi }_C$	0.1983	0.0047	0.9996

shows the accuracy results of EME model under different dynamic paths in this paper. Analysis indicators include the sum of squares due to error ($SSE$), root mean squared error ($RMSE$), and coefficient of determination ($R^2$). The closer the $SSE$ and $RMSE$ are to 0, the better the model selection and fitting, and the more successful the data prediction. The closer $R^2$ is to 1, the better the model fits the data. The results show that the path for dynamic parameters fitting will affect the accuracy of the EME model. Path New has higher accuracy than Path 1 & 2. Facing different engines, a well-designed dynamic path can improve the accuracy of the EME model. The difficulty is that the dynamic path selection is random, and the optimal path for different engines is also difficult to obtain. On the contrary, the combination of Path 1 and Path 2 is the general operating range of turbofan engines. From the results of Table 3, the EME model identified by Path 1 & 2 has the satisfactory accuracy. Therefore, the study adopts Path 1 & 2 to carry out dynamic parameter identification of EME model, and carries out follow-up work based on this path.

Figure 8b shows the input variables changing along Path 1 (0–90 s) and Path 2 (90–180 s). Thus, the scheduling variables $\left(\alpha \right)$, derivatives ($\dot{X}$), and deviations ($\mathsf{\Delta }X=X-X\left(\alpha \right)$, $\mathsf{\Delta }U=U-U\left(\alpha \right)$) in the EME model of the turbofan engine can be obtained. That is, $\mathsf{\Delta }{W}_f$, $\mathsf{\Delta }{n}_{LC}$, ${\dot{n}}_{HC}$, and ${\dot{n}}_{LC}$ in Equation (11) can be calculated by the NCL and EM simulations. The coefficient matrix of the EME model is identified as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{a}_{12}\left(\alpha \right)=k_{a21}{\alpha }_1+k_{a22}{\alpha }_2+k_{a23}\hfill \\ a_{22}\left(\alpha \right)=k_{a41}{\alpha }_1+k_{a42}{\alpha }_2+k_{a43}\hfill \\ b_{11}\left(\alpha \right)=k_{b11}{\alpha }_1+k_{b12}{\alpha }_2+k_{b13}\hfill \\ b_{21}\left(\alpha \right)=k_{b31}{\alpha }_1+k_{b32}{\alpha }_2+k_{b33}\hfill \\ c_{12}\left(\alpha \right)=k_{c21}{\alpha }_1+k_{c22}{\alpha }_2+k_{c23}\hfill \\ d_{11}\left(\alpha \right)=k_{d11}{\alpha }_1+k_{d12}{\alpha }_2+k_{d13}\hfill \\ \alpha =\left[\begin{array}{cc}{\alpha }_1\hfill & {\alpha }_2\hfill \end{array}\right]=\left[n_{HC},A_n\right]\hfill \end{array}.\right.\end{array}$$

(13)

The identification results of Equation (13) are shown in

Table 4. The coefficients of the EME $\left(\alpha =[n_{HC},A_n]\right)$ model.

${\dot{\mathit{n}}}_{\mathit{HC}}$		${\dot{\mathit{n}}}_{\mathit{LC}}$		$\mathit{EPR}$
Coef 1	Result 1	Coef 2	Result 2	Coef 3	Result 3
$k_{a21}$	−0.2362	$k_{a41}$	−1.4657	$k_{c21}$	−0.0875
$k_{a22}$	1.5140	$k_{a42}$	0.2270	$k_{c22}$	−0.0401
$k_{a23}$	−2.1259	$k_{a43}$	−2.0621	$k_{c23}$	−0.9159
$k_{b11}$	0.0148	$k_{b31}$	−0.1033	$k_{d11}$	0.0749
$k_{b12}$	0.0023	$k_{b32}$	−0.0175	$k_{d12}$	−0.0227
$k_{b13}$	1.7453	$k_{b33}$	1.3592	$k_{d13}$	0.3058

.

So far, all the parameters of the EME model have been identified. Naturally, it is necessary to validate the accuracy of the EME model. During the verification, the input variables $W_f$ and $A_n$ are still shown in Figure 8. In addition, the comparison and error analysis of the turbofan engine’s NCL model and EME model are shown in Figure 9. The results given by Figure 9 show that the $n_{LC}$ and $EPR$ of the EME ($\alpha =\left[n_{HC},A_n\right]$) have high steady-state accuracy, while the $n_{HC}$ has low steady-state accuracy. Obviously, the steady-state accuracy of $n_{HC}$ and $n_{LC}$ is determined by the identification results of $n_{HCe}$ and $n_{LCe}$, respectively. Meanwhile, the simulations also indicate that the $n_{LC}$ and $EPR$ meet high dynamic accuracy, in which the dynamic time of the EME model is similar to that of the NCL model. In particular, there appears a dynamic error with absolute value around 0.02 for the simulations of $n_{LC}$ at 10 s and 70 s in Figure 9b. In Figure 9c, the results of $EPR$ show a dynamic error with absolute value around 0.05. The large dynamic error of the $EPR$ may be caused by the negative regulation characteristics. However, the $n_{HC}$ shown in Figure 9a is not satisfactory. In the simulation, large overshoot occurs in the EME model at 20 s, 60 s, and 70 s, and the overshoot at 70 s is close to 0.02.

3.3. Effects of Scheduling Variables on Model Accuracy

In Section 3.2, there are some errors in the steady-state and dynamic accuracy of EME ($\alpha =\left[n_{HC},A_n\right]$) model, especially in $n_{HC}$. Thankfully, the previous studies have shown that the optimal identification form of the EME model always exists in the nonlinear plants with analytic representation. Although MIMO turbofan engine with strong nonlinearity cannot be expressed analytically, the accuracy of the EME model can also be improved by selecting different combinations of the scheduling variables.

$W_f$ and $A_n$ are selected as $\alpha$ to complete the EME model of the turbofan engine in this section. Thus, the form of $\alpha$ can be rewritten as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\alpha }_1=W_f=W_f\left(\alpha \right)\hfill \\ {\alpha }_2=A_n=A_n\left(\alpha \right)\hfill \\ \alpha =\left[{\alpha }_1,{\alpha }_2\right]\hfill \end{array}.\right.\end{array}$$

(14)

According to Equation (10), $\mathsf{\Delta }{W}_f=W_f-W_{fe}\left(\alpha \right)=0$, and $\mathsf{\Delta }{A}_n=A_n-A_{ne}\left(\alpha \right)=0$. Thus, the EME model of the turbofan engine shown in Equation (9) can be replaced by:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\dot{n}}_{HC}=a_{11}\left(\alpha \right)\mathsf{\Delta }{n}_{HC}+a_{12}\left(\alpha \right)\mathsf{\Delta }{n}_{LC}\hfill \\ {\dot{n}}_{LC}=a_{21}\left(\alpha \right)\mathsf{\Delta }{n}_{HC}+a_{22}\left(\alpha \right)\mathsf{\Delta }{n}_{LC}\hfill \\ EPR=EP{R}_e\left(\alpha \right)+c_{11}\left(\alpha \right)\mathsf{\Delta }{n}_{HC}+c_{12}\left(\alpha \right)\mathsf{\Delta }{n}_{LC}\hfill \\ n_{HC}=n_{HCe}\left(\alpha \right)+\mathsf{\Delta }{n}_{HC}\hfill \\ n_{LC}=n_{LCe}\left(\alpha \right)+\mathsf{\Delta }{n}_{LC}\hfill \\ \alpha =\left[{\alpha }_1,{\alpha }_2\right]=\left[W_f,A_n\right]\hfill \end{array}.\right.\end{array}$$

(15)

Firstly, we identify the EM of the turbofan engine. When $W_f$ and $A_n$ are selected as scheduling variables $\alpha$, the EMs of the turbofan engine are $n_{HCe}$, $n_{LCe}$, and $EP{R}_e$ which meet the form of Equation (12). Based on the method given in Section 3.2, the identification results of steady-state EM are shown in

Table 5. The coefficients of the EM $\left(\alpha =[W_f,A_n]\right)$.

${\mathit{n}}_{\mathit{HCe}}$		${\mathit{n}}_{\mathit{LCe}}$		${\mathit{EPR}}_{\mathit{e}}$
Coef 1	Result 1	Coef 2	Result 2	Coef 3	Result 3
$k_{h00}$	0.0047	$k_{l00}$	0.0044	$k_{p00}$	0.2207
$k_{h01}$	0.0261	$k_{l01}$	0.0832	$k_{p01}$	−0.2133
$k_{h02}$	−0.0989	$k_{l02}$	−0.0001	$k_{p02}$	0.0389
$k_{h03}$	0.0896	$k_{l03}$	0.1935	$k_{p03}$	−0.0500
$k_{h04}$	−0.0207	$k_{l04}$	−0.0678	$k_{p04}$	0.0054
$k_{h10}$	1.0655	$k_{l10}$	0.2134	$k_{p10}$	0.6443
$k_{h11}$	−1.0625	$k_{l11}$	0.2511	$k_{p11}$	0.5970
$k_{h12}$	1.2794	$k_{l12}$	0.6322	$k_{p12}$	−0.7356
$k_{h13}$	−0.3306	$k_{l13}$	−0.2856	$k_{p13}$	0.1915
$k_{h20}$	−1.1696	$k_{p20}$	0.8536	$k_{f20}$	0.5171
$k_{h21}$	1.9950	$k_{p21}$	0.0144	$k_{f21}$	−1.2565
$k_{h22}$	−0.7915	$k_{p22}$	−0.4243	$k_{f22}$	0.5088
$k_{h30}$	1.0608	$k_{p30}$	−0.9268	$k_{f30}$	−0.5161
$k_{h31}$	−0.7422	$k_{p31}$	0.0930	$k_{f31}$	0.4477
$k_{h40}$	−0.3170	$k_{p40}$	0.3604	$k_{f40}$	0.1315

. Meanwhile, the steady-state results and the absolute error analysis of EM with NCL model are shown in Figure 10, Figure 11 and Figure 12. The identification results of EM ($\alpha =\left[n_{HC},A_n\right]$) are presented in Figure 10a, Figure 11a, and Figure 12a. Among them, the black points are the steady-state results of the NCL model in Section 3.1, and the red asterisks are fitting results. The absolute error of the EM with the NCL model’s steady-state results is given in Figure 10b, Figure 11b, and Figure 12b.

Clearly, the operating range of the $n_{HCe}$ is shown in Figure 10, and the maximum absolute error is −0.0120 $\left(W_f=1,A_n=1\right)$. However, the strong nonlinear region of the EM in EME ($\alpha =\left[n_{HC},A_n\right]$) is not significant in EME ($\alpha =\left[W_f,A_n\right]$). Figure 11 shows that the absolute error range of $n_{LCe}$ is about −0.006–0.010. The extremum of the absolute errors are −0.0058 ($W_f=0,A_n=0.3333$) and 0.0098 ($W_f=0.1111,A_n=0.4444$) in Figure 11b. In addition, Figure 12b indicates that the absolute error range of $EP{R}_e$ is about −0.0060–0.0060. The extremum of the absolute errors are −0.0061 ($W_f=0.3333,A_n=0.5556$) and 0.0040 ($W_f=0.3333,A_n=0.2222$) in Figure 12b. In general, the accuracy of the EM in EME ($\alpha =\left[W_f,A_n\right]$) is higher than in EME ($\alpha =\left[n_{HC},A_n\right]$).

Secondly, the coefficient matrix $A\left(\alpha \right)$, $B\left(\alpha \right)$, $C\left(\alpha \right)$ and $D\left(\alpha \right)$ in Equation (11) is identified by the dynamic simulation results of the turbofan engine NCL model. Clearly, the specific fitting method is given in Section 3.2. That is, $\mathsf{\Delta }{n}_{HC}$, $\mathsf{\Delta }{n}_{LC}$, ${\dot{n}}_{HC}$, and ${\dot{n}}_{LC}$ in Equation (15) can be calculated by the results of the NCL model and the EM. The coefficients of the EME model can be identified as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{a}_{11}\left(\alpha \right)=k_{a11}{\alpha }_1+k_{a12}{\alpha }_2+k_{a13}\hfill \\ a_{12}\left(\alpha \right)=k_{a21}{\alpha }_1+k_{a22}{\alpha }_2+k_{a23}\hfill \\ a_{21}\left(\alpha \right)=k_{a31}{\alpha }_1+k_{a32}{\alpha }_2+k_{a33}\hfill \\ a_{22}\left(\alpha \right)=k_{a41}{\alpha }_1+k_{a42}{\alpha }_2+k_{a43}\hfill \\ c_{11}\left(\alpha \right)=k_{c11}{\alpha }_1+k_{c12}{\alpha }_2+k_{c13}\hfill \\ c_{12}\left(\alpha \right)=k_{c21}{\alpha }_1+k_{c22}{\alpha }_2+k_{c23}\hfill \\ \alpha =\left[{\alpha }_1,{\alpha }_2\right]=\left[W_f,A_n\right]\hfill \end{array}.\right.\end{array}$$

(16)

Therefore, the identification results of Equation (16) are shown in

Table 6. The coefficients of the EME $\left(\alpha =[W_f,A_n]\right)$ model.

${\dot{\mathit{n}}}_{\mathit{HC}}$		${\dot{\mathit{n}}}_{\mathit{LC}}$		$\mathit{EPR}$
Coef 1	Result 1	Coef 2	Result 2	Coef 3	Result 3
$k_{a11}$	0.5015	$k_{a31}$	3.7503	$k_{c11}$	−0.5078
$k_{a12}$	−0.0845	$k_{a32}$	1.0956	$k_{c12}$	−0.1344
$k_{a13}$	−2.2704	$k_{a33}$	−1.4782	$k_{c13}$	0.3204
$k_{a21}$	0.6230	$k_{a41}$	−4.2386	$k_{c21}$	0.0782
$k_{a22}$	2.2073	$k_{a42}$	0.9056	$k_{c22}$	−0.7508
$k_{a23}$	−2.6501	$k_{a43}$	−2.0944	$k_{c23}$	0.9199

.

Because the validation is an important part of modeling, so the validation of the EME model is carried out as follows. The input variables $W_f$ and $A_n$ vary according to the rule shown in Figure 8. Through the simulation, the comparison and error analysis of the turbofan engine’s NCL model and EME model are shown in Figure 13. The results demonstrate the accuracy of the EME model along the design Path 1 & 2.

Figure 13 indicates that $n_{HC}$, $n_{LC}$, and $EPR$ in EME ($\alpha =\left[W_f,A_n\right]$) all have high steady-state accuracy, which is determined by the identification results of the EM. Simultaneously, the results also have high dynamic accuracy in EME ($\alpha =\left[W_f,A_n\right]$). There are large dynamic errors for $n_{HC}$, $n_{LC}$, and $EPR$ at 60s, and the values are 0.025, 0.0275, and −0.055, respectively. The EME ($\alpha =\left[W_f,A_n\right]$) model can reflect the overshoot characteristics of $n_{HC}$ in the range of 0–0.035. Besides, the overshoot is less than the EME ($\alpha =\left[n_{HC},A_n\right]$) at 60 s in Figure 13a. The accuracy of the EME model is affected by $\alpha$, which is confirmed in the simulation. The analysis indicators in Table 3 shows that the EME ($\alpha =\left[W_f,A_n\right]$) model has higher steady-state accuracy than the EME ($\alpha =\left[n_{HC},A_n\right]$) model, and the identification result of the $n_{HC}$ is more satisfactory. However, direct measurement of inputs, such as real engine fuel, will limit the feasibility of inputs as scheduling variables for the EME model.

3.4. Effects of Noise on Model Accuracy

The EME model is a turbofan engine mathematical model based on data identification. These data can be obtained by dynamic simulation of the NCL model, on the one hand, and can be obtained by a large number of real engine tests, on the other hand. However, the cost of engine testing is too high to be used as a data source for the study of EME modeling methods. In fact, there is a big difference between NCL model simulation data and engine test data, among which noise is an important influencing factor. To analyze the influence of noise on EME model identification, we add random noise to the input $W_f$ and $A_n$ of the dynamic Path 1 & 2 in Section 3.2, and complete the EME model identification according to the same process. The output results of the NCL model and the EME model are shown in Table 3.

First, it can be seen from Figure 14 that the input noise has a greater impact on the engine pressure ratio $EPR$ than two rotors $n_{HC}$ and $n_{LC}$. Secondly, under the effect of input noise, the EME model still has a strong model fitting ability. From the perspective of steady-state results, the EME model can ensure stability; meanwhile, the steady-state accuracy of $n_{HC}$ and $n_{LC}$ is kept about within 2%. Although the pressure ratio $EPR$ error is relatively large, the EME model ensures the steady-state accuracy of the engine and shows a certain filtering ability.

4. Validation and Stability Analysis of the EME Model

4.1. Validation at Dynamic Off-Design Conditions

This paper proposes that the dynamic coefficients of the EME model are identified by Path 1 & 2 shown in Figure 15a. Despite differences in accuracy, the identification results of the EME model shown in Section 3.2 and Section 3.3 are satisfactory in control-oriented modeling. However, verifying the accuracy of the EME model on the design path alone is insufficient. Therefore, this section considers the validation of the EME ($\alpha =\left[n_{HC},A_n\right]$) model and the EME ($\alpha =\left[W_f,A_n\right]$) model at dynamic off-design conditions. Path 3, which is different from Path 1 and Path 2, is given in Figure 15a. Figure 15b shows the corresponding input variables for the validation.

The input variables given in Figure 15b are used in the NCL, EME ($\alpha =\left[n_{HC},A_n\right]$), and EME ($\alpha =\left[W_f,A_n\right]$) model of the turbofan engine for simulation. The output variables $n_{HC}$, $n_{LC}$ and $EPR$ of the turbofan engine are shown in Figure 16. At the same off-design conditions, the validation results of different turbofan engine models, and the absolute errors between the EME ($\alpha =\left[n_{HC},A_n\right]$) model to the NCL model and the EME ($\alpha =\left[W_f,A_n\right]$) model to the NCL model are shown in Figure 16.

Firstly, the steady-state $n_{HC}$, $n_{LC}$ and $EPR$ of the EME model with different scheduling variables have different accuracy. In particular, Figure 16a shows that there is a difference between the EME ($\alpha =\left[n_{HC},A_n\right]$) and EME ($\alpha =\left[W_f,A_n\right]$) model in the steady-state accuracy of $n_{HC}$. The EME ($\alpha =\left[W_f,A_n\right]$) model has higher steady-state accuracy at the off-design conditions given in Path 3, which is consistent with the results obtained from Path 1 and 2. In addition, Figure 16b indicates that the steady-state accuracy of $n_{LC}$ at off-design conditions is almost the same, except that the steady-state accuracy of the EME ($\alpha =\left[n_{HC},A_n\right]$) is lower than the EME ($\alpha =\left[W_f,A_n\right]$) at 70–80 s. The steady-state accuracy of $EPR$ at off-design conditions is similar to that shown in Figure 16c, except that the EME ($\alpha =\left[n_{HC},A_n\right]$) has lower steady-state accuracy than the EME ($\alpha =\left[W_f,A_n\right]$) at 40–50 s. Secondly, Figure 16 shows that the dynamic time of the EME ($\alpha =\left[n_{HC},A_n\right]$) and EME ($\alpha =\left[W_f,A_n\right]$) is almost the same, which means that the two EME models can well meet the dynamic process of the turbofan engine. Meanwhile, Figure 16a indicates that there is an overshoot of $n_{HC}$ at 30 s, which is caused by the strong nonlinearity of the engine.

In short, the EME model has satisfactory accuracy at both design and off-design conditions. The model meets the demand of real-time simulation. These are significant for the control-oriented model of the turbofan engine. Although it is difficult to find the EME model with the smallest error for the strongly nonlinear system including turbofan engines that cannot be expressed analytically, we can identify the EME model that meet the requirements of the accuracy by selecting different combinations of the scheduling variables.

4.2. Linearization and Stability Analysis

As a control-oriented model of the MIMO turbofan engine, it is significant for the EME model to meet the requirements of the control system design. At present, the most common and well-developed control methods are linear control, while the nonlinear control methods have difficulties in solving and practical applications. Moreover, many nonlinear control methods involve the linearization of the complex nonlinear plants. Therefore, a wide-range linear model with high accuracy can greatly reduce the difficulty of a control system design for nonlinear plants. Nevertheless, the EME model is a nonlinear parameter-varying system, and some coefficients of the linearization for the EME model are missing in $A\left(\alpha \right)$, $B\left(\alpha \right)$, $C\left(\alpha \right)$, and $D\left(\alpha \right)$ due to the selection of the scheduling variables. $a_{11}\left(\alpha \right)$, $a_{21}\left(\alpha \right)$, $b_{12}\left(\alpha \right)$, $b_{22}\left(\alpha \right)$, $c_{11}\left(\alpha \right)$, and $d_{12}\left(\alpha \right)$ are missing in Equation (11), whereas $b_{11}\left(\alpha \right)$, $b_{12}\left(\alpha \right)$, $b_{21}\left(\alpha \right)$, $b_{22}\left(\alpha \right)$, $d_{11}\left(\alpha \right)$, and $d_{12}\left(\alpha \right)$ are missing in Equation (15). Thus, it is important to obtain a complete linearization result of the EME model. Fortunately, Property 2 of the EME model introduced in Section 2.2 can be used to complete the linearization results, which contains the constraints of the EME model. For Equation (11), the sufficient and necessary conditions Equation (6) given in Property 2 can be expanded to:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{a}_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial n_{HCe}}+a_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial n_{HCe}}+b_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial n_{HCe}}+b_{12}\left(\alpha \right)\frac{\partial A_{ne}\left(\alpha \right)}{\partial n_{HCe}}=0\hfill \\ a_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial A_{ne}}+a_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial A_{ne}}+b_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial A_{ne}}+b_{12}\left(\alpha \right)\frac{\partial A_{ne}\left(\alpha \right)}{\partial A_{ne}}=0\hfill \\ a_{21}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial n_{HCe}}+a_{22}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial n_{HCe}}+b_{21}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial n_{HCe}}+b_{22}\left(\alpha \right)\frac{\partial A_{ne}\left(\alpha \right)}{\partial n_{HCe}}=0\hfill \\ a_{21}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial A_{ne}}+a_{22}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial A_{ne}}+b_{21}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial A_{ne}}+b_{22}\left(\alpha \right)\frac{\partial A_{ne}\left(\alpha \right)}{\partial A_{ne}}=0\hfill \\ c_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial n_{HCe}}+c_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial n_{HCe}}+d_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial n_{HCe}}+d_{12}\left(\alpha \right)\frac{\partial A_{ne}\left(\alpha \right)}{\partial n_{HCe}}=\frac{\partial EP{R}_e\left(\alpha \right)}{\partial n_{HCe}}\hfill \\ c_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial A_{ne}}+c_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial A_{ne}}+d_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial A_{ne}}+d_{12}\left(\alpha \right)\frac{\partial A_{ne}\left(\alpha \right)}{\partial A_{ne}}=\frac{\partial EP{R}_e\left(\alpha \right)}{\partial n_{HCe}}\hfill \\ \alpha =\left[n_{HCe},A_{ne}\right]=\left[n_{HC},A_n\right]\hfill \end{array}.\right.\end{array}$$

(17)

Considering the orthogonal properties of the scheduling variables, we have:

$$\begin{array}{c}\hfill \frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial n_{HCe}}=1,\frac{\partial A_{ne}\left(\alpha \right)}{\partial n_{HCe}}=0,\frac{\partial A_{ne}\left(\alpha \right)}{\partial A_{ne}}=1,\frac{\partial A_{ne}\left(\alpha \right)}{\partial n_{HCe}}=0.\end{array}$$

(18)

Substituting Equation (18) into Equation (17), the missing coefficients in Equation (11) can be completed as:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{a}_{11}\left(\alpha \right)=-a_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial n_{HCe}}-b_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial n_{HCe}}\\ b_{12}\left(\alpha \right)=-a_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e\left(\alpha \right)}}{\partial A_{ne}}-b_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial A_{ne}}\\ a_{21}\left(\alpha \right)=-a_{22}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial n_{HCe}}-b_{21}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial n_{HCe}}\\ b_{22}\left(\alpha \right)=-a_{22}\left(\alpha \right)\frac{\partial n_{LCe}\left(\alpha \right)}{\partial A_{ne}}-b_{21}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial A_{ne}}\\ c_{11}\left(\alpha \right)=\frac{\partial EP{R}_e\left(\alpha \right)}{\partial n_{HCe}}-c_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial n_{HCe}}-d_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial n_{HCe}}\\ d_{12}\left(\alpha \right)=\frac{\partial EP{R}_e\left(\alpha \right)}{\partial A_{ne}}-c_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial A_{ne}}-d_{11}\left(\alpha \right)\frac{\partial W_{fe}\left(\alpha \right)}{\partial A_{ne}}\\ \alpha =\left[n_{HCe},A_{ne}\right]=\left[n_{HC},A_n\right]\end{array}.\right.\end{array}$$

(19)

Similarly, the missing coefficients in Equation (15) can be completed as:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{b}_{11}\left(\alpha \right)=-a_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial W_{fe}}-a_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial W_{fe}}\\ b_{12}\left(\alpha \right)=-a_{11}\left(\alpha \right)\frac{\partial n_{HCe}\left(\alpha \right)}{\partial A_{ne}}-a_{12}\left(\alpha \right)\frac{\partial n_{LCe}\left(\alpha \right)}{\partial A_{ne}}\\ b_{21}\left(\alpha \right)=-a_{21}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial W_{fe}}-a_{22}\left(\alpha \right)\frac{\partial n_{L{C}_{\theta }}\left(\alpha \right)}{\partial W_{fe}}\\ b_{22}\left(\alpha \right)=-a_{21}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial A_{ne}}-a_{22}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial A_{ne}}\\ d_{11}\left(\alpha \right)=\frac{\partial EP{R}_e\left(\alpha \right)}{\partial W_{fe}}-c_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial W_{fe}}-c_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial W_{fe}}\\ d_{12}\left(\alpha \right)=\frac{\partial EP{R}_e\left(\alpha \right)}{\partial n_{HCe}}-c_{11}\left(\alpha \right)\frac{\partial n_{H{C}_e}\left(\alpha \right)}{\partial A_{ne}}-c_{12}\left(\alpha \right)\frac{\partial n_{L{C}_e}\left(\alpha \right)}{\partial A_{ne}}\\ \alpha =\left[W_{fe},A_{ne}\right]=\left[W_f,A_n\right]\end{array}.\right.\end{array}$$

(20)

Obviously, after the EME ($\alpha =\left[n_{HC},A_n\right]$) model is complemented by Equation (19) and the EME ($\alpha =\left[W_f,A_n\right]$) model is complemented by Equation (20), the linear model of the MIMO turbofan engine at any steady-state operating point can be obtained. As the turbofan engine is a stable system which can be seen from the simulation in Section 3.1, the EME model should be stable, and it is verified in the time domain simulation shown in Figure 9 and Figure 13. However, the analysis from the time domain is not rigorous enough. Stability analysis is necessary for the linearized results of the EME model.

The sufficient and necessary condition for the stability of linear systems is that all the roots of the characteristic equations of closed-loop systems have negative real parts. This condition means all the $A\left(\alpha \right)$ matrix from the linearization results of the EME model with the state space form have negative real parts. Thus, the linearization of the EME ($\alpha =\left[n_{HC},A_n\right]$) model and the EME ($\alpha =\left[W_f,A_n\right]$) model are completed at the 100 steady-state points given in Figure 2, Figure 3 and Figure 4, and the real part of the eigenvalues of $A\left(\alpha \right)$ matrix is shown in Figure 17 and Figure 18.

Since all the eigenvalues of $A\left(\alpha \right)$ have negative real parts, the linearization results of the EME models ($\alpha =\left[n_{HC},A_n\right],\alpha =\left[W_f,A_n\right]$) are stable within the range of the scheduling variables. So far, the accuracy and stability of the EME model are proved by the time domain simulation and linearization results, which indicates that the EME model can be used as real-time simulation model for control system design.

5. Conclusions

This paper investigates a control-oriented model of nonlinear MIMO turbofan engines using the equilibrium manifold method. To obtain the EME model of turbofan engine with input variables $W_f$ and $A_n$ and output variables $n_{HC}$, $n_{LC}$, and $EPR$, extending the equilibrium manifold modeling method from the SISO system to the MIMO system is studied. The main conclusions of this article are:

• The two EME models that contain the input variables ($W_f,A_n$) and the output variables ($n_{HC},n_{LC},EPR$) are established. The results of the steady state EM ($\alpha =\left[n_{HC},A_n\right]$) show that the absolute error range of $n_{LCe}$ (−0.015–0.017) is better than $EP{R}_e$ (−0.022–0.028) and $W_{fe}$ (−0.02–0.024). The absolute error range of EM ($\alpha =\left[W_f,A_n\right]$) are about −0.006–0.010 in $n_{LCe}$, −0.0060–0.0060 in $EP{R}_e$, and −0.012–0.010 in $n_{HCe}$. Obviously, the accuracy of EME ($\alpha =\left[W_f,A_n\right]$) is higher than EME ($\alpha =\left[n_{HC},A_n\right]$). Because of the simple structure, the EME model also meets the needs of real-time simulation. Meanwhile, the EME model still shows high accuracy under the influence of noise.
• The EME model meets the accuracy requirement of the simulation and the stability in the entire range of operating conditions. The validity of the EME model is verified by simulation at off-design points. At the same time, the results also show that the identification of the dynamic coefficients for the EME model can be completed through only considering simple paths. Depending on the property of the EME model, the linearization of the model is completed. By calculating the eigenvalues, both EME models ($\alpha =\left[n_{HC},A_n\right]$,$\alpha =\left[W_f,A_n\right]$) have stable linearization result, thus ensuring that the EME model can be regarded as a control-oriented model.

Furthermore, the dynamic coefficients of the EME model have low dependence of the excitation condition in identification. It means that the identification of the dynamic coefficients needs few data, which is different from the general method of model identification. This finding provides a possibility for the EME model to solve the performance degradation of the turbofan engine in future.