# Direct Surge Margin Control for Aeroengines Based on Improved SVR Machine and LQR Method

Authors: Haibo Zhang; Fengyong Sun
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A novel scheme of high stability engine control (HISTEC) on the basis of an improved linear quadratic regulator (ILQR), called direct surge margin control, is derived for super-maneuver flights. Direct surge margin control, which is different from conventional control scheme, puts surge margin into the engine closed-loop system and takes surge margin as controlled variable directly. In this way, direct surge margin control can exploit potential performance of engine more effectively with a decrease of engine stability margin which usually happened in super-maneuver flights. For conquering the difficulty that aeroengine surge margin is undetectable, an approach based on improved support vector regression (SVR) machine is proposed to construct a surge margin prediction model. The surge margin modeling contains two parts: a baseline model under no inlet distortion states and the calculation for surge margin loss under supermaneuvering flight conditions. The previous one is developed using neural network method, the inputs of which are selected by a weighted feature selection algorithm. Considering the hysteresis between pilot input and angle of attack output, an online scrolling window least square support vector regression (LSSVR) method is employed to firstly estimate inlet distortion index and further compute surge margin loss via some empirical look-up tables.


---

## Body

## 1. Introduction

To begin with, we will provide a brief background on high stability control for aeroengines, which is an aerodynamic concept different from the stability concept in control theory and an unstable state means that the fan or compressor of an engine goes into a surge state. Super-maneuverability is one of the essential techniques of modern fighters. Yet, in a post-stall state, when engines work under the conditions of high angle of attack, the inlet distortion becomes severe, leading to unstable operations such as weakened engine performance, reduced steady operating surge margin, and even surge in serious conditions. For this issue, national aeronautics and space administration (NASA) initialed and led a famous high stability engine control (HISTEC) research project in 1993, in which an engine stability control in super-maneuver flights was highlighted, and the main idea is that the pressure ratio of an F-100 turbofan engine could be adaptively controlled to regulate cooperating working point or the surge margin of its fan and compressor when an F-15 ACTIVE (advanced control technology for integrated vehicles) aircraft entered into super-maneuver states. Thus, sufficient stability margin of engines can be ensured in super maneuver state [1–3]. Similarly, Wang et al. [4, 5] proposed an engine stability control law through the recovery of surge margin by compensating losses of the fan pressure ratio when a severe inlet distortion occurs. The literatures [6, 7] reported a direct surge margin control, in which a baseline and its loss for surge margin in severe distorted states need to compute, but the calculation for the baseline just employs a simpler definition that related to the fan pressure ratio; as a result it is impossible to reflect the influence in air flow. Although these above two schemes are able to ensure the engine stability, since the surge margin cannot be measured directly, it may be enlarged a lot in the maneuver process for the engine potential not to be fully exploited, resulting in much smaller engine thrust which cannot help the aircraft to accomplish the fast and right modulations in angular position.Around this interesting topic of aeroengine stability control, a lot of articles, which concerns active component control, have been reported [8–12]. In recent years, an approach of more realistic active stall/surge control has been proposed as an effective approach to realize high stability control [13–15]. This technology mainly takes advantage of the flow characteristics which is denoted as stall inception in high frequency phase of compressor and implements a specific mathematical model between pressure correlation and surge margin. Thus, when the engine enters into a severe distortion state, some emergency fuel control with high frequency can be applied to the engine system so that the temperature and pressure of the combustor will drop quickly; therefore the stall margin will be restored to the safe area. Essentially, this control approach is within the scope of limit protection control near the surge border, the defections of which are still found in a worse coordination with the main control loop, as well as an insufficient exploration in performance potentials.Obviously, those above methods have a comment prominent drawback which is that surge margin cannot be measured directly or be estimated accurately. Therefore, if the real-time accurate prediction of engine surge margin could be achieved, it can be used as the direct virtual parameter for accurate control, so that a fixed surge margin or a certain distance, between the cooperating work line and the surge border line in the compressor, can be always kept in super-maneuver flights. If this imagination could be realized, engines might sustain a more efficient and stable operation in super maneuver states. Near recently from a NASA report [16], it is right an expected solution to engine high stability control to apply some nonlinear methods like artificial mapping and compound kalman filter [17].As said above, the establishment of engine surge margin prediction model has always been a challenge in the field of surge control. Our research is carried out on a two-spool mixing exhaust turbofan engine with afterburning, the power of which lies in the same class of F-100 engine. First, through a novel feature selection algorithm, the most suitable measurable variables with the strongest correlation with surge margin are selected as inputs of a surge margin baseline model. Furthermore, based on a BP neural network, a surge margin baseline model is set up, which is able to predict the surge margin in undistorted states. In addition, an online scrolling window LSSVR model (OSW-LSSVR), which has a time series of the angle of attack and relative rudder angles as its inputs, is proposed and designed for the real-time prediction of inlet distortion index in distortion states. Due to a definite relationship between inlet distortion index and surge margin loss, the surge margin loss can be estimated using the OSW-LSSVR model. Ultimately, with the combination of the output of surge margin baseline model and the loss in distortion states, the surge margin under engine inlet distortion states can be accurately estimated. As for the design of high stability engine control, an improved linear quadratic regulation (ILQR) robust control is chosen here due to its good robustness for adapting a large variation in engine power states and envelope points [18].


## 2. Main Idea for Direct Surge Margin Control

### 2.1. A Novel Direct Surge Margin Control

An integrated aircraft and engine dynamic model depicted in Figure1 is introduced to implement some necessary validations [6, 7]. This comprehensive model is composed by a dynamic aircraft model and its flight controller, a component-level turbofan engine model and its controller which will be designed as the proposed controller, and a conventional one for comparisons. For the aircraft model, the dynamics in longitudinal plane is only employed for simplicity, which can simulate level flight, climbing, accelerating, descending and super maneuver tasks. For the engine model, A8 and Wfb are selected as the control parameter, whereas Smf and Nc are the relative controlled variables. Therefore, the engine command r includes the references value of Smf and Nc, and e is a deviation vector with respect to command signal. A data set about the relationship between α and DC60 that is a type of distortion index are originated from flight test data of a F/A-18A aircraft in the Flight Research Center in NASA Dryden [19]. As discussed above, if the surge margin for an engine can be predicted accurately, a direct engine surge margin control, or a high stability engine control, easily can be further realized. In our research, a novel scheme is proposed (see Figure 1), where k represents the current moment and St is a sequential engine state values. The main idea is as follows.(a)
Considering the hysteresis between pilot input and rudder output, the current manipulations will affect the outputs ind-step time delay (or τ seconds in continuous sense) and the prediction for α is very important in a super maneuver task in which α changes very violently [1–3] to have a great impact on the surge margin. Therefore, two parameters are needed to control this high stability control scheme. One is S~mf(k+d), that is a predicted surge margin S~mf at the d-step ahead time, to adapt the variation in Smf due to the changes in α. The other need to control is Nc which regulates the power of the engine to supply enough thrust F.(b)
So when a super maneuver of the fighter is carried out, how to estimateS~mf(k+d) is a core problem in this scheme. For calculating S~mf(k+d), a time series of angle of α and the elevator δz are taken as inputs of the angle of attack prediction model to predict the d-step ahead value α~(k+d). Necessary engine state values in St series are obtained through the relative sensors to mapping the baseline value of Smf. After that, the α~(k+d) and St sequence are gathered as the inputs of fan surge margin prediction model; consequently the predicted surge margin S~mf(k+d) would be gotten in real time. For showing clearly, an equivalent pilot input delay system is labeled with a dotted rectangle in Figure 1.Figure 1
Structure for direct engine surge margin control.Remark 1.
In our research, the simulation step is preset as 20 ms; sod=τ/0.02=0.5/0.02=25 is predefined here based on [1–3].


### 2.2. Necessity to Estimate Surge Margin

For a turbofan engine, a parameter or a concept is defined as surge margin to show a physical meaning of a relative distance from the cooperate working point to the surge border (see Figure2), and the following equation is generally used to represent the surge margin of a fan, and similar for a compressor [20]:
(1)Sm=πs/Wcor,s-πo/Wcor,oπo/Wcor,o=πsπo·Wcor,oWcor,s-1.Figure 2
Surge margin change in a fan characteristic map.When engine inlet distortions happen due to high attack angle states, the surge border of the engine will be shifted down or the operation line will be shifted up, leading to a surge margin loss. Take the example of the shifting-down of surge border, and we have:(2)ΔSm=Sm-S~m=πs/Wcor,s-πo/Wcor,oπo/Wcor,o-πl/Wcor,s-πo/Wcor,oπo/Wcor,o=πs/Wcor,s-π1/Wcor,lπo/Wcor,o,
where ΔSm shows surge margin loss and S~m represents the surge margin when entering distortion; the subscript 1 shows the lowered stable boundary value.As can be seen from the above equations,Sm and ΔSm are parameters that cannot be measured. Under distorted states, if the undistorted Sm and the surge margin loss ΔSm might be accurately predicted, the distorted surge margin S~m will be expressed as
(3)S~m=Sm-ΔSm.For calculatingΔSm, an indirect parameter, called pressure ratio loss of the cooperating work point, is usually utilized according to the [21, 22]. Due to the shifting down of the boundary line, the pressure ratio loss ΔPrs is expressed as
(4)ΔPrs=πs/Wcor,s-π1/Wcor,lπs/Wcor,s.With the combination of (1) and (2) and the definition of ΔPrs, the surge margin loss is easily expressed as
(5)ΔSm=πs/Wcor,s-π1/Wcor,lπo/Wcor,o=πs/Wcor,s-π1/Wcor,lπs/Wcor,sπs/Wcor,sπo/Wcor,o=ΔPrs(1+Sm).


### 2.3. Surge Margin Prediction Model

Based on the current and historical information of input parameters, dynamic solutions for the baseline valueSm and the loss value ΔSm can be acquired or estimated by a novel engine surge margin model for super-maneuver flight, which consists of a surge margin baseline model and surge margin loss model. The modeling process is shown in Figure 3 and described as follows.Figure 3
Engine surge margin modeling in distortion states.


#### 2.3.1. Surge Margin Baseline Model

In the range of super-maneuver flight envelope (height of 3~7 km, march number of 0.3~0.7), a surge margin baseline model based on a BP neural network, which is under no inlet distortion of the fighter, is trained to predict surge margin baseline value. Consider a lot of and redundant measurable parameters correlated with surge margin, a feature selection algorithm, for choosing adaptable engine sensors to estimate surge margin, is proposed by a weighted LSSVR method. With a contribution criterion of each parameter to surge margin, a group of most related and affordable parameters are selected as the model inputs.


#### 2.3.2. Surge Margin Modeling in Distortion States

On the basis of surge margin baseline model, surge margin modeling in distortion states can be implemented in the following steps.First, an online scrolling window LSSVR model (OSW-LSSVR) is proposed to predict thed-step ahead value α~(k+d). The OSW-LSSVR model is capable of updating the prediction model based on the time series of angle of attack and elevators to adapt significant changes in power states and the extensive envelope range.Second, the predictedα~(k+d) is used to identify DC60 at the d steps ahead moment [21], and a certain relationship [19, 22], some empirical curves or look-up tables, between DC60 and the pressure ratio loss ΔPrs, is employed to get the relative value at the d steps ahead moment. It follows that(6)DC~60(k+d)=DC60(Ma,α~(k+d)),(7)ΔP~rs(k+d)=(DC~60(k+d)).Next, from (5) engine surge margin loss at k+d moment is estimated as
(8)ΔS~m(k+d)≈ΔP~rs(k+d)(1+Sm(k)).Finally, the surge margin in distortion states atk+d moment can be obtained as
(9)S~m(k+d)=Sm(k)-ΔS~m(k+d).Remark 2.
S
m
(
k
) is used to calculate ΔSm(k+d) in (8) if d steps are too short for phugoid dynamics in flight motions to influence the Sm baseline model. The detailed designation and validation for surge margin baseline model, angle of attack prediction model and direct engine surge margin control are, respectively, clarified as follows.


## 2.1. A Novel Direct Surge Margin Control

An integrated aircraft and engine dynamic model depicted in Figure1 is introduced to implement some necessary validations [6, 7]. This comprehensive model is composed by a dynamic aircraft model and its flight controller, a component-level turbofan engine model and its controller which will be designed as the proposed controller, and a conventional one for comparisons. For the aircraft model, the dynamics in longitudinal plane is only employed for simplicity, which can simulate level flight, climbing, accelerating, descending and super maneuver tasks. For the engine model, A8 and Wfb are selected as the control parameter, whereas Smf and Nc are the relative controlled variables. Therefore, the engine command r includes the references value of Smf and Nc, and e is a deviation vector with respect to command signal. A data set about the relationship between α and DC60 that is a type of distortion index are originated from flight test data of a F/A-18A aircraft in the Flight Research Center in NASA Dryden [19]. As discussed above, if the surge margin for an engine can be predicted accurately, a direct engine surge margin control, or a high stability engine control, easily can be further realized. In our research, a novel scheme is proposed (see Figure 1), where k represents the current moment and St is a sequential engine state values. The main idea is as follows.(a)
Considering the hysteresis between pilot input and rudder output, the current manipulations will affect the outputs ind-step time delay (or τ seconds in continuous sense) and the prediction for α is very important in a super maneuver task in which α changes very violently [1–3] to have a great impact on the surge margin. Therefore, two parameters are needed to control this high stability control scheme. One is S~mf(k+d), that is a predicted surge margin S~mf at the d-step ahead time, to adapt the variation in Smf due to the changes in α. The other need to control is Nc which regulates the power of the engine to supply enough thrust F.(b)
So when a super maneuver of the fighter is carried out, how to estimateS~mf(k+d) is a core problem in this scheme. For calculating S~mf(k+d), a time series of angle of α and the elevator δz are taken as inputs of the angle of attack prediction model to predict the d-step ahead value α~(k+d). Necessary engine state values in St series are obtained through the relative sensors to mapping the baseline value of Smf. After that, the α~(k+d) and St sequence are gathered as the inputs of fan surge margin prediction model; consequently the predicted surge margin S~mf(k+d) would be gotten in real time. For showing clearly, an equivalent pilot input delay system is labeled with a dotted rectangle in Figure 1.Figure 1
Structure for direct engine surge margin control.Remark 1.
In our research, the simulation step is preset as 20 ms; sod=τ/0.02=0.5/0.02=25 is predefined here based on [1–3].


## 2.2. Necessity to Estimate Surge Margin

For a turbofan engine, a parameter or a concept is defined as surge margin to show a physical meaning of a relative distance from the cooperate working point to the surge border (see Figure2), and the following equation is generally used to represent the surge margin of a fan, and similar for a compressor [20]:
(1)Sm=πs/Wcor,s-πo/Wcor,oπo/Wcor,o=πsπo·Wcor,oWcor,s-1.Figure 2
Surge margin change in a fan characteristic map.When engine inlet distortions happen due to high attack angle states, the surge border of the engine will be shifted down or the operation line will be shifted up, leading to a surge margin loss. Take the example of the shifting-down of surge border, and we have:(2)ΔSm=Sm-S~m=πs/Wcor,s-πo/Wcor,oπo/Wcor,o-πl/Wcor,s-πo/Wcor,oπo/Wcor,o=πs/Wcor,s-π1/Wcor,lπo/Wcor,o,
where ΔSm shows surge margin loss and S~m represents the surge margin when entering distortion; the subscript 1 shows the lowered stable boundary value.As can be seen from the above equations,Sm and ΔSm are parameters that cannot be measured. Under distorted states, if the undistorted Sm and the surge margin loss ΔSm might be accurately predicted, the distorted surge margin S~m will be expressed as
(3)S~m=Sm-ΔSm.For calculatingΔSm, an indirect parameter, called pressure ratio loss of the cooperating work point, is usually utilized according to the [21, 22]. Due to the shifting down of the boundary line, the pressure ratio loss ΔPrs is expressed as
(4)ΔPrs=πs/Wcor,s-π1/Wcor,lπs/Wcor,s.With the combination of (1) and (2) and the definition of ΔPrs, the surge margin loss is easily expressed as
(5)ΔSm=πs/Wcor,s-π1/Wcor,lπo/Wcor,o=πs/Wcor,s-π1/Wcor,lπs/Wcor,sπs/Wcor,sπo/Wcor,o=ΔPrs(1+Sm).


## 2.3. Surge Margin Prediction Model

Based on the current and historical information of input parameters, dynamic solutions for the baseline valueSm and the loss value ΔSm can be acquired or estimated by a novel engine surge margin model for super-maneuver flight, which consists of a surge margin baseline model and surge margin loss model. The modeling process is shown in Figure 3 and described as follows.Figure 3
Engine surge margin modeling in distortion states.


### 2.3.1. Surge Margin Baseline Model

In the range of super-maneuver flight envelope (height of 3~7 km, march number of 0.3~0.7), a surge margin baseline model based on a BP neural network, which is under no inlet distortion of the fighter, is trained to predict surge margin baseline value. Consider a lot of and redundant measurable parameters correlated with surge margin, a feature selection algorithm, for choosing adaptable engine sensors to estimate surge margin, is proposed by a weighted LSSVR method. With a contribution criterion of each parameter to surge margin, a group of most related and affordable parameters are selected as the model inputs.


### 2.3.2. Surge Margin Modeling in Distortion States

On the basis of surge margin baseline model, surge margin modeling in distortion states can be implemented in the following steps.First, an online scrolling window LSSVR model (OSW-LSSVR) is proposed to predict thed-step ahead value α~(k+d). The OSW-LSSVR model is capable of updating the prediction model based on the time series of angle of attack and elevators to adapt significant changes in power states and the extensive envelope range.Second, the predictedα~(k+d) is used to identify DC60 at the d steps ahead moment [21], and a certain relationship [19, 22], some empirical curves or look-up tables, between DC60 and the pressure ratio loss ΔPrs, is employed to get the relative value at the d steps ahead moment. It follows that(6)DC~60(k+d)=DC60(Ma,α~(k+d)),(7)ΔP~rs(k+d)=(DC~60(k+d)).Next, from (5) engine surge margin loss at k+d moment is estimated as
(8)ΔS~m(k+d)≈ΔP~rs(k+d)(1+Sm(k)).Finally, the surge margin in distortion states atk+d moment can be obtained as
(9)S~m(k+d)=Sm(k)-ΔS~m(k+d).Remark 2.
S
m
(
k
) is used to calculate ΔSm(k+d) in (8) if d steps are too short for phugoid dynamics in flight motions to influence the Sm baseline model. The detailed designation and validation for surge margin baseline model, angle of attack prediction model and direct engine surge margin control are, respectively, clarified as follows.


## 2.3.1. Surge Margin Baseline Model

In the range of super-maneuver flight envelope (height of 3~7 km, march number of 0.3~0.7), a surge margin baseline model based on a BP neural network, which is under no inlet distortion of the fighter, is trained to predict surge margin baseline value. Consider a lot of and redundant measurable parameters correlated with surge margin, a feature selection algorithm, for choosing adaptable engine sensors to estimate surge margin, is proposed by a weighted LSSVR method. With a contribution criterion of each parameter to surge margin, a group of most related and affordable parameters are selected as the model inputs.


## 2.3.2. Surge Margin Modeling in Distortion States

On the basis of surge margin baseline model, surge margin modeling in distortion states can be implemented in the following steps.First, an online scrolling window LSSVR model (OSW-LSSVR) is proposed to predict thed-step ahead value α~(k+d). The OSW-LSSVR model is capable of updating the prediction model based on the time series of angle of attack and elevators to adapt significant changes in power states and the extensive envelope range.Second, the predictedα~(k+d) is used to identify DC60 at the d steps ahead moment [21], and a certain relationship [19, 22], some empirical curves or look-up tables, between DC60 and the pressure ratio loss ΔPrs, is employed to get the relative value at the d steps ahead moment. It follows that(6)DC~60(k+d)=DC60(Ma,α~(k+d)),(7)ΔP~rs(k+d)=(DC~60(k+d)).Next, from (5) engine surge margin loss at k+d moment is estimated as
(8)ΔS~m(k+d)≈ΔP~rs(k+d)(1+Sm(k)).Finally, the surge margin in distortion states atk+d moment can be obtained as
(9)S~m(k+d)=Sm(k)-ΔS~m(k+d).Remark 2.
S
m
(
k
) is used to calculate ΔSm(k+d) in (8) if d steps are too short for phugoid dynamics in flight motions to influence the Sm baseline model. The detailed designation and validation for surge margin baseline model, angle of attack prediction model and direct engine surge margin control are, respectively, clarified as follows.


## 3. How to Set Up a Surge Margin Baseline Model

As discussed above, the surge margin baseline modeling should firstly identify the model inputs. Based on the LSSVR algorithm [23, 24], a weighted feature selection algorithm is proposed so as to screen a group of most valuable parameters that can be measured as inputs of the baseline model. Then, a neural network with nonlinear mapping capability is used for the surge margin baseline model design. Three parts, containing the preliminary LSSVR algorithm, surge margin feature selection and baseline model, are deduced and stated respectively as follows.


### 3.1. LSSVR Algorithm

To solve a nonlinear regression problem for a training data set{(xi,yi)}i=1,…,M, where xi∈Rm is the input with m-dimension and yi∈R is its corresponding output.A least square support vector regression (LSSVR) for this training data set can be transformed to an optimization problem as(10)minw,eiJ(w,ei)=12wTw+γ2∑i=1Mei2s.t.yi=wTφ(xi)+b+ei,i=1,…,M,
where w represents the model complexity, b is the offset, e=[e1,e2,…,eM] represents prediction residual vector, γ∈R+ is a regularization parameter, and φ(·) is a nonlinear mapping which can transform the input data into a high-dimensional feature space. In order to solve the above problem, a Lagrange function without constraints can be constructed as
(11)L(w,b,ei,αi)=J-∑i=1Mαi{wTφ(xi)+b+ei-yi},
where αi is Lagrange multiplier.Through Karush-Kuhn-Tucker (KKT) derivation, a Wolfe dual optimization problem for (10) can be expressed as
(12)minb,αi{L(b,αi)=[12∑i,j=1Mαiαjk(xi,xj)+12γ∑i=1Mαi2-∑i=1M∑i,j=1Mαiαjk(xi,xj)αiyi+b∑i=1Mαi]},
where k(xi,xj) is a kernel function chosen as exp{-∥xi-xj∥2/2ν2}.After that, (12) can be reformulated for convenience:
(13)min{L(b,α)=[12[αTb][K11T0][αb]-[αTb][Y0]]}.Then the optimal solution for (14) is easily obtained as
(14)[α*b*]=[K11T0]-1[Y0]=R[Y0],
where Y=[y1,y2,…,yM]T, α=[α1,α2,…,αM]T, 1=[11,12,…,1M]T, and K is the kernel matrix in which
(15)Kij=k(xi,xj)+δijγ=φT(xi)φ(xj)withδij={1,i=j0,i≠j.After obtaining the solutionα and b by (14), for any new testing sample x∈Rm, a regression or a predictor for x is gotten as
(16)f(x)=∑i=1Mαi*k(xi,x)+b*.


### 3.2. Feature Selection for Different Sensors

An aeroengine is not allowable to install too many sensors; otherwise it will degrade overall performance and increase computational complexity of the control systems. Therefore, a problem must be solved before setting up a surge margin baseline model; that is, which measurable variables will be selected as input variables for the model. In this section, an improved criterion for ranking variables is induced below [25, 26].Considering that each measurable component inx has different affordable fractions such as maintainability and price for variant sensors, a weighted matrix is necessary to be introduced here to quantify the influence of the affordable aspects for each parameter. Therefore, every sample or input should be transformed into the following expression:
(17)x′=Λx,
where Λ=ΛT>0 is defined as an affordable weighted matrix and proper dimension, and the problem (16) is changed as follows:
(18)min{L(b,α)=[12[αTb][K′11T0][αb]-[αTb][Y0]]},
where K′ is a weighted kernel matrix in which
(19)Kij′=k(xi′,xj′)+δijγwithδij={1,i=j0,i≠j.So from (14), the optimal value L* of (17) is easily gotten as
(20)L*=-12[YT0][K′11T0]-1[Y0]=-12[YT0][α*′b*′]=-12YTα*′.In the calculation process, if the optimal value ofL, L* is obtained, the surge margin prediction has the best precision. If the ith variable is removed, L-*(i)=-1/2YTα-(i), where α-(i) is the solution of (18) without the ith input variable. Thus, a criterion for ranking variables can be put forward:
(21)Υ(i)=L-*L*=YTα-*′(i)YTα*′.
If the value of Υ(i) is smaller than the value of Υ(j)  (j≠i), the ith variable is considered to make less contribution to the optimal value L-* than the jth variable. As for the aeroengine, there are more than 20 measurable variables for the selection. Then, each Υ(i) for different sensors should be computed based on (21), so all the variables can be ranked and some of them might be selected as the most contributors to surge margin baseline model. At last, the inputs St of Smf baseline model are chosen as Nf, Wfb, T22, P2, and P3 in Table 1, and the ones of Smc baseline model are Nc, T3, πc,T22, and  P3 (see Table 2).Table 1
Effects by feature selection for fan surge margin model.
Order
1
2
3
4
5
Others
Sensors
N
f
W
fb
T
22
P
3
P
2
⋯
Υ
(
i
)
9.2179
1.5517
1.5298
1.1396
1.1178
⋯Table 2
Effects by feature selection for  compressor surge margin model.
Order
1
2
3
4
5
Others
Sensors
N
c
T
3
π
c
T
22
P
3
⋯
Υ
(
i
)
12.3796
2.1351
1.4773
1.3526
1.2584
⋯


### 3.3.Sm Baseline Model

After the determination for the input parameters of the surge margin baseline model, a three-layer BP neural network [27, 28], shown in Figure 4, can be employed for the model designing, and its mathematical expression is
(22)Y=NN(X).Figure 4
The topological structure for BP model.Figure4 shows the topological structure of the three-layer BP neural network, in which there are p input points for input layer, q hidden points for hidden layer, and r out layer points. wij is denoted as connecting weighted coefficient between input and hidden layers, and wjk denoted as connecting weighted coefficient between the input and hidden layers. BP model can be described as follows.For the hidden layer, the followed equations can be acquired for thejth point:
(23)H(j)=fH(netj),netj=∑i=1pwijX(i),
where fH(·) is denoted as an excitation function for hidden layer, i=1 to p and j=1 to q.And for the output layer, the following expressions can also be gotten for themth point:
(24)netm=∑j=1qwjmH(j),Y(m)=gO(netm),
where gO(·) is denoted as an excitation function for output layer and m=1 to r. f(·) and g(·) are both chosen as sigmoid functions.If an evaluation function is defined asE=1/2∑m=1r(d(m)-Y(m))2, where d(m) is an expected output for the network, the adjusting rule for those weighted coefficients is deduced as follows:
(25)Δwjk=-η∂E∂wjk,Δwij=-η∂E∂wij,
where the scalar η>0 is a factor about convergent rate. For more details about BP algorithm, one can refer to [26].Overfitting often faces designers of intelligent identification methods like neural network, and it means that a more complex structure does not mean a more accurate mapping for all the inputs. So, the number of middle layer neurons needs to modulate with care to acquire the better choice, when the number of middle layer neurons in theSmf model is adjusted as 8, 14, and 16, and the relative test error is calculated as 7.58e−3, 6.38e−3, and 6.68e−3 separately in Figure 5. Similarly, the number in the Smc model is modulated as 7, 12, and 13, and the relative test error is 6.73e−4, 5.69e−4, and 5.79e−4, respectively (see Figure 6). So we choose 14 as the best number of middle layer neurons for the Smf model, and 12 for the Smc model.Precision test of different middle layer neurons of surge margin model.
(a)
Test result ofSmf model with 8 middle layer neurons
(b)
Test result ofSmf model with 16 middle layer neurons
(c)
Test result ofSmf model with 14 middle layer neuronsPrecision test of different middle layer neurons of surge margin model.
(a)
Test result ofSmc model with 7 middle layer neurons
(b)
Test result ofSmc model with 13 middle layer neurons
(c)
Test result ofSmf model with 12 middle layer neuronsFigure7 presents a comparison between the predicted and real Smf and Smc in Matlab environment. The relative testing errors are within 5% and 1%, respectively, indicating surge margin models have satisfactory prediction precision. Moreover, for testing the accuracy and real-time ability in VC++ environment which is our destination simulation environment, Figure 8 illustrates the simulation results between the outputs from real nonlinear plant and its prediction model for comparisons, and these results discover that the prediction model is capable of tracking the relative surge margin satisfactorily (dynamic error within 5% for Smf and 1% for, Smc), where “pre” means the predicted one.Precision test for surge margin prediction model in Matlab environment.
(a)
Test results ofSmf model
(b)
Test results ofSmfFigure 8
Precision test of surge margin prediction model in VC++ environment.


## 3.1. LSSVR Algorithm

To solve a nonlinear regression problem for a training data set{(xi,yi)}i=1,…,M, where xi∈Rm is the input with m-dimension and yi∈R is its corresponding output.A least square support vector regression (LSSVR) for this training data set can be transformed to an optimization problem as(10)minw,eiJ(w,ei)=12wTw+γ2∑i=1Mei2s.t.yi=wTφ(xi)+b+ei,i=1,…,M,
where w represents the model complexity, b is the offset, e=[e1,e2,…,eM] represents prediction residual vector, γ∈R+ is a regularization parameter, and φ(·) is a nonlinear mapping which can transform the input data into a high-dimensional feature space. In order to solve the above problem, a Lagrange function without constraints can be constructed as
(11)L(w,b,ei,αi)=J-∑i=1Mαi{wTφ(xi)+b+ei-yi},
where αi is Lagrange multiplier.Through Karush-Kuhn-Tucker (KKT) derivation, a Wolfe dual optimization problem for (10) can be expressed as
(12)minb,αi{L(b,αi)=[12∑i,j=1Mαiαjk(xi,xj)+12γ∑i=1Mαi2-∑i=1M∑i,j=1Mαiαjk(xi,xj)αiyi+b∑i=1Mαi]},
where k(xi,xj) is a kernel function chosen as exp{-∥xi-xj∥2/2ν2}.After that, (12) can be reformulated for convenience:
(13)min{L(b,α)=[12[αTb][K11T0][αb]-[αTb][Y0]]}.Then the optimal solution for (14) is easily obtained as
(14)[α*b*]=[K11T0]-1[Y0]=R[Y0],
where Y=[y1,y2,…,yM]T, α=[α1,α2,…,αM]T, 1=[11,12,…,1M]T, and K is the kernel matrix in which
(15)Kij=k(xi,xj)+δijγ=φT(xi)φ(xj)withδij={1,i=j0,i≠j.After obtaining the solutionα and b by (14), for any new testing sample x∈Rm, a regression or a predictor for x is gotten as
(16)f(x)=∑i=1Mαi*k(xi,x)+b*.


## 3.2. Feature Selection for Different Sensors

An aeroengine is not allowable to install too many sensors; otherwise it will degrade overall performance and increase computational complexity of the control systems. Therefore, a problem must be solved before setting up a surge margin baseline model; that is, which measurable variables will be selected as input variables for the model. In this section, an improved criterion for ranking variables is induced below [25, 26].Considering that each measurable component inx has different affordable fractions such as maintainability and price for variant sensors, a weighted matrix is necessary to be introduced here to quantify the influence of the affordable aspects for each parameter. Therefore, every sample or input should be transformed into the following expression:
(17)x′=Λx,
where Λ=ΛT>0 is defined as an affordable weighted matrix and proper dimension, and the problem (16) is changed as follows:
(18)min{L(b,α)=[12[αTb][K′11T0][αb]-[αTb][Y0]]},
where K′ is a weighted kernel matrix in which
(19)Kij′=k(xi′,xj′)+δijγwithδij={1,i=j0,i≠j.So from (14), the optimal value L* of (17) is easily gotten as
(20)L*=-12[YT0][K′11T0]-1[Y0]=-12[YT0][α*′b*′]=-12YTα*′.In the calculation process, if the optimal value ofL, L* is obtained, the surge margin prediction has the best precision. If the ith variable is removed, L-*(i)=-1/2YTα-(i), where α-(i) is the solution of (18) without the ith input variable. Thus, a criterion for ranking variables can be put forward:
(21)Υ(i)=L-*L*=YTα-*′(i)YTα*′.
If the value of Υ(i) is smaller than the value of Υ(j)  (j≠i), the ith variable is considered to make less contribution to the optimal value L-* than the jth variable. As for the aeroengine, there are more than 20 measurable variables for the selection. Then, each Υ(i) for different sensors should be computed based on (21), so all the variables can be ranked and some of them might be selected as the most contributors to surge margin baseline model. At last, the inputs St of Smf baseline model are chosen as Nf, Wfb, T22, P2, and P3 in Table 1, and the ones of Smc baseline model are Nc, T3, πc,T22, and  P3 (see Table 2).Table 1
Effects by feature selection for fan surge margin model.
Order
1
2
3
4
5
Others
Sensors
N
f
W
fb
T
22
P
3
P
2
⋯
Υ
(
i
)
9.2179
1.5517
1.5298
1.1396
1.1178
⋯Table 2
Effects by feature selection for  compressor surge margin model.
Order
1
2
3
4
5
Others
Sensors
N
c
T
3
π
c
T
22
P
3
⋯
Υ
(
i
)
12.3796
2.1351
1.4773
1.3526
1.2584
⋯


## 3.3.Sm Baseline Model

After the determination for the input parameters of the surge margin baseline model, a three-layer BP neural network [27, 28], shown in Figure 4, can be employed for the model designing, and its mathematical expression is
(22)Y=NN(X).Figure 4
The topological structure for BP model.Figure4 shows the topological structure of the three-layer BP neural network, in which there are p input points for input layer, q hidden points for hidden layer, and r out layer points. wij is denoted as connecting weighted coefficient between input and hidden layers, and wjk denoted as connecting weighted coefficient between the input and hidden layers. BP model can be described as follows.For the hidden layer, the followed equations can be acquired for thejth point:
(23)H(j)=fH(netj),netj=∑i=1pwijX(i),
where fH(·) is denoted as an excitation function for hidden layer, i=1 to p and j=1 to q.And for the output layer, the following expressions can also be gotten for themth point:
(24)netm=∑j=1qwjmH(j),Y(m)=gO(netm),
where gO(·) is denoted as an excitation function for output layer and m=1 to r. f(·) and g(·) are both chosen as sigmoid functions.If an evaluation function is defined asE=1/2∑m=1r(d(m)-Y(m))2, where d(m) is an expected output for the network, the adjusting rule for those weighted coefficients is deduced as follows:
(25)Δwjk=-η∂E∂wjk,Δwij=-η∂E∂wij,
where the scalar η>0 is a factor about convergent rate. For more details about BP algorithm, one can refer to [26].Overfitting often faces designers of intelligent identification methods like neural network, and it means that a more complex structure does not mean a more accurate mapping for all the inputs. So, the number of middle layer neurons needs to modulate with care to acquire the better choice, when the number of middle layer neurons in theSmf model is adjusted as 8, 14, and 16, and the relative test error is calculated as 7.58e−3, 6.38e−3, and 6.68e−3 separately in Figure 5. Similarly, the number in the Smc model is modulated as 7, 12, and 13, and the relative test error is 6.73e−4, 5.69e−4, and 5.79e−4, respectively (see Figure 6). So we choose 14 as the best number of middle layer neurons for the Smf model, and 12 for the Smc model.Precision test of different middle layer neurons of surge margin model.
(a)
Test result ofSmf model with 8 middle layer neurons
(b)
Test result ofSmf model with 16 middle layer neurons
(c)
Test result ofSmf model with 14 middle layer neuronsPrecision test of different middle layer neurons of surge margin model.
(a)
Test result ofSmc model with 7 middle layer neurons
(b)
Test result ofSmc model with 13 middle layer neurons
(c)
Test result ofSmf model with 12 middle layer neuronsFigure7 presents a comparison between the predicted and real Smf and Smc in Matlab environment. The relative testing errors are within 5% and 1%, respectively, indicating surge margin models have satisfactory prediction precision. Moreover, for testing the accuracy and real-time ability in VC++ environment which is our destination simulation environment, Figure 8 illustrates the simulation results between the outputs from real nonlinear plant and its prediction model for comparisons, and these results discover that the prediction model is capable of tracking the relative surge margin satisfactorily (dynamic error within 5% for Smf and 1% for, Smc), where “pre” means the predicted one.Precision test for surge margin prediction model in Matlab environment.
(a)
Test results ofSmf model
(b)
Test results ofSmfFigure 8
Precision test of surge margin prediction model in VC++ environment.


## 4. How to Set Up an Angle of Attack Prediction Model

As discussed above, the inlet distorted index caused by flight angle of attack determines the loss of the surge margin. For adapting the changes of angle of attack that cannot be predefined due to time delay, an OSW-LSSVR model is proposed to estimate thed step ahead value and it can pull new samples into an online training set, called a scrolling window, to adapt large variations in engine states and envelope points. The OSW-LSSVR algorithm and angle of attack prediction model are, respectively, described as follows.


### 4.1. OSW-LSSVR Algorithm

In order to realize an online LSSVR, after obtaining a new training data pair(xn,yn), a new predictor is needed to be reconstructed according to (14). For reducing the computational complexity to solve (14) directly, an iterative strategy [29, 30] is adopted here.A conceptP is firstly introduced as an index set, and corresponding training data of which are utilized to construct the predictor at nth iteration. |·| represents the cardinality. And the corresponding training data set or scrolling window set is defined as
(26)xP={xiP},|P|≤l-,
where l∈N is denoted as the storage capacity of training samples in the scrolling window.So for thenth iteration, it follows that
(27)Rn=[KPPkPnkPnTknn]-1,
where KPP is a relative weighted kernel matrix in which
(28)KPP,ij=k(xiP,xjP)+δijγwithδij={1,i=j,0,i≠j,knn=k(xn,xn)+1γ,kPn=(k(x1P,xn)+1γ⋮k(xiP,xn)+1γ⋮k(x|P|P,xn)+1γ).Then,Rn can be obtained as
(29)Rn=[Rn-100T0]+λ[β-1][βT-1],
where β=Rn-1[1kPn], λ=(knn-[1kPnT]TβT)-1,  x1P∈xP.Assuming thata and b at the n-1th iteration are computed with the equation [αPn-1bn-1]=Rn-1[YP0], the formulation of computing a and b at the nth iteration can be expressed as
(30)[αPnαnbn]=Rn[YPyn0]=[Rn-1[YP0]0]+λ(βT[YP0]-yn)[β-1]=[αPn-1bn-10]+λ(βT[YP0]-yn)[β-1].From (29) and (30), R, a, and b can be efficiently updated, and the predictor at the nth irritation is
(31)fn(x)=∑i∈PαP,ink(xiP,x)+bn.Thenth predictor in (31) is utilized to estimate the output value of xn, denoted as y~n. If a criterion followed is defined as:
(32)|yn-y~n|<ε,
where ε is a threshold which can control the tradeoff between the prediction accuracy and the parsimoniousness.Then,xn will be discarded. Otherwise, it will be chosen as a new support vector to the scrolling window.Remark 3.
Different from the above online LSSVR proposed in [30] in which only considering current samples. For meeting a d step ahead prediction model, some special improvement must be made and described as follows.(a)
As the output value in the future moment is unknown, the criterion proposed above cannot be realized when directly using|α(k+d)-α~(k+d)|<ε. Provided that the current attack angle has closer features with the one needing to predict at the d-step ahead moment, the criterion at k+d moment is replaced with the one at current moment. In other words, |α(k)-α~(k)|<ε is taken as the criterion of prediction model. For showing the precision more clearly, a relative deviation e(k) is taken as a new threshold rule as follows:
(33)e(k)=|α(k)-α~(k)α(k)|×100%<ε′.(b)
For enhancing the real-time ability for the online LSSVR, thel- for a scrolling window should be kept constant, but in [30] the value is not limited. Meanwhile, a simple but efficient judge logic is adopted to determine whether the new sample or support vector can be accepted as the support vector. If the number of training samples is still not more than l-, xn will be taken as support vector directly when xn is added into the scrolling window. Otherwise, xminP would be given away before xn is put into the scrolling window, where
(34)xminP=xiP|min(|yiP-y^iP|<ε).(c)
Asd step ahead information is impossible to get, some offline samples need to be gathered in advance for starting the online training process. In this study case, as said above not less than 25 samples are needed.


### 4.2. Angle of Attack Prediction Model

With the above OSW-LSSVR algorithm, the online angle of attack prediction model can be devised, and a 3-order prediction model with OSW-LSSVR is set up as follows:(35)Y′(k)=AOA(X′(k-d)),
where X′=[α(k-d-2),α(k-1-d),α(k-d),δz(k-d)], and Y′=[α~(k),α~(k-1)].Thus, for predicting the value, thek+d moment is easily gotten as
(36)Y′(k+d)=AOA(X′(k)),
where X′(k)=[α(k-2),α(k-1),α(k),δz(k)], and Y′(k+d)=[α~(k+d),α~(k+d-1)].As discussed above, the threshold determination is conducted through (33) to decide whether the scrolling window should be updated. If e(k) is smaller than ε′, it means that the prediction has a satisfactory precision, and α(k+d) can be accurately predicted without updating the scrolling window and vice versa. Modulated by trial and error, parameters in the OSW-LSSVR model are as follows.Rolling window capacityl-=80, threshold value ε′=0.0048, kernel parameter v=1.6, and penalty factor γ=224. Where the selection principle of l- is to keep the rolling window capacity as small as possible, so as to optimize the instantaneity of the online model.As shown in Figures9 and 10, in the envelop point of H=4 km, Ma=0.5, and H=5km, Ma=0.55, the angle of attack prediction model is verified, respectively. For simulating the dynamic change of angle of attack in super maneuver flight, the elevator is sufficiently excited to get a variation range of angle of attack from −10° to +50°. The correlation curves are shown in Figures 9(a) and 10(a), including about 1990 groups of samples for the actual and predicted one. The model can quickly and accurately predict the angle of attack at k+d moment. And also seen clearly from Figures 9(b) and 10(b), the relative testing errors are within 1%, which means a much better prediction effect compared to the effects reported in the literatures [6, 7].Model precision test atH=4 km, Ma=0.5.
(a)
Test result ofα(k+d) prediction model
(b)
Relative test error ofα(k+d)Model precision test atH=5 km, Ma=0.55.
(a)
Test result ofα(k+d) prediction model
(b)
Relative test error ofα(k+d)


## 4.1. OSW-LSSVR Algorithm

In order to realize an online LSSVR, after obtaining a new training data pair(xn,yn), a new predictor is needed to be reconstructed according to (14). For reducing the computational complexity to solve (14) directly, an iterative strategy [29, 30] is adopted here.A conceptP is firstly introduced as an index set, and corresponding training data of which are utilized to construct the predictor at nth iteration. |·| represents the cardinality. And the corresponding training data set or scrolling window set is defined as
(26)xP={xiP},|P|≤l-,
where l∈N is denoted as the storage capacity of training samples in the scrolling window.So for thenth iteration, it follows that
(27)Rn=[KPPkPnkPnTknn]-1,
where KPP is a relative weighted kernel matrix in which
(28)KPP,ij=k(xiP,xjP)+δijγwithδij={1,i=j,0,i≠j,knn=k(xn,xn)+1γ,kPn=(k(x1P,xn)+1γ⋮k(xiP,xn)+1γ⋮k(x|P|P,xn)+1γ).Then,Rn can be obtained as
(29)Rn=[Rn-100T0]+λ[β-1][βT-1],
where β=Rn-1[1kPn], λ=(knn-[1kPnT]TβT)-1,  x1P∈xP.Assuming thata and b at the n-1th iteration are computed with the equation [αPn-1bn-1]=Rn-1[YP0], the formulation of computing a and b at the nth iteration can be expressed as
(30)[αPnαnbn]=Rn[YPyn0]=[Rn-1[YP0]0]+λ(βT[YP0]-yn)[β-1]=[αPn-1bn-10]+λ(βT[YP0]-yn)[β-1].From (29) and (30), R, a, and b can be efficiently updated, and the predictor at the nth irritation is
(31)fn(x)=∑i∈PαP,ink(xiP,x)+bn.Thenth predictor in (31) is utilized to estimate the output value of xn, denoted as y~n. If a criterion followed is defined as:
(32)|yn-y~n|<ε,
where ε is a threshold which can control the tradeoff between the prediction accuracy and the parsimoniousness.Then,xn will be discarded. Otherwise, it will be chosen as a new support vector to the scrolling window.Remark 3.
Different from the above online LSSVR proposed in [30] in which only considering current samples. For meeting a d step ahead prediction model, some special improvement must be made and described as follows.(a)
As the output value in the future moment is unknown, the criterion proposed above cannot be realized when directly using|α(k+d)-α~(k+d)|<ε. Provided that the current attack angle has closer features with the one needing to predict at the d-step ahead moment, the criterion at k+d moment is replaced with the one at current moment. In other words, |α(k)-α~(k)|<ε is taken as the criterion of prediction model. For showing the precision more clearly, a relative deviation e(k) is taken as a new threshold rule as follows:
(33)e(k)=|α(k)-α~(k)α(k)|×100%<ε′.(b)
For enhancing the real-time ability for the online LSSVR, thel- for a scrolling window should be kept constant, but in [30] the value is not limited. Meanwhile, a simple but efficient judge logic is adopted to determine whether the new sample or support vector can be accepted as the support vector. If the number of training samples is still not more than l-, xn will be taken as support vector directly when xn is added into the scrolling window. Otherwise, xminP would be given away before xn is put into the scrolling window, where
(34)xminP=xiP|min(|yiP-y^iP|<ε).(c)
Asd step ahead information is impossible to get, some offline samples need to be gathered in advance for starting the online training process. In this study case, as said above not less than 25 samples are needed.


## 4.2. Angle of Attack Prediction Model

With the above OSW-LSSVR algorithm, the online angle of attack prediction model can be devised, and a 3-order prediction model with OSW-LSSVR is set up as follows:(35)Y′(k)=AOA(X′(k-d)),
where X′=[α(k-d-2),α(k-1-d),α(k-d),δz(k-d)], and Y′=[α~(k),α~(k-1)].Thus, for predicting the value, thek+d moment is easily gotten as
(36)Y′(k+d)=AOA(X′(k)),
where X′(k)=[α(k-2),α(k-1),α(k),δz(k)], and Y′(k+d)=[α~(k+d),α~(k+d-1)].As discussed above, the threshold determination is conducted through (33) to decide whether the scrolling window should be updated. If e(k) is smaller than ε′, it means that the prediction has a satisfactory precision, and α(k+d) can be accurately predicted without updating the scrolling window and vice versa. Modulated by trial and error, parameters in the OSW-LSSVR model are as follows.Rolling window capacityl-=80, threshold value ε′=0.0048, kernel parameter v=1.6, and penalty factor γ=224. Where the selection principle of l- is to keep the rolling window capacity as small as possible, so as to optimize the instantaneity of the online model.As shown in Figures9 and 10, in the envelop point of H=4 km, Ma=0.5, and H=5km, Ma=0.55, the angle of attack prediction model is verified, respectively. For simulating the dynamic change of angle of attack in super maneuver flight, the elevator is sufficiently excited to get a variation range of angle of attack from −10° to +50°. The correlation curves are shown in Figures 9(a) and 10(a), including about 1990 groups of samples for the actual and predicted one. The model can quickly and accurately predict the angle of attack at k+d moment. And also seen clearly from Figures 9(b) and 10(b), the relative testing errors are within 1%, which means a much better prediction effect compared to the effects reported in the literatures [6, 7].Model precision test atH=4 km, Ma=0.5.
(a)
Test result ofα(k+d) prediction model
(b)
Relative test error ofα(k+d)Model precision test atH=5 km, Ma=0.55.
(a)
Test result ofα(k+d) prediction model
(b)
Relative test error ofα(k+d)


## 5. Direct Engine Surge Margin Control

Based on the above surge margin prediction model, a direct surge margin control can be implemented, and the structure of it has been shown in Figure2. A quasi-PID or ILQR method [31] is used to design this novel controller. Here, we just consider the surge margin of the fan as an application. And the ILQR robust control method, direct surge margin controller design, and the digital simulation are, respectively, introduced as follows.


### 5.1. ILQR Robust Control

Provided is that an engine dynamics with surge margin estimation can be formulated as follows:(37)x˙=Ax+G1(x)+Bu,y=Cx+G2(x)+ω,
where x, y, and u are denoted as system state, output, and input vector, respectively, and G1(x), G2(x) are defined as nonlinear compensation items for system state and output dynamics accordingly. Note that G1(x) and G2(x) represent differences between a nonlinear engine model and its linear one, due to variations in states and flight conditions, and ω represents the estimation and measured errors produced in surge margin prediction.Take the augmented state vector asx-=[xT∫oteTdτ]T=[xT∫ot(r-y)Tdτ]T, while r is a set-point command (37) which is easily transformed to
(38)x-˙=A-x-+G-1(x-)+B-u-,y-=C-x-+G-2(x-)+ω-,
where u-=u, A-=[A0-C0], B-=[B-D], C-=[C0], G-1(x-)=[G1(x)-G2(x)], G-2(x-)=G2(x), and ω-=ω.Generally for aeroengines, a bounded condition is guaranteed as [32]
(39)∥[G1(x)G2(x)+ω]∥∞<+∞.Based on the lemma proposed in [31], for (38) a relative LQR regulator can be easily designed so that all states can converge to a bounded area, that is,
(40)∥x-∥∞=∥[xT∫0teTdτ]T∥∞<+∞.And it means thate∈L∞. Furthermore, e∈L2 and e∈C1 can also be guaranteed for engines. Therefore based on Barbalat lemma, the output can track the command with limited conditions followed by
(41)limt→∞e=limt→∞(r-y)=0,whenE(ω)=0.Otherwise,(42)limt→∞e=limt→∞(r-y)<σ,when|E(ω)|=σ.Consider a well-known LQR problem for the system (37) as
(43)minu-{∥Q1/2x-R1/2u-∥2=∫0∞(x-TQx-+u-Tu-)dt},
where Q=QT≥0 and R=RT>0 are proper dimension weighted matrices. Thus, an ILQR controller can be derived as u-=K-x- where K-=R-1B-TP and P satisfies the following Riccati equation
(44)A-TP+PA--PB-R-1B-TP+Q=0.K
- can be divided by x and ∫0tedτ as a block matrix K-=[KxKe]; so we obtain the following quasi-PID control type [33]:
(45)u=K-x-=[KxKe][x∫0tedτ]=Kxx+Ke∫0tedτ.
And the structure of an ILQR controller is shown in Figure 11.Figure 11
ILQR control structure.Remark 4.
As well known, LQR method has better robustness with infinite magnitude margin and phase margin over 60°. However, LQR controllers cannot have the ability to eliminate the steady state error in control problems to trace some command signals. For nonlinear plants such as aeroengines with significant nonlinearity and wide variant states range, simple robust linear methods like LQR is not able to maintain a good controllability [34, 35]. So an improved type of the above ILQR method is utilized here. In an ILQR control design, the tracing errors are augmented into the state vector so as to realize a convergence to the references signals. So the ILQR control method can be regarded as an improved LQR type; therefore it has not only has robustness to adapt large variation of engine power states and envelope changes, but also has capability to eliminate the steady state error when engine command signals change.


### 5.2. Design of Direct Surge Margin Control

Firstly, at the design point ofH=0 km, Ma=0, and Pla=70°, a state-space model of the engine dynamics is established. And the actuators dynamics of u=[Wfb,A8]T is also considered in it. Thus, a nonlinear model augmented with the actuator is gotten as follows:(46)x˙=Ax+G1(x)+Bu,y=Cx+G2(x)+ω,
where x=[Nf,Nc,Wfb,A8]T, u=[Wfb,A8]T, and y=[Nf,Smf]T. Through some conventional identification method, the system matrices in (44) can be obtained as follows:
(47)A=[-2.32641.33501.41600.90590.3056-3.90530.16130.592900-100000-5],B=[000010005],C=[0.01210.9597-0.00420.00178.3219-13.79044.1979-0.6527].For reducing the number of parameters needing to modulate, positive diagonal matrices are chosen forQ=QT≥0 and R=RT>0, in which the former is related to output performance of the controlled system, and the latter is related to control input performance of the controlled system. The ranges of two matrices are roughly given to meet that of a stabilized closed-loop system. Then, they can be modulated by trial and error according to the dynamic and static performance of the closed system. Ultimately, the two weighted matrices are selected as Q=diag[1,1,0.1,0.2,10,2] and R=diag[0.5,0.5], where diag[] means diagonal matrices:
(48)Ke=[29.509512.03537.2756-0.3396],Kx=[-26.75818.2457-5.87000.2909-0.6300-2.77120.0010-0.4108].The control effects of the above controller are verified firstly in the design point, as shown in Figure12, and the commands are set as Smf=20%, Nc=99%. Figures 12(a) and 12(b) present the response curves of Smf and Nc, and Figures 12(c) and 12(d) show the response curves of Wfb and A8. As can be seen from these figures, when a step command signal is given at t=5 s, the closed system can response quickly to reach the control objective within 4 seconds. Besides, there are no static errors for two controlled variables (see Figures 12(a) and 12(b)).Effects of the ILQR control at the design point.
(a)
Response curves of fan surge margin
(b)
Response curves of compressor corrected speed
(c)
Response curves of nozzle area
(d)
Response curves of main fuel flowFor validating the robustness of this ILQR controller, some test results are also given in other envelope points in the scope of the super maneuver envelope range. It can be seen clearly that at the envelope point of (H=4 Km, Ma=0.6) and (H=5 Km, Ma=0.5), the closed system designed by ILQR method still has excellent dynamic and static performance. Figures 13(a) and 13(b) show a fast, small overshoot, and convergent tracking with respect to command signals, whereas Figures 13(c) and 13(d) illustrate acceptable and feasible control input performances. These results give a proof in good robustness of this control method, and it is qualified for application in the direct surge margin control under variant power states and at different envelope point.Control effects in other envelope points.
(a)
Response curves of fan surge margin
(b)
Response curves of compressor corrected speed
(c)
Response curves of nozzle area
(d)
Response curves of main fuel flow


## 5.1. ILQR Robust Control

Provided is that an engine dynamics with surge margin estimation can be formulated as follows:(37)x˙=Ax+G1(x)+Bu,y=Cx+G2(x)+ω,
where x, y, and u are denoted as system state, output, and input vector, respectively, and G1(x), G2(x) are defined as nonlinear compensation items for system state and output dynamics accordingly. Note that G1(x) and G2(x) represent differences between a nonlinear engine model and its linear one, due to variations in states and flight conditions, and ω represents the estimation and measured errors produced in surge margin prediction.Take the augmented state vector asx-=[xT∫oteTdτ]T=[xT∫ot(r-y)Tdτ]T, while r is a set-point command (37) which is easily transformed to
(38)x-˙=A-x-+G-1(x-)+B-u-,y-=C-x-+G-2(x-)+ω-,
where u-=u, A-=[A0-C0], B-=[B-D], C-=[C0], G-1(x-)=[G1(x)-G2(x)], G-2(x-)=G2(x), and ω-=ω.Generally for aeroengines, a bounded condition is guaranteed as [32]
(39)∥[G1(x)G2(x)+ω]∥∞<+∞.Based on the lemma proposed in [31], for (38) a relative LQR regulator can be easily designed so that all states can converge to a bounded area, that is,
(40)∥x-∥∞=∥[xT∫0teTdτ]T∥∞<+∞.And it means thate∈L∞. Furthermore, e∈L2 and e∈C1 can also be guaranteed for engines. Therefore based on Barbalat lemma, the output can track the command with limited conditions followed by
(41)limt→∞e=limt→∞(r-y)=0,whenE(ω)=0.Otherwise,(42)limt→∞e=limt→∞(r-y)<σ,when|E(ω)|=σ.Consider a well-known LQR problem for the system (37) as
(43)minu-{∥Q1/2x-R1/2u-∥2=∫0∞(x-TQx-+u-Tu-)dt},
where Q=QT≥0 and R=RT>0 are proper dimension weighted matrices. Thus, an ILQR controller can be derived as u-=K-x- where K-=R-1B-TP and P satisfies the following Riccati equation
(44)A-TP+PA--PB-R-1B-TP+Q=0.K
- can be divided by x and ∫0tedτ as a block matrix K-=[KxKe]; so we obtain the following quasi-PID control type [33]:
(45)u=K-x-=[KxKe][x∫0tedτ]=Kxx+Ke∫0tedτ.
And the structure of an ILQR controller is shown in Figure 11.Figure 11
ILQR control structure.Remark 4.
As well known, LQR method has better robustness with infinite magnitude margin and phase margin over 60°. However, LQR controllers cannot have the ability to eliminate the steady state error in control problems to trace some command signals. For nonlinear plants such as aeroengines with significant nonlinearity and wide variant states range, simple robust linear methods like LQR is not able to maintain a good controllability [34, 35]. So an improved type of the above ILQR method is utilized here. In an ILQR control design, the tracing errors are augmented into the state vector so as to realize a convergence to the references signals. So the ILQR control method can be regarded as an improved LQR type; therefore it has not only has robustness to adapt large variation of engine power states and envelope changes, but also has capability to eliminate the steady state error when engine command signals change.


## 5.2. Design of Direct Surge Margin Control

Firstly, at the design point ofH=0 km, Ma=0, and Pla=70°, a state-space model of the engine dynamics is established. And the actuators dynamics of u=[Wfb,A8]T is also considered in it. Thus, a nonlinear model augmented with the actuator is gotten as follows:(46)x˙=Ax+G1(x)+Bu,y=Cx+G2(x)+ω,
where x=[Nf,Nc,Wfb,A8]T, u=[Wfb,A8]T, and y=[Nf,Smf]T. Through some conventional identification method, the system matrices in (44) can be obtained as follows:
(47)A=[-2.32641.33501.41600.90590.3056-3.90530.16130.592900-100000-5],B=[000010005],C=[0.01210.9597-0.00420.00178.3219-13.79044.1979-0.6527].For reducing the number of parameters needing to modulate, positive diagonal matrices are chosen forQ=QT≥0 and R=RT>0, in which the former is related to output performance of the controlled system, and the latter is related to control input performance of the controlled system. The ranges of two matrices are roughly given to meet that of a stabilized closed-loop system. Then, they can be modulated by trial and error according to the dynamic and static performance of the closed system. Ultimately, the two weighted matrices are selected as Q=diag[1,1,0.1,0.2,10,2] and R=diag[0.5,0.5], where diag[] means diagonal matrices:
(48)Ke=[29.509512.03537.2756-0.3396],Kx=[-26.75818.2457-5.87000.2909-0.6300-2.77120.0010-0.4108].The control effects of the above controller are verified firstly in the design point, as shown in Figure12, and the commands are set as Smf=20%, Nc=99%. Figures 12(a) and 12(b) present the response curves of Smf and Nc, and Figures 12(c) and 12(d) show the response curves of Wfb and A8. As can be seen from these figures, when a step command signal is given at t=5 s, the closed system can response quickly to reach the control objective within 4 seconds. Besides, there are no static errors for two controlled variables (see Figures 12(a) and 12(b)).Effects of the ILQR control at the design point.
(a)
Response curves of fan surge margin
(b)
Response curves of compressor corrected speed
(c)
Response curves of nozzle area
(d)
Response curves of main fuel flowFor validating the robustness of this ILQR controller, some test results are also given in other envelope points in the scope of the super maneuver envelope range. It can be seen clearly that at the envelope point of (H=4 Km, Ma=0.6) and (H=5 Km, Ma=0.5), the closed system designed by ILQR method still has excellent dynamic and static performance. Figures 13(a) and 13(b) show a fast, small overshoot, and convergent tracking with respect to command signals, whereas Figures 13(c) and 13(d) illustrate acceptable and feasible control input performances. These results give a proof in good robustness of this control method, and it is qualified for application in the direct surge margin control under variant power states and at different envelope point.Control effects in other envelope points.
(a)
Response curves of fan surge margin
(b)
Response curves of compressor corrected speed
(c)
Response curves of nozzle area
(d)
Response curves of main fuel flow


## 6. Control Simulations and Validations

To validate effectiveness of the new method, some super-maneuver flight tasks are simulated with a flight condition in low altitudes and velocities, and for maintaining a stable and efficient operation, the fan surge margin command is usually expected to be 15%. The direct surge margin control is compared to the conventional control as shown in Figures14 and 15. The conventional engine controller is a bivariate controller designed by ILQR, and the control inputs are A8 and Wfb, and the controlled variables are πt and Nc. Figures 14(a) and 15(a) depict the response curves of fan surge margin Smf. “-no” indicates the conventional control response, “-DSC” represents the response of the direct surge margin control, and “-pre” represents the estimation of the angle of attack. Figures 14(b) and 15(b) illustrate different curves of the angle of attack α, the elevator angle δz, inlet distorted index DC60, fan and compressor surge margin loss, compressor surge margin Smc, relative corrected speed Nf, relative corrected speed Nc, main fuel flow Wfb, and throat area of nozzle A8.Effects of direct surge margin control atH=4 km, Ma=0.6.
(a)
Response curves of fan surge margin
(b)
Response curves of each parameter of the engineEffects of direct surge margin control atH=5 km, Ma=0.5.
(a)
Response curves of fan surge margin
(b)
Response curves of each parameter of the engineFigure14 demonstrates the simulation results under the condition of H=4 km, Ma=0.6. At t=0 s, the conventional control and the direct surge margin control are, respectively, adopted to begin a high attack angle task. Then, the elevator δz is pulled along with some predefined path file (a negative degree means an upward shift of the elevator). When the elevator shifts downwards, the angle of attack α is increased, so DC60 is accordingly increased and vice versa, and the variation of DC60 is monotonically related with the surge margin loss of fan and compressor. Moreover, from the attack angle response curves, one can see that a real time and precisely prediction for α at d steps ahead moment is verified here.As also can be observed from Figure14(a), while using the conventional control, Smf  changes very violently during the super maneuver flight, and in some time even less than 10% indicting that the cooperating work point is very close to the surge border. So, it fails to meet the safe requirements for engine in the super maneuver flight. On the contrary, the novel control can keep the surge margin Smf always around 15%, that is, an ideal engine stability control effect. In order to keep the surge margin at 15% overall the super maneuver flight, a combined regulation is applied via Wfb and A8. Since a direct correlation [1–3] between Smf and A8, the modulation of A8 serves as a more important role. So one can obviously observe that the bigger is the A8, the larger is the surge margin. But for the changes of Wfb, a bit similar tendency is found either in the new control or the conventional one, as the main physical contribution of it is the proper controlling of Nc, which can ensure that there is no great thrust loss in the super maneuver flight.For validating the robustness of this new control law, simulation results atH=5 km, Ma=0.5 are also given in Figure 15. Similarly through a fast adjustment of elevator angle δz, a maneuvering flight is originated from an level flight. Directly using the predicted surge margin, Figure 15 clearly presents that the new control scheme can successfully achieve a real-time adjustment of the cooperating point of the engine. In the entire super maneuver mission, the fan surge margin is controlled slightly around 15%; that is to say, the cooperating work point can always keep a safe distance from the surge border. Nevertheless, when using that conventional method, since surge margin is not well considered or better estimated, safety margin of the engine cannot be sustained in the large maneuver flight.


## 7. Conclusions

Considering the special stability need in super-maneuver flights, a novel direct surge margin controller is proposed and designed based on a surge margin prediction model and a ILQR robust control method which is capable of guaranteeing a good robustness and performance in a large variation of power states and flight conditions for the closed-loop system.In terms of the surge margin prediction model, a weighted surge margin feature selection algorithm based on LSSVR is proposed to get a group of optimal measurable inputs for its baseline model, and an angle of attack model with an online scrolling window LSSVR method is proposed to estimate the inlet distortion index in super maneuver states. Furthermore, the surge margin in distortion states could be computed in real time by aid of the surge margin baseline and angle of attack prediction.All these contributions are verified to effectively realize a high stability control with the flight conditions of low altitude, low speed, and high angle of attack.



---

*Source: 870215-2013-12-25.xml*