An improved multivariable generalized predictive control algorithm for direct performance control of gas turbine engine

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.1 (1-17)

Aerospace Science and Technology ••• (••••) ••••••

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

An improved multivariable generalized predictive control algorithm for direct performance control of gas turbine engine Xin Zhou, Feng Lu, Wenxiang Zhou, Jinquan Huang ∗ Jiangsu Province Key Laboratory of Aerospace Power System, College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

a r t i c l e

i n f o

Article history: Received 24 March 2019 Received in revised form 30 September 2019 Accepted 17 November 2019 Available online xxxx Keywords: Gas turbine engine Direct performance control Multivariable generalized predictive control Quantum-behaved particle swarm optimization Limit protection

a b s t r a c t Providing thrust for the aircrafts is the primary task of the gas turbine engine. Accurate and safe thrust controller design has always been the focus of the research. In this paper, an improved control algorithm for direct performance control of the gas turbine engine is proposed to realize the on-line estimation and tracking of the performance parameters. The controlled plant, identification module, controller and diagnosis module are integrated as a complete system, which has two loops. First, the adaptive model calculates the engine’s performance parameters (thrust and compressor surge margin) as the feedback of the controller; then a Quantum-behaved Particle Swarm Optimization (QPSO) based multivariable generalized predictive controller (MGPC) is adopted as the main controller. To solve complex nonlinear optimization problems with constraints, the basic QPSO algorithm is improved from two aspects, global average best position and initial global best position. In addition, a penalty factor is added to the cost function to realize the limit protection of the rotor speed and exhaust gas temperature (EGT). Compared with the traditional control methods, the main contribution of this paper to realize the direct performance control without the conversion by the rotor speed, and the control system is responsible for limit protection and fault diagnosis in the meanwhile. Finally, simulations are performed to investigate the response of the performance parameters, the effects of the controller parameters, the ability of fault tolerance of the controller, and robustness against measurement noise and model uncertainties. The results show that the proposed scheme is robust and can achieve accurate regulation of the performance parameters and the limited outputs within the safe range, whether there is any fault or not. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction As the power unit of the airplane, the most important role of the gas turbine engine is to provide the thrust for the aircraft. Traditional thrust control systems of the gas turbine engine are usually based on sensors [1]. The performance parameters (such as the thrust) and the parameters that are significant to the engine safe operation (such as the surge margin) are selected as the output controlled variables. But these parameters are unmeasurable and hard to estimate accurately during the flight. Instead, some measureable state parameters that can reflect the magnitude and variation of performance parameters are used. The most common is to control the engine according to the rotor speed and pressure ratio. Litt has done a lot of relative research [1,2], in which fan speed command was modified to obtain a desired thrust based on throttle position. Csank and May proposed a speed and pres-

*

Corresponding author. E-mail address: [email protected] (J. Huang).

https://doi.org/10.1016/j.ast.2019.105576 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

sure ratio control system design process of a generic commercial turbofan engine, including a gain-scheduled controller and a limit protection module [3]. Montazeri-Gh et al. designed a fuzzy fuel controller using genetic algorithm for aero-engine thrust regulation based on the rotor speed [4]. However, these control methods have some problems; in fact, the mathematical relationship between the observed variables and the control variables is not fixed, but may change with the deterioration of the engine [5,6]. The thrust - speed conversion table cannot be applied to the whole life cycle of the engine; thus the pilots need to manually adjust the power level angle (PLA), which will increase their workload. In addition, since the performance parameters are unmeasurable, a large margin should be reserved in the design process to ensure the safety of the gas turbine engine. In order to solve these problems, the traditional engine control structure is gradually transformed to a model-based control system. The performance parameters of the engine are calculated by the adaptive model according to the measurements from the sensors; thus these performance parameters can be used as the feedback to realize the direct performance control. Intelligent engine control (IEC),

JID:AESCTE

AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.2 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

2

Nomenclature A8 CARIMA CLM E EGT H HPT FG f J Kp Ki Kd LPT Ma MGPC N Ny Nu P PSO QPSO R RLS RMSE

nozzle throat area controller auto regression integrated moving average component level model expectation exhausted gas temperature altitude high-pressure turbine thrust function cost function proportionality coefficient integral coefficient differential coefficient low-pressure turbine Mach number multivariable generalized predictive control rotor speed prediction domain control domain pressure particle swarm optimization quantum-behaved particle swarm optimization measurement noise covariance matrix of the nonlinear filter recursive least squares root mean square error

as an advanced control concept ranked in the top four, one of the core contents of its research is direct thrust control [7]. Litt proposed a double-loop control system structure, the inner loop controls fan speed while the outer loop automatically adjusts the engine’s fan speed command to maintain the thrust at the desired level; this is based on pilot input and the estimated thrust, even as the engine deteriorates with use [8,9]. Ring and Henriksson developed the concept of the thrust control and simplified the thrust control loop to a single input single output system with the estimator [10]. Chen et al. added another limit protection controller besides the thrust tracking controller [11]. These methods are still based on the relationship between the thrust and the rotor speed. The inner loop controller serves as the main controller, and the estimated thrust is only used to modify the speed reference. It cannot be called direct thrust control in the full sense. Besides, for the limit protection of the surge margin, temperature and pressure that are related to the engine safe operation, additional controllers need to be designed; thus the complexity of the system will be increased. A two-variable integrated control system including a diagnostic module and a modeling module will be proposed in this paper. The controlled variables are thrust and compressor surge margin (S M Comp ); the adaptive model is modified by a nonlinear filter online; the unmeasured performance parameters calculated by the model are regarded as that of the real engine; the pilots provide thrust and S M Comp commands directly for the controller to produce the response. As for the advanced control algorithm of the gas turbine engine, many scholars have done the research [12–15]. Predictive control is an excellent candidate for gas turbine engines, due to its inherent capacity for control of constrained nonlinear systems. The performance improvement of constrained nonlinear model predictive control (NMPC) with state and parameter estimation over traditional control architectures is investigated and applied to a model turbofan aircraft engine by Brunell et al. [16]. Du et al.

SE SM SW T W U UKF Y

 ξ

η

efficiency scalers of the rotating components surge margin flow capacity scalers of the rotating components temperature flow capacity control vector unscented Kalman filter output vector weight coefficient of the control input uncorrelated random process noise efficiency

Subscripts cor Comp d f L H r 22 3 43 5 8

corrected value compressor design point fuel low-pressure shaft high-pressure shaft reference compressor inlet compressor exit high-pressure turbine exit low-pressure turbine exit nozzle throat

applied a novel multivariate constrained predictive control technique based on linear state space model to a commercial turbofan engine [17]. An online, fully NMPC for a gas turbine using the reduced-order, internal model was presented by Wiese et al., and can be extended to any gas turbine system [18]. A robust estimator and a robustly-tuned generalized predictive controller are incorporated to present a robust adaptive scheme for a linear time varying (LTV) system of a launch vehicle [19]. Zheng et al. proposed an engine NMPC, integrated with flight predictive information to obtain better engine response speed and reduces the rotor transient droop [20]. Among various predictive control algorithms, the generalized predictive control (GPC) has attracted many extensive interests and been successfully applied in industry during the last few years [21]. Generalized predictive control, also known as multi-step predictive self-tuning control, maintains the advantage of online identification, output prediction and minimum output variance of the minimum variance self-tuning control, while absorbing the rolling optimization strategy in predictive control. It combines the superiority of both adaptive control and predictive control [22,23]. The GPC algorithm designed for multiple input multiple output (MIMO) predictive models results in the MGPC algorithm. MGPC obtains the optimal control vector by solving the corresponding constraint optimization problem, which has inspired many fast optimization algorithms. Sedraoui et al. proposed an improvement of the particle swarm optimization (PSO) algorithm to iteratively resolve the cost problem of the MGPC under multiple constraints [24]. Alkorta et al. presented the design method of a new linear MGPC for sensorless induction motor drives [25]. Smoczek and Szpytko developed a novel control approach based on a MGPC with particle swarm optimizer for limiting the transient process of the plant [26]. And a system identification and GPC method of a class of MIMO models are studied by Ye and Lou et al. [27]. Recently, more and more intelligent optimization algorithms are used to solve constrained MGPC problems. A novel global-convergence-guaranteed QPSO algorithm is proposed by Sun

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.3 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

and Xu et al. [28]. In this paper, this basic QPSO algorithm will be improved to solve the constraint MGPC problem to realize the direct performance control for the gas turbine engines in real time. The main contribution of this paper is to propose an improved MGPC controller and applied it to the direct thrust control of the gas turbine engine. The performance parameters calculated by the adaptive component-level model are taken as the controlled variables directly. And the prediction model is a Controller Auto Regression Integrated Moving Average (CARIMA) model which is identified online by the least squares algorithm with forgetting factor. Meanwhile, the penalty factor is added to handle the limitations when designing the cost function, namely, both the performance tracking and limit protection are performed by one controller. To verify the validity of the proposed algorithm, we compare the performance of the improved MGPC and PID controller, then simulate its robustness and fault tolerance ability. This paper is organized as follows. Section 2 gives a brief introduction about the unconstrained MGPC algorithm; the novel constrained MGPC algorithm for gas turbine engine direct performance control is described in Section 3; Section 4 introduces QPSO based MGPC control system integration; some simulation results and analyses are given in Section 5, while Section 6 concludes the paper. 2. Unconstrained MGPC algorithm









A z−1 Y (k) = z−d B z−1 U (k) + ξ(k)/

(1)

where U (k) ∈ Rn×1 := [u 1 (k), u 2 (k), · · · , un (k)] T , Y (k) ∈ Rm×1 := [ y 1 (k), y 2 (k), · · · , ym (k)] T are the inputs and outputs of the system, respectively. d denotes the delay of the system. ξ(k) is the uncorrelated random process noise and  = 1 − z−1 the difference operator. The coefficient matrices have the following forms [22]:



A ( z−1 ) = I n×n + A 1 z−1 + A 2 z−2 + · · · + A na z−na B ( z−1 ) = B 0 + B 1 z−1 + B 2 z−2 + · · · + B nb z−nb

(2)

where na , nb , denote the highest orders of A ( z−1 ), B ( z−1 ). Without loss of generality, assuming the delay of the system is d = 1, Equation (1) can be simplified to:









¯ z−1 Y (k) = B z−1 U (k − 1) + ξ(k) A ¯ (z−1 ) = A (z−1 ) = A (z−1 )(1 − z−1 ), and where A

(3)

U (k) = U (k) −

U (k − 1) is the difference of U (k). The goal of MGPC is to make the actual outputs track the reference inputs as fast as possible and to reduce the error of the system. Consider the following quadratic cost function [29]:

J=E

 Ny 

Nu  2  2 Y (k + j ) − Y r (k + j ) + γ j U (k + j − 1)

j =1

Fig. 1. Sequence of control structure diagram.

reference trajectory usually converted by a first-order lag model (also called a first-order smoothing model) as shown Equation (5).





j =1

(4) where N y is the prediction domain, N u the control domain, Y r the reference trajectory and γ j the weight of the j-th control input. The process of MGPC mainly includes the following stages. 1) Smoothing In predictive control, in order to smoothly convert the predicted outputs to the expected values at a certain response speed, the

Y r (k) = Y (k) Y r (k + j ) = α Y r (k + j − 1) + (1 − α )ω, j = 1, 2, . . . , N y

(5)

where α ∈ [0, 1) is the softening factor, ω the expected value and Y r the reference values entered into the controller. 2) Multistep Prediction To derive the optimal output prediction sequence, the following Diophantine equation is considered [30]:



Generalized predictive control, also known as multi-step predictive self-tuning control, was firstly proposed by Clarke in 1987 [22,23]. A CARIMA model for an n inputs and m outputs multivariable process can be expressed by [29]:

3

¯ ( z −1 ) + z − j G j ( z −1 ) I = E j ( z −1 ) A F j ( z −1 ) = E j ( z −1 ) B ( z −1 ) = L j ( z −1 ) + z − j H j ( z −1 )

(6)

where

⎧ E j ( z−1 ) = I + E j ,1 z−1 + E j ,2 z−2 + · · · + E j ,nej z−nej ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ G j ( z−1 ) = G j ,0 + G j ,1 z−1 + · · · + G j ,n z−n g j gj

⎪ L j ( z−1 ) = L j ,0 + L j ,1 z−1 + · · · + L j ,nlj z−nlj ⎪ ⎪ ⎪ ⎪ ⎩ H j ( z−1 ) = H j ,0 + H j ,1 z−1 + · · · + H j ,nhj z−nhj , ne j , n g j , nlj , nhj are the highest orders of corresponding matrix. Premultiplying Equation (3) by E j ( z−1 ) and using Equation (6), we have:

Y (k + j ) = E j ξ(k + j ) + F j U (k + j − 1) + G j Y (k)

(7)

where Y (k) = [Y (k), Y (k − 1), · · · , Y (k − na )] T . As E j ξ(k + j ) is uncorrelated with other terms so that the optimal predictor can be expressed as:













Yˆ (k + j |k) = L j z−1 U + G j z−1 Y (k) + H j z−1 U (k − j ) (8) where U = [U (k), U (k + 1), · · · , U (k + N u − 1)] T is the current and future control increment vector, and U (k − j ) = [U (k − 1), U (k − 2), · · · , U (k − nb )] T the previous control increment vector. In the sequences mentioned above, the time ranges correspond to Y (k), Yˆ (k + j |k), U , U (k − j ) are shown in Fig. 1. 3) Optimization Solution The principle of solving the optimal control law is to minimize the cost function J . The control increment sequences are obtained by solving ∂ J /∂ U = 0 [30].

 −1    U (k) = L T L + L T Y r (k) − G Y (k) − H U (k − j ) where

(9)

JID:AESCTE

AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.4 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

4

Fig. 2. Schematic representation of a gas turbine engine.

⎡

L 1,0

⎢ ⎢ L 2,0 ⎢ L=⎢ ⎢ .. ⎣ . L N2 ,N1

⎡

···

0

L 2,1

···

0

.. .

..

.. .

···

···

.

L N2 ,N u

⎤

···

H 1,nb

H 2,1

···

H 2,nb ⎥

.. .

..

H N 2 ,0

H N 2 ,1

···

G 1,0

G 1,1

···

G 1,na

G 2,1

···

G 2,na ⎥

.. .

..

G N 2 ,1

···

⎢ ⎢ H 2,0 ⎢ H =⎢ . ⎢ . ⎣ .

⎢ ⎢ G 2,0 ⎢ G =⎢ . ⎢ . ⎣ . G N 2 ,0

⎥ ⎥ ⎥ ⎥ ⎦

.. .

.

.

where f refers to the nonlinear function, Y ∈ R2×1 is the output vector which are the thrust (F G) and the surge margin of the compressor (S M Comp ), and U ∈ R2×1 is the control input vector which are fuel flow (W f ) and nozzle throat area ( A 8 ). The predictive model shown as Equation (1) can be obtained by linearizing the nonlinear process mentioned above through online identification technology.

⎥ ⎥ ⎥ ⎥, ⎥ ⎦

H 1,1

H 1,0

⎡

0

⎤

3.1.2. On-line identification and output predictor The process of CARIMA modeling is equivalent to that of linearizing the engine at the current operating point. With the component characteristics of the engine, the thermodynamic parameters of the design point and engine rotational inertia, the inputoutput dynamic relationship of a double-shaft engine under the condition of small deviation can be deduced to have the following second-order transfer function form [31]:

and

H N 2 ,nb

.. .

⎤ ⎥ ⎥ ⎥ ⎥ ⎦

G (s) =

G N 2 ,na

are the solution of Diophantine equations, and = γ I (nN u ) is the weight coefficient matrix of control variables. According to the rolling optimization and feedback correction principle, the current control input is:

U (k) = U (k − 1) + U (k)

 − 1 = U (k − 1) + diT L T L + L T   × Y r (k) − G Y (k) − H U (k − j )

where

⎡

1 0

diT = ⎣ 0 0

...

... 0 0 0 ... 0

(10)

⎤

... 0 1 0 ... 0 ⎦ ... ... ... ... n×2N



A z

B z

−1

−1

 

 = I 2×2 +  =

0

0

a122

b011

b012

b021

b022





+





a111

z

−1

+

b111

b112

b121

b122

a211

0

0

a222



 z −2 (13)

z

−1

where the first number of the subscripts of the matrix elements refers to the order of the coefficient matrix, and the latter two mark the position of the elements in the matrix. Let:

u

3.1. Process model and output prediction 3.1.1. Process description In this paper, a MGPC algorithm for direct performance control is proposed with a dual-rotor turbofan engine as the controlled object. Fig. 2 shows the schematic diagram of the turbofan engine. As can be seen, the main components of the engine are: an inlet, a fan, a compressor, a combustor, a High-Pressure Turbine (HPT), a Low-Pressure Turbine (LPT), a mixing chamber, a bypass, and a nozzle. The operation process of the engine to be controlled can be assumed to a two-input two-output nonlinear system described by the following discrete-time model:



Y (t ) = f Y (t − 1), . . . Y (t − N y ); U (t − 1), . . . , U (t − N u )

+ ξ(t )

(12)

The bivariate CARIMA model of the gas turbine engine is assumed as a second-order system. The parameters of the model are estimated online using the Recursive Least Squares (RLS) algorithm with forgetting factor. To simplify the solution of Diophantine equation and reduce the parameters need to be estimated, the output coefficient matrix is set as a diagonal matrix:



3. Constrained MGPC algorithm for gas turbine engine direct performance control



δ y (s) b1 s + b2 = 2 δ u (s) a1 s + a2 s + a3

(11)

 θ1 = a111 a211 b011 b111 b012  θ2 = a122 a222 b021 b121 b022  ϕ1 = − y 1 (k − 1) − y 1 (k − 1) T u 2 (k − 1) u 2 (k − 2)  ϕ2 = − y 2 (k − 1) − y 2 (k − 1) T u 2 (k − 1) u 2 (k − 2)

b112 b122

T T

u 1 (k − 1) u 1 (k − 2) u 1 (k − 1) u 1 (k − 2) (14)

There is:

Y = ψT  + ξ  T ϕ1 where ψ = 0

0

ϕ2T



 ,  =



(15)

θ1 . In this way, model identifiθ2

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.5 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

cation is transformed into a standard least squares identification problem. The model parameters can be solved recursively by Equation (16).

⎧ ˆ ˆ ˆ ⎪ ⎨ (k) = (k − 1) + K (k)[ y (k) − ψ(k)(k − 1)] K (k) = P (k − 1)ψ T (k)[ψ(k) P (k − 1)ψ T (k) + λ I ]−1 ⎪ ⎩ P (k) = [ I − K (k)ψ(k)] P (k − 1)/λ

(16)

After solving the Diophantine equation with model parameters, the predictive outputs can be calculation by Equation (8). 3.2. Cost function with penalty factor The working temperature and pressure of the gas turbine are usually very high. In order to ensure its safe work, various restrictions should be met. In this paper, a cost function with penalty factor for nonlinear optimization algorithm is proposed to improve the MGPC controller. The advantages of this cost function are as follows: first, the output limitation management is realized; second, the burden of optimization calculation can be reduced by reducing the constraints; third, it is equivalent to imposing the terminal inequality constraints to the cost function, which ensures the stability of the controller. The target of control is to adjust the thrust and S M Comp to the desired value with W f and A 8 as inputs, while ensuring that the limited outputs: low-pressure rotor speed (N L ), high-pressure rotor speed (N H ), and EGT are within the preset range. For this purpose, the following quadratic cost function is proposed:

min J = E

 Ny 

  2 λj  S M Comp (k + j ) − S M Comp ,r (k + j ) j =1 Nu  

2 γ j U (k + j − 1)



j =1

0

⎢ ⎣ 0

0

⎡

⎢ = ⎣ B 21

=

⎥

⎢ ⎥ ⎣ NH ⎦

A 33 B 12

⎤

lim



⎥

B 22 ⎦ B 32

NL

EGT

W f  A8

lim

⎤



⎡

ξ1

⎤

⎢ ⎥ + ⎣ ξ2 ⎦ ξ3

(18) lim

ˆ L lim N



ˆ H lim N

= L lim 3 U (k) + G lim 3 Y lim 3 (k) + H lim 3 U (k − j ) (19)

(17)

U min (k + i ) ≤ U (k + i ) ≤ U max (k + i ) The cost function shown in Equation (17) mainly consists of three parts. The first part is to ensure the minimum variance between the output of the plant and a desired trajectory at the given future sampling points; the second part is to minimize the energy of the control quantity and in the third part, the penalty function is used to ensure that the parameters related to the safe operation of the engine do not exceed the limit. In Equation (17), F G and  S M Comp are the predictive outputs of the CARIMA model, F G r and S M Comp,r are the corresponding

T

desired tracking trajectories. U =  W f  A 8 is the control increment vector. N y is the prediction horizon and N u the control

T

⎡

= L lim 1−2 U (k) + G lim 1−2 Y lim 1−2 (k) + H lim 1−2 U (k − j ) Yˆ lim 3 = [ E Gˆ T ]

U min (k + i ) ≤ U (k + i ) ≤ U max (k + i )



0 B 11

⎤

0 ⎦

A 22



+ max{Yˆ lim 1 − Y lim 1_ max , 0}   S 1 × max{Yˆ lim 1 − Y lim 1_ max , 0} T  + max{Yˆ lim 2 − Y lim 2_ max , 0}   × S 2 max{Yˆ lim 2 − Y lim 2_ max , 0}



0

Yˆ lim 1−2

T



s.t .

A 11

Similar to the above, the model used to solve the Diophantine equations is identified using RLS. And the coefficient matrices of Diophantine are decomposed for rotor speed and EGT respectively, so that the limited predictive outputs are [29,30]:

j =1

+

⎡

B 31

Ny

+

and U max are the maximum of control vector and control increment vector. The weight coefficients S 1 and S 2 are diagonal positive definite matrices, which are usually 100 times larger than other weight coefficients at least. λ is the weight coefficient of the S M Comp relative to the F G, and γ the weight coefficient of control input. The penalty function dealing with the upper limit of the outputs refers to the larger value of each term of Yˆ lim − Y lim _ max and 0. If Yˆ lim < Y lim _ max , the corresponding output does not exceed the limit. Therefore the larger of the two is 0, which means the penalty function doesn’t work, and the optimization process is equivalent to that without constraints. Once the predictive outputs Yˆ lim exceed the corresponding limits Y lim _ max , even if the difference is small, this penalty function term will be large multiplied by the weight coefficient matrix S, which leads to the increase of J . The nonlinear optimization algorithm aims to solve the control law that makes J minimum. Thus, the output is maintained near the limit line eventually, ensuring the output prediction not to exceed the limit during the whole process. The CARIMA model for limited outputs prediction is constructed by selecting N L , N H , and EGT as outputs and W f and A 8 as inputs [22,29].

2

F G (k + j ) − F G r (k + j )

5

ˆ L Nˆ H horizon. Yˆ lim 1 = N and Yˆ lim 2 = E G T are the predictive outputs that need to be limited. Y lim 1_ max and Y lim 2_ max are upper limit of corresponding outputs. U min and U min are the minimum of control vector and control increment vector respectively, U max

where L lim , G lim , H lim refer to the decomposed matrices. At this point, the cost function can be solved by recursive optimization. 4. QPSO based MGPC control system integration 4.1. Structure of direct performance control system The constraint MGPC algorithm is applied in the gas turbine engine direct performance control system. The controller should have the ability to adjust the engine performance parameters to the reference value with excellent dynamic and static characteristics while satisfying the safe limitations. On the one hand, the filter estimates the health state of the gas turbine engine in real-time to modify the nonlinear adaptive model to obtain its unmeasurable performance parameters; on the other hand, least squares algorithm with forgetting factor is used to identify the CARIMA predictive model at the current working point to solve the control law. The structural block diagram of the control system is shown in Fig. 3.

JID:AESCTE

AID:105576 /FLA

6

[m5G; v1.261; Prn:27/11/2019; 15:03] P.6 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

Fig. 3. Structure of the control system.

The control system is mainly composed of four modules: 1) Plant: the real engine, which works according to the given input while providing the corresponding measurements; 2) Diagnostic module, including a nonlinear adaptive model and a UKF filter. The filter estimates the health parameters of gas turbine engine on the basis of the measurement residual between model and real engine to modify the adaptive model. The corrected model should track the output of real engine precisely so that the thrust and the S M Comp of the model can be regarded as that of real engine; 3) Identification module: the CARIMA model of engine is identified using RLS with forgetting factors online. Essentially, the model is linearized at the current operating point; 4) MGPC module: solve the optimal control law through multistep prediction and rolling optimization. Since the Identification module and MGPC controller haven been covered in detail earlier, the component level model and the Diagnostic Module will be briefly introduced here. 4.1.1. Component level model (CLM) The model used in this paper is built by the component method, that is, each component model has several inputs/outputs that form a separate module, regardless of the changes of the internal parameters of the component. Then, according to the matching conditions of each component, common operation equations are constructed and solved by Newton-Raphson approach. The component-level-model is written using C++ language and packaged with dynamic link library (DLL) for simulation in MATLAB environment. The modeling procedures are based on Ref. [32, 33], which provides the calculation formulas of each component, the selection of the initial guess vector, the construction and solution of the common operation equations, and the operation results of the model. 4.1.2. Diagnostic module The Diagnosis Module includes a component-level adaptive model and a nonlinear filter. Assuming the nonlinear aero-thermodynamic model of a turbofan engine is given by Equation (20):

xk+1 = f (xk , uk ) + w k yk = h(xk , uk ) + v k

(20)

where k is the time index, y the measurements, x the state vector representing for the engine health condition, and u the control

vector. The noise terms w and v represent the process inaccuracies and measurement inaccuracies in the model. During service of the engine, the components may experience faults and degeneration inevitably. To describe these changes, the health parameters are introduced to quantify the abrupt change and gradual degeneration degree of each component or the difference between the nominal model and the actual engine. The health parameters are selected as the variation coefficients of the efficiency or the flow capacity of the rotating components:

S Ei =

ηi , ηi∗

SWi =

Wi W i∗

(21)

where subscripts i = 1, 2, 3, 4 indicates the fan, the compressor, the HPT and the LPT, respectively. ηi and W i are the actual value of efficiency and flow capacity of the corresponding components under the current working condition, while ηi∗ and W i∗ are the ideal values of the efficiency and flow capacity of the corresponding components under the fault-free condition. For the controlled plant in this paper, eight health parameters are defined, i.e. x = [ S E 1 , S W 1 , S E 2 , S W 2 , S E 3 , S W 3 , S E 4 , S W 4 ] T . Turbofan engine is a typical nonlinear system. Considering the accuracy and computational complexity of the filtering algorithm, a UKF is designed in this article. UKF uses the statistical linearization technique to estimate the posterior mean and variance of random variables through the unscented transformation [34,35]. The filter estimates the health parameters of the components according to the residual sequence between the actual engine measurements and the nonlinear model outputs. Thus, the problem of the adaptive modeling is transformed into the problem of estimating the state parameters of the nonlinear system. The specific method and procedures of estimating the health parameters by the UKF could be referred to Ref. [36]. 4.2. Improved QPSO algorithm for control law solving It is hard to give a control law in an analytical form for the nonlinear programming problem mentioned above, instead, the intelligent optimization algorithm can be used. PSO, first introduced by Clerc and Kennedy [37], is a stochastic optimization algorithm based on cooperation and competition between groups and individuals. However, the global optimal solution cannot be always guaranteed since the particle velocity is limited. A novel global-

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.7 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

Algorithm 1: Diversity guided QPSO. 1) Initialize X i (0) on [ X min , X max ], let p_best i (0) = X i (0), and then set the global best position according to Equation (23).

1 M | A|

  M  N M  1   X i, j ,  ( X i , j − X¯ j )2 , X¯ j = i =1

j =1

√

M



C (t ) = C 1 (t ), C 2 (t ), . . . , C n (t )

3) Calculate the diversity of particle swarm according to: di versit y ( S X ) =

the closer a particle is to the optimal solution, the larger its weight coefficient. The mean best position C , therefore, is calculated as:



2) Find the fitness function value of each particle.



=

i =1

M 1 



where | A | = ( X max − X min ) N is the length of longest the diagonal in the search space. 4) Compare diversity (S X ) with the given threshold: if diversity (S X ) < the threshold g_best j = g_best j + δ A ε , ε ∼ N (0, 1) else Skip to step 5 end 5) Update the best individual position and best global position according to Equations (22) and (23). 6) Update the mean best position C according to C (t ) = (C 1 (t ), C 2 (t ), M . . . , C n (t )) = M1 i =1 p_best i (t ) 7) Update X i , j according to X i , j (t + 1) = p i , j (t ) ± β|C j (t ) − X i , j (t )| ln u i , j (t ) ∼ U (0, 1)

1 , u i , j (t )

convergence-guaranteed PSO—QPSO is developed [28], which introduces the concept of quantum evolution into the PSO algorithm. Due to the clustering of particle swarms, the algorithm will gradually lose the diversity and fall into local optimization, Sun and Xu et al. proposed a diversity guided QPSO [38,39]. In this section, the diversity guided QPSO is improved to solve the control law for the gas turbine engine direct performance control. In an N-dimensional target search space, particle swarm consists of M particles X = { X 1 , X 2 , . . . , X M }, the current position in the search space is X i (t ) = { X i ,1 (t ), X i ,2 (t ), . . . , X i , N (t )}, i = 1, 2, . . . , M, and X i , j (t ) ∈ [ X max , X min ], j = 1, 2, . . . , N. The personal best position is p_best i (t ) = { p_best i ,1 (t ), p_best i ,2 (t ), . . . , p_best i , N (t )} and the global best position is g_best i (t ) = { g_best i ,1 (t ), g_best i ,2 (t ), . . . , g_best i , N (t )}. The best individual position of particle i is determined by the Equation (22).

p_best i (t ) =

X i (t ),

f [ X i (t )] < f [ p_best i (t − 1)]

p_best i (t − 1),

f [ X i (t )] ≥ f [ p_best i (t − 1)] (22)

where f is the cost function. And best global position is determined by Equation (23).

g = arg min

1≤ i ≤ M

 



f p_best i (t )

(23)

g_best (t ) = P g (t ) The diversity guided QPSO is summarized as Algorithm 1. To improve the global search capability of the algorithm and shorten the calculation time, the following improvements will be made to diversity guided QPSO. (1) QPSO with Exponential Weighted Mean Best Position The average best position of the particle is used to evaluate the range of the search space and guide the evolution of the particle. In previous algorithm, C is the arithmetic mean of the best positions of all particles, which means that each particle has the same effect on it. In fact, different particles have different fitness values. In view of this, the arithmetic mean best position should be replaced with weighted mean best position. It is reasonable to rank the particles in descendent order according to their fitness value, and then the weight coefficients will be assigned to each particle,

ωi p_best i,1 (t ),

i =1

M 1 



M 1 



ωi p_best i,2 (t ), . . . ,

i =1

ωi p_best i,n (t )

(24)

i =1

where ωi indicates the weight coefficient of the i-th particle and  = iM=1 ωi . Exponential weight coefficients can arrange the importance of the particles in proportion. The following weight coefficients with exponential distribution is proposed:

ω2 = αω1 , ω3 = αω2 , . . . ωn = αωn−1 where

8) Repeat steps 2 through 5 until the termination condition is met (the optimal value is found or the maximum number of iterations is reached).



7

(25)

α is the exponential factor.

(2) Improvement on Initial Global Best Position In PSO and QPSO, the particles are initialized randomly since the region where the optimal solution will fall into cannot be known in the incipient stage when solving the optimal problem. However, for the practical issue of gas turbine engine direct performance control, the plant will only exceed the safe limits under some extraordinary conditions, for example, the engine encounters a component failure or the speed command is too high. In most working condition, the penalty term in the cost function will not come into effect, which means the optimal control law can be calculated as Equation (26).

∂J =0 ∂ U

(26)

The solution is:

− 1  U (k) = L 1T L 1 + L 2T  L 2 +

   · L 1T F G r (k) − G 1 F G (k) − H 1 U (k − j )   − L 2T  G 2 S MC (k) + H 2 U (k − j )

(27)

where  is the weight matrix of S M Comp relative to F G, and is the weight matrix of control input. L, G, H are solutions of Diophantine equation, subscript 1 indicates the part corresponding to the thrust, and subscript 2 that of S M Comp . In the case that the engine do not exceed its limits, Equation (27) is the optimal solution of this control problem. In view of this, the following modifications are made: the positions of M particles are initialized randomly as before, and the global best position is set according to Equation (28).



 −1

g_best (0) = U (k) = L 1T L 1 + L 2T  L 2 +

·







L 1T F G r (k) − G 1 F G (k) − H 1 U (k − j )

(28)



 − L 2T  G 2 S MC (k) + H 2 U (k − j )

Due to the diversity guided strategy, the particles barely premature. The modification on the mean best position and the initial global best position can highlight the contribution of the particles close to the optimal position and guide the particles evolving towards the optimal solution more quickly. The improved QPSO algorithm can comprehensively deal with the optimization problems no matter the engine exceed the limits or not, the process can be summarized as Algorithm 2.

JID:AESCTE

AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.8 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

8

Algorithm 2: Improved diversity guided QPSO. 1) Initialize X i (0) on [ X min , X max ], let p_best i (0) = X i (0), and then set the global best position according to Equation (28). 2) Find the fitness function value of each particle. 3) Calculate the diversity of particle swarm as in Algorithm 1. 4) Compare diversity (S X ) with the given threshold: if diversity (S X ) < the threshold g_best j = g_best j + δ A ε , ε ∼ N (0, 1) else Skip to step 5 end 5) Update the best individual position and best global position according to Equations (22) and (23). 6) Update the mean best position C according to Equation (24). 7) Update X i , j according to X i , j (t + 1) = p i , j (t ) ± β|C j (t ) − X i , j (t )| ln u 1(t ) , i, j

u i , j (t ) ∼ U (0, 1) 8) Repeat steps 2 through 5 until the termination condition is met (the optimal value is found or the maximum number of iterations is reached).

Fig. 4. Engine’s operating envelope.

Table 1 Available measurements instruction. Measurement

Acronyms

Standard deviation

Low pressure spool speed

NL

0.0015

High pressure spool speed

NH

0.0015

Fan exit temperature

T 22

0.002

Fan exit pressure

P 22

0.0015

Compressor exit temperature

T3

0.002

Compressor exit pressure

P3

0.0015

High-pressure turbine exit temperature

T 43

0.002

High-pressure turbine exit pressure

P 43

0.0015

Low-pressure turbine exit temperature

T5

0.002

Low-pressure turbine exit pressure

P5

0.0015

5. Simulation and analyses In this section, we will verify that the MGPC-based gas turbine engine direct performance control system designed above can handle the control problems under various conditions, and adjust the thrust and S M Comp to the reference value quickly and accurately. Besides, the engine ought to work within the safe range even when failure occurs. The simulation were carried out under design operating point (H = 0 m, Ma = 0, N H ,cor = 100%) to validate the effectiveness of this method. The inputs of engine are W f and A 8 , the controlled variables are F G and S M Comp , and measurements that filter needs are N L , N H , T 22 , P 22 , T 3 , P 3 , P 43 , T 43 , T 5 , and P 5 . The instruction and the standard deviation of the measurements are listed in Table 1.

Fig. 5. Outputs comparison curves.

5.1. CARIMA model verification The accuracy of the predictive model is the key to the MGPC controller. To verify the effectiveness of the CARIMA model, the complete operation of the engine within the flight envelope including take-off - acceleration - cruise - landing is simulated in this section. Fig. 4 gives the engine’s operation envelop. During the flight, parameters of CARIMA model are identified in real time and the output between the engine and the model are compared. This part of simulation only focuses on the model identification, thus the system is open-loop without a controller. The outputs comparison curves and CARIMA model parameters are shown in Fig. 5 and Fig. 6, respectively. From the simulation results, it can be seen that the CARIMA model can track the engine output with high accuracy within the envelope. When engine working condition changes, the model parameters can be adjusted adaptively to predict the thrust and S M Comp accurately. The average relative error between the model and the actual engine output thrust is 0.2836%, and the average

Fig. 6. CARIMA model parameters.

relative error of S M Comp is 2.375%, which is slightly larger due to the effect of the noise. Both of the errors can meet the requirement of controller design.

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.9 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

9

Fig. 7. Dynamic response of the MGPC controller at the design point of the gas turbine engine.

5.2. Response of performance parameters The control system was built for simulation according to Fig. 3. In order to verify the effectiveness of the controller, different thrust references are chosen to examine the ability of performance tracking and limit protection. The relative parameters are set as: S 1 = S 2 = 500, N u = 3, N y = 5,  = 1.2I 2×2 and = 1.2I 2×2 . This set of parameters is determined after several simulation comparisons, which is more empirical rather than mathematically. The specific selection basis includes setting time and overshoot. The simulation last for 30 seconds and consists of several processes of step response and ramp response. The upper limits of outputs are N L ≤ 105%N L ,d , N H ≤ 105%N H ,d and E G T ≤ 1600K , the subscript d indicates the design point. And considering engineering experience, after the normalization of the control variables, ±5% of their nominal value is taken as the limit of the control increment. The requirements of the control system are to ensure that the gas turbine engine can track the thrust reference, maintain the S M Comp unchanged and not exceed the safe limits. In addition, the improved MGPC controller will be compared with the widely used PID controller to illustrate its superiority. In the PID controller, the two outputs is controlled by W f and A 8 respectively. The coefficients are determined by trial-and-error search as follows. Thrust control loop: K p = 0.5, K i = 0.1, K d = 0.1; S M Comp control loop: K p = 0.5, K i = 0.01, K d = 0.01. The simulation results are shown in Fig. 7, (a) is the response of the performance parameter and (b) the trend of control variables. In Fig. 7, the thrust and control variables are normalized to the design point and S M Comp is the actual value. As can be seen from Fig. 7, the improved MGPC controller can adjust the actual thrust tracking the reference in a large range. There is barely overshoot during the transient process and the setting time is quite short; namely, the dynamic performance of the controller is satisfying. S M Comp is kept stable, ensuring the safety of the compressor during the acceleration and deceleration process and preventing the occurrence of surge. The changes of control variables are also within the reasonable range. The program running time of fault diagnosis and thrust estimation is 14.919s, the running time of CARIMA model identification and controller calculation is 9.891s, both of which can meet the real-time requirements. Compared with PID, the improved MGPC controller has shorter setting time, smaller overshoot and better dynamic performance. As can be seen from Fig. 7 (b), PID controller is more obvi-

ously affected by the measurement noise, which is manifested in more intense vibration of control variables. Besides, the improved MGPC has a better ability in handling the bivariate coupling problems than PID. The response curves show that when the thrust is stepped, the S M Comp will suffer a slight fluctuation and then stabilize at the reference value. However, the amplitude of S M Comp fluctuation controlled by PID is larger. To verify the effectiveness of the penalty factor in the cost function and test its ability of limit protection, another experiment is carried out. The maximum thrust reference is increased to 8% to make the EGT of the gas turbine engine exceed 1600K without limitation. Fig. 8 provides the comparison of performance tracking with and without penalty factor in cost function, (a) is the response of thrust parameters, (b) the trend of control variables, (c) the response of EGT and (d) the value of the cost function. As can be seen from Fig. 8. (a) and (c), when the limit protection module is not involved in the controller, it can achieve accurate tracking of thrust, however, the limit output EGT would obviously exceed its upper limit during the transient process. When the output limit is taken into account, the EGT operates at 1600K during the transient process and the gas turbine engine performance must be sacrificed to some degree to ensure the safety. The maximum thrust that the engine can provide is about 106% that of the design point. In this case, the value of the cost function will significantly increase due to the existence of the thrust error, as shown in Fig. 8 (d). From Fig. 8 (b), compared with A 8 , the difference between two groups of W f is more significant, which means W f is subjected to stricter restrictions due to its greater impact on the thrust. This example illustrates that the penalty factor can deal with outputs limit well. The PID controller doesn’t have the capability of limit protection, thus its response curves is similar to that of unconstraint MGPC except that the response is slower, the fluctuation is more severe, and the EGT will also be overtemperature. The common solution is to design additional limit protection controllers, such as the Min-Max controller, but this will undoubtedly increase the complexity of the system. Besides, PID is unable to adjust the controller parameters (K p , K i and K d ) adaptively, requiring the designer to adjust for different operation points manually. While MGPC is a kind of self-adaptive control and can calculate the controller parameters within the operating envelope autonomously. Further, the influence of limitation on the steady-state performance, mainly referring to the steady-state error of the thrust, is analyzed. As can be seen from the previous simulation, EGT is most

JID:AESCTE

AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.10 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

10

Fig. 8. Comparison of performance tracking with and without penalty factor in cost function.

likely to exceed the limitation among the three limit outputs. To maintain the actual EGT at EGT max , a steady state error in thrust is often caused. Fig. 9 presents maximum thrust level that the gas turbine engine can achieve with different EGT limits, where (a) is the response of thrust and (b) the response of EGT. The reference thrust is 7% and the upper limits of EGT are 1580K, 1600K, 1650K respectively. As can be seen from Fig. 9, thrust can be increased by 7% when the EGTmax is set as 1650K. When the EGTmax is limited to 1580K and 1600K, the engine runs close to the limit line of EGT, and the controller can only increase the thrust by 3% and 6% approximately, which means the limitation determine the maximum thrust that engine can provide. 5.3. The influence of controller parameters The prediction domain N y aims to optimize the cost function of the outputs in the future N y step size. When N y approaches infinity, it is equivalent to the global optimization algorithm, and the optimization can be completed offline all at once. Increasing N y can enhance the stability of the control system, but the dy-

namic response will slow down. The control domain N u should be less than or equal to N y and has the opposite effect on the system to N y . The influence of N u and N y on control quality and computational burden is studied in this section. 25 pairs of N u and N y are selected to study the influence of controller parameters from the following four aspects: calculation burden (characterized by computer running time), root mean square error (RMSE) between actual and reference thrust, overshoot and setting time (T s ). The simulation lasts for 10 seconds, and mainly assess the step response of the controller. When the system have run for 4 seconds, the thrust command will be increased by 5% while the S M Comp command keeps unchanged. All the data are obtained by averaging the results of 10 independent experiments. The computer running time, RMSE, overshoot and setting time of the simulation are listed in Table 2. To clearly reveal the change trend of these indexes, these discrete points are processed by a three-dimensional polynomial fitting mothed to get the curved surface graphs, as shown in Fig. 10. In Fig. 10, (a), (b), (c), (d) indicate computer running time, RMSE, overshoot and T s , respectively. The length of N u and N y directly determines the number of decision variables in QPSO optimization problem, which is related

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.11 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

11

Fig. 9. The influence of limitation on the steady-state error of the thrust. Table 2 Step response performance under different controller parameters. (N u , N y )

T /s

RMSE

Overshoot

T s /s

(N u , N y )

T /s

RMSE

Overshoot

T s /s

(1, 2)

1.8889

0.0068

0.2066

1.6480

(3, 9)

3.4013

0.0039

0.1155

0.704

(1, 3)

2.2802

0.0045

0.1258

1.1700

(3, 10)

3.4496

0.0038

0.1025

0.6622

(1, 4)

2.5796

0.0042

0.1389

0.7180

(4, 5)

2.8251

0.0042

0.1364

0.916 0.762

(1, 5)

2.6552

0.0039

0.1134

0.6600

(4, 6)

2.8465

0.0040

0.1091

(1, 6)

2.6204

0.0034

0.0938

0.6040

(4, 8)

2.8762

0.0040

0.1084

0.73

(2, 3)

2.1168

0.0047

0.1449

0.974

(4, 9)

3.4596

0.0040

0.1161

0.718

(2, 4)

2.2614

0.0044

0.1424

0.832

(4, 10)

3.5064

0.0039

0.1136

0.734

(2, 5)

2.4132

0.0042

0.1252

0.764

(5, 6)

2.8882

0.0040

0.1058

1.012

(2, 6)

3.0088

0.0041

0.0975

0.634

(5, 7)

2.9498

0.0041

0.1001

0.736

(2, 8)

3.1193

0.0036

0.1051

0.648

(5, 8)

3.0031

0.0041

0.1279

0.7911

(3, 5)

2.6457

0.0041

0.1289

0.828

(5, 9)

3.5630

0.0041

0.1252

0.74

(3, 6)

2.7490

0.0040

0.1037

0.764

(5, 10)

3.8771

0.0041

0.1142

0.77

(3, 7)

2.9351

0.0041

0.1117

0.756

to the computation burden. It can be analyzed from Table 2 and Fig. 10 (a) that, the computational burden of MGPC increases with the time domain; however, the quality of control will not always improve. When N u and N y increase to a certain value, the decline of RMSE, overshoot and setting time all become slow. In addition, due to the increasing complexity of the system, the above indicators may become worse in some cases. Synthetically considering the computational burden of the system and the requirements of the controller, N u = 3 and N y = 5 are selected for the other simulation in this paper. 5.4. Fault tolerant control In the actual control of gas turbine engines, it is usually necessary to sacrifice some performance in order to ensure their safety and stability after the failures occur. However, in some cases, engines still need to maintain the original operation state, which requires that controller has the ability of fault tolerance. Therefore, to make the control system designed in this paper have a wider range of applications, it is necessary to investigate whether it can guide the engine to operate efficiently and safely when and after the failures occur. The realization of fault-tolerant control requires the establishment of a fault diagnosis unit that can detect faults quickly and accurately, which is performed by UKF. Failures on high-pressure

shaft components and low-pressure shaft components have different effects on engines’ performance. In this paper, three typical failure cases are selected for simulation to study the practicability of the controller under different commands. They are the abrupt faults of the efficiency of the fan, the compressor and the HPT. As for other failures, these three cases can be referred to for treatment: the flow capacity faults and efficiency faults of the same component are similar, and faults in the LPT are similar to that in the HPT. According to the Ref. [40], the failure level can be classified as Table 3. The deterioration of the performance parameters of the engine is not obvious under mild faults; under severe faults, the safe operation of the engine may not be guaranteed and maintenance may be required. Therefore, the fault scenario is set according to the standard of medium faults in the simulation, that is, the health parameter decreases 3%. After the failure occurs, the assumptions of command change is mainly refers to Ref. [31]. There will be little effect on the operation of the core engine when the fan fails, and the requirements of the increasing thrust can be met within a reasonable range. The turbine is prone to overtemperature after failure, thus the engine needs to operate in a lower power state. The following simulation assumptions are made: 1. When the compressor fails, the reference will stay unchanged, 2. When the

JID:AESCTE

AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.12 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

12

Fig. 10. The curved surface fitted by the discrete points. Table 3 Gas path components failure classification. Failure

Mild

Medium

Severe

Fan or compressor efficiency Fan or compressor flow capacity Turbine efficiency Turbine flow capacity

−1.50% −1.50% −1.50% +1.50%

−3% −3% −3% +3%

−5% −5% −5% +5%

fan fails, the reference can be increased. 3. When the HPT fails, the reference will be reduced. It should be pointed that, these assumptions are only to verify the generality of the control method. The strategies and schemes after failure should be determined according to the specific flight mission requirements and failure conditions in practical applications. Example 1. Supposing the gas turbine engine is operating at the design point. At 6s, an abrupt medium fault occurs on the fan efficiency (S E 1 decrease 3%); at 10s, the thrust reference is raised by 2% and at 20s, the reference resumes. The dynamic response of the system under fault conditions is investigated. Fig. 11 (a) to

(d) display the health parameters estimated by the diagnostic system, the change of performance parameters, control variables and the limit outputs, respectively. In Fig. 11, the x-axes are the scales of the simulation time. The y-legend in figure (a) marks the absolute value of health parameters. In figure (b), (c), (d), y-axes of the thrust, W f , A 8 , N L and N H mean the normalized value of these parameters relative to those of the design point; y-axes of S M Comp and EGT mark their actual values. After the abrupt failure occurs on S E 1 during the steady running of the engine, the thrust decreases by about 1.5% without fault tolerant strategies, and S M Comp changes little. Under the faulttolerant control, the controller maintains the thrust at the pre-fault level when the command has not changed. At 10s, the control variables change rapidly in the initial phase of the response to increase the thrust; later, due to the limitation of the EGT, the response speed slows down lightly; however, desired thrust is achieved within 1s and there is no steady-state error, which ensures not only the speed of dynamic response, but also the safety of the engine. The program running times of fault diagnosis and controller calculation are 11.813s and 13.457s respectively.

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.13 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

13

Fig. 11. Simulation results of fault tolerant control with a failure on the fan’s efficiency.

Example 2. Supposing the gas turbine engine is operating at the design point. At 6s, an abrupt medium fault occurs on the compressor efficiency (S E 2 decrease 3%). The controller is required to restore the thrust at the pre-fault level. The simulation results are shown in Fig. 12. The definition of the legends is as above. It can be seen from Fig. 12, the thrust decreases by about 1.4% and the S M Comp decreases to 10% without fault tolerant strategies when the fault occurs. Under fault-tolerant control, the nonlinear model and the prediction model are adjusted in real time to adapt to the current engine health condition; thus the controller can quickly adjust the control variables and restore the thrust to the pre-fault level. The two controlled variables cannot be fully recovered at the same time due to the failures. The users can adjust the weight of S M Comp in the cost function, namely λ in Equation (17), according to their demands to select the variable to be recovered. If the thrust is required to restore as far as possible, λ can be chosen to be smaller or the opposite if it is for safety. In this example, the primary target is to restore the thrust, so λ is set to 0.1. We only need to ensure that the S M Comp will not fall to the surge boundary. Therefore, the S M Comp does not track the reference values, but decreases slightly in Fig. 12(b). In addition, the effect of the penalty factor in the cost function keeps the EGT under 1600K throughout the regulation process. The program running times of

fault diagnosis and controller calculation are 14.845s and 12.933s respectively. Example 3. Supposing the gas turbine engine is operating at the design point. At 6s, an abrupt medium fault occurs on the HPT efficiency (S E 3 decrease 3%). Since EGT is extremely prone to be overtemperature after the HPT fails, the thrust reference is reduced by 3% artificially to ensure safety. Fig. 13 provides the simulation results. The definition of the legends is as above. It should be pointed that in this example, the reference is modified only for verification. In practical applications, the pilots need should adjust the reference thrust depending on the current fight conditions and failure conditions. Fig. 13 shows that the S M Comp decreases to 1.4% without fault tolerant strategies when the fault occurs on the HPT. Although the thrust decreases by only 1%, the EGT increases to 1593K and is about to reach the upper bound. As can be seen from the simulation results, after modifying the thrust reference, the controller responses rapidly, and the thrust tracks the reference completely within 1s without steady-state error. However, EGT has not decreased significantly, indicating that the overtemperature caused by turbine component failure cannot be compensated by modifying the thrust instructions. And due to the action of the temperature

JID:AESCTE

AID:105576 /FLA

14

[m5G; v1.261; Prn:27/11/2019; 15:03] P.14 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

Fig. 12. Simulation results of fault tolerant control with a failure on the compressor’s efficiency.

limiter, the EGT is kept below the limit line, which ensures the safety working. In addition, modifying the thrust reference makes the S M Comp be partially restored, reaching 11.5%. The program running times of fault diagnosis and controller calculation are 10.770s and 16.849s respectively. 5.5. Robustness evaluations The robustness of the MGPC controller mainly refers to the ability of the system to remain within the allowable range when there is a mismatch between the prediction model and the actual process. At present, the robustness analysis of model predictive control is mostly illustrated by simulation results, which is mainly due to the difficulty of theoretical research on the algorithm construction and control structure of MPC [41]. In this section, some simulation will be carried out to evaluate the robustness of the controller under the conditions of measurement noise and modeling uncertainty. It should be pointed that all the analyses and simulation conducted in this paper have taken the measurement noise into account, and according to relevant research, the measurement noise can be approximated as Gaussian noise [8]. The noise processing is implemented by the nonlinear filter UKF and a reasonable covariance matrix can smooth the noise and help to estimate health

parameters accurately. In the following simulation, the noise covariance matrix R of UKF is fixed as R = 0.00152 and the noise level of the actual sensors will be continuously increased. Supposing that the standard deviations of sensor noise are 0.001, 0.0015, 0.003, and 0.005 in turn, the step response performance of the system is investigated. The simulation results are shown in Fig. 14. As can be seen from Fig. 14, the differences between the thrust and S M Comp responses are not obvious when disturbed by noise of different amplitudes, even if R is fixed. The thrust can quickly reach the reference value within 1s with almost no overshoot, while S M Comp remains unchanged. Only when the noise standard deviation increases to 0.005, the response fluctuation starts to become significant but can still remain near the reference value, indicating that the controller has a certain anti-noise capability. The uncertainty mainly comes from predictive CARIMA model. In this paper, CARIMA model is identified in real-time using the RLS. The accuracy of model identification has been verified in Section 5.1, and the anti-interference capability of the controller will be further investigated here. Supposing the gas turbine engine is operating at the design point. During 4s to 8s, continuous disturbances of different amplitudes are applied to the parameters of the CARIMA model. The elements of matrix polynomial A ( z−1 ) of the identified CARIMA model are multiplied by 200% and 500% respectively to see whether the system can be maintained at the original

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.15 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

15

Fig. 13. Simulation results of fault tolerant control with a failure on the HPT’s efficiency.

minimum. When  = 500%, the cost function increases obviously but its amplitude is still around 10−3 . The response curve shows that the outputs can maintain the original state except the fluctuation is more significant. The system is not out of balance due to model disturbance, which indicates that the controller has the ability of anti-interference. 6. Conclusion

Fig. 14. Step response at different noise levels.

state. Fig. 15 shows the simulation results, where (a) is the outputs curves and (b) the cost function. As can be seen from Fig. 15, when  = 200%, the outputs can follow the reference closely and the cost function is also kept at a

In this study, we proposed an improved MGPC algorithm which has been applied to realize the direct performance control of gas turbine engine. The unmeasurable performance parameters are estimated by the adaptive CLM and the parameters of the predictive model are online identified using the RLS algorithm with forgetting factors. In order to meet the requirements of limit protection, a penalty factor is added to the cost function. The nonlinear constraint optimization problem at each step is solved using the diversity-guided QPSO algorithm which is improved from two aspects of global average best position and initial global best position. The experiments proved the effectiveness of the proposed controller in terms of direct performance tracking, limit protection, fault tolerant control, as well as robustness against measurement noise and model uncertainties. The results are quite promising and can be a subject of interesting investigations. In the follow-

JID:AESCTE

AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.16 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

16

Fig. 15. Response under different disturbances.

ing research, further improvement of the operation speed of the controller will become a focus, especially under fault and/or degeneration conditions. In addition, since direct performance control has been realized, the proposed controller may also be applied to performance seeking control of gas turbine engines under different flight missions. Declaration of competing interest There is no conflict of interest.

[10]

[11]

[12]

[13]

Acknowledgements [14]

This work is supported by National Science and Technology Major Project (2017-V-0004-0054), the Fundamental Research Funds for the Central Universities (NO. 56XAA18034) and Chinese Scholarship Council (20190630019).

[15]

[16]

References [1] J.S. Litt, T.S. Sowers, S. Garg, A Retro-Fit Control Architecture to Maintain Engine Performance With Usage, Report No.: NASA/TM-2007-214977, Glenn Research Center, National Aeronautics and Space Administration, Cleveland, Ohio, 2007. [2] J.S. Litt, K.I. Parker, S. Chatterjee, Adaptive Gas Turbine Engine Control for Deterioration Compensation Due to Aging, Report No.: NASA/TM-2003-212607, Glenn Research Center, National Aeronautics and Space Administration, Cleveland, Ohio, 2003. [3] J. Csank, R.D. May, J.S. Litt, T.H. Guo, Control design for a generic commercial aircraft engine, in: 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Joint Propulsion Conferences, 2010 Jul 25-28, Nashville, TN, AIAA, 2010, pp. 2010–6629. [4] M. Montazeri-Gh, A. Safari, Tuning of fuzzy fuel controller for aero-engine thrust regulation and safety considerations using genetic algorithm, Aerosp. Sci. Technol. 15 (3) (2011) 183–192, https://doi.org/10.1016/j.ast.2010.10.004. [5] M. Henriksson, C. Breitholtz, Estimation of thrust and mass flow in a jet engine, in: Proceedings of the 2004 IEEE International Conference on Control Applications, 2004 Sept 2-4, Taipei, Taiwan, P.R.C., vol. 1, IEEE, 2004, pp. 229–235. [6] Z. Gastineau, G. Happawana, Robust model-based control for jet engines, in: 34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Joint Propulsion Conferences, Cleveland, OH, U.S.A., AIAA, 1998, p. 3752. [7] S. Adibhatla, H. Brown, Z. Gastineau, S. Adibhatla, H. Brown, Z. Gastineau, Intelligent engine control (IEC), in: 28th Joint Propulsion Conference and Exhibit, Nashville, TN, U.S.A., AIAA, 1992, p. 3484. [8] J.S. Litt, J. Turso, N. Shah, et al., A demonstration of a retrofit architecture for intelligent control and diagnostics of a turbofan engine, in: Infotech@ Aerospace Conferences, 2005 Sept 26-29, Arlington, Virginia, AIAA, 2005, p. 6905. [9] J.S. Litt, T. Sowers, Evaluation of an outer loop retrofit architecture for intelligent turbofan engine thrust control, in: 42nd AIAA/ASME/SAE/ASEE Joint

[17]

[18]

[19]

[20]

[21] [22] [23]

[24]

[25]

[26]

Propulsion Conference & Exhibit, 2006 Jul 9-12, Sacramento, California, AIAA, 2006, p. 5103. D. Ring, Henriksson M. Thrust, Control for a turbofan engine using estimation, in: ASME Turbo Expo 2006: Power for Land, Sea, and Air, 2006 May 8-11, Barcelona, Spain, ASME, 2006, pp. 901–910. C. Chen, Y. Sui, J. Zhao, Coordinated switching control of thrust tracking and safety protection for aero-engines, in: Proceedings of the 33rd Chinese Control Conference, 2014 Jul 28-30, Nanjing, P.R.C., IEEE, 2014, pp. 4050–4055. E. Mohammadi, M. Montazeri-Gh, Active fault tolerant control with selfenrichment capability for gas turbine engines, Aerosp. Sci. Technol. 56 (2016) 70–89, https://doi.org/10.1016/j.ast.2016.07.003. K. Zhu, J. Zhao, G.M. Dimirovski, H ∞ tracking control for switched LPV systems with an application to aero-engines, IEEE/CAA J. Autom. Sin. 5 (3) (2018) 699–705, https://doi.org/10.1109/JAS.2016.7510025. F.Q. Nian, H. Ma, G.H. Yang, Model reference adaptive control for aeroengine, in: The 27th Chinese Control and Decision Conference, 2015 May 23-25, Qingdao, P.R.C., IEEE, 2015, pp. 6563–6568. L.L. Tang, J.Q. Huang, M.X. Pan, Switching LPV control with double-layer LPV model for aero-engines, Int. J. Turbo Jet-Engines 34 (4) (2017) 313–320, https:// doi.org/10.1515/tjj-2016-0014. B.J. Brunell, D.E. Viassolo, R. Prasanth, Model adaptation and nonlinear model predictive control of an aircraft engine, in: ASME Turbo Expo 2004: Power for Land, Sea, and Air, 2004 Jun 14-17, Vienna, Austria, ASME, 2004, pp. 673–682. X. Du, Y.Q. Guo, X.L. Chen, Multivariable constrained predictive control and its application to a commercial turbofan engine, Adv. Mater. Res. 909 (2014) 281–287, https://doi.org/10.4028/www.scientific.net/AMR.909.281. A.P. Wiese, M.J. Blom, C. Manzie, et al., Model reduction and MIMO model predictive control of gas turbine systems, Control Eng. Pract. 45 (2015) 194–206, https://doi.org/10.1016/j.conengprac.2015.09.015. K. Salahshoor, A. Khaki-Sedigh, P. Sarhadi, An indirect adaptive predictive control for the pitch channel autopilot of a flight system, Aerosp. Sci. Technol. 45 (2015) 78–87, https://doi.org/10.1016/j.ast.2015.04.016. Q.G. Zheng, Z.G. Xu, H.B. Zhang, et al., A turboshaft engine NMPC scheme for helicopter autorotation recovery maneuver, Aerosp. Sci. Technol. 76 (2018) 421–432, https://doi.org/10.1016/j.ast.2018.01.034. J. Richalet, Industrial applications of model based predictive control, Automatica 29 (5) (1993) 1251–1274, https://doi.org/10.1016/0005-1098(93)90049-Y. D.W. Clarke, C. Mohtadi, Properties of generalized predictive control, Automatica 25 (6) (1989) 859–875, https://doi.org/10.1016/0005-1098(89)90053-8. D.W. Clarke, C. Mohtadi, P.S. Tuffs, Generalized predictive control—Part I. The basic algorithm, Automatica 23 (2) (1987) 137–148, https://doi.org/10.1016/ 0005-1098(87)90087-2. M. Sedraoui, S. Abdelmalek, S. Gherbi, Multivariable generalized predictive control using an improved particle swarm optimization algorithm, Informatica 35 (3) (2011) 363–374. P. Alkorta, O. Barambones, J.A. Cortajarena, et al., Efficient multivariable generalized predictive control for sensorless induction motor drives, IEEE Trans. Ind. Electron. 61 (9) (2014) 5126–5134, https://doi.org/10.1109/TIE.2013.2281172. J. Smoczek, J. Szpytko, Particle swarm optimization-based multivariable generalized predictive control for an overhead crane, IEEE/ASME Trans. Mechatron. 22 (1) (2017) 258–268, https://doi.org/10.1109/TMECH.2016.2598606.

JID:AESCTE AID:105576 /FLA

[m5G; v1.261; Prn:27/11/2019; 15:03] P.17 (1-17)

X. Zhou et al. / Aerospace Science and Technology ••• (••••) ••••••

[27] Q. Ye, X. Lou, L. Sheng, Generalized predictive control of a class of MIMO models via a projection neural network, Neurocomputing 234 (2017) 192–197, https://doi.org/10.1016/j.neucom.2016.12.067. [28] J. Sun, W. Xu, B. Feng, A global search strategy of quantum-behaved particle swarm optimization, in: IEEE Conference on Cybernetics and Intelligent Systems, 2004 Dec 1-3, Singapore, Singapore, IEEE, 2004, pp. 111–116. [29] P. Codron, P. Boucher, Multivariable generalized predictive control with multiple reference models: a new choice for the reference model, in: ECC’93. European Control Conf., June 28-July, 1993. [30] P. Zelinka, B. RohaˇI-Ilkiv, A.G. Kuznetsov, Experimental verification of stabilising predictive control, Control Eng. Pract. 7 (5) (1999) 601–610, https://doi.org/10. 1016/S0967-0661(99)00019-2. [31] H. Yao, Full Authority Digital Electronic Control System for Aero-Engine, Aviation Industry Press, Beijing, 2014, pp. 218–223, Chinese. [32] W.X. Zhou, Research on Object-Oriented Modeling and Simulation for Aeroengine and Control System, Ph.D. thesis, Nanjing University of Aeronautics and Astronautics, Nanjing, China, 2006. [33] F. Lu, W. Zheng, J. Huang, M. Feng, Life cycle performance estimation and inflight health monitoring for gas turbine engine, J. Dyn. Syst. Meas. Control 138 (9) (2016) 091009, https://doi.org/10.1115/1.4033556. [34] S.J. Julier, J.K. Uhlmann, H.F. Durrant-Whyte, A new approach for filtering nonlinear systems, in: Proceedings of the American Control Conference, 1995 Jul 21-23, Seattle, Washington, IEEE, 1995, pp. 1628–1632.

17

[35] S. Julier, J. Uhlmann, H.F. Durrant-Whyte, A new method for the nonlinear transformation of means and covari-ances in filters and estimators, IEEE Trans. Autom. Control 45 (3) (2000) 477–482, https://doi.org/10.1109/9.847726. [36] Lu Feng, et al., Gas path on-line fault diagnostics using a nonlinear integrated model for gas turbine engines, Int. J. Turbo Jet-Engines 31 (3) (2014) 261–275, https://doi.org/10.1515/tjj-2014-0001. [37] M. Clerc, J. Kennedy, The particle swarm - explosion, stability, and convergence in a multidimensional complex space, IEEE Trans. Evol. Comput. 6 (1) (2002) 58–73, https://doi.org/10.1109/4235.985692. [38] J. Sun, W. Xu, W. Fang, A diversity-guided quantum-behaved particle swarm optimization algorithm, in: Asia-Pacific Conference on Simulated Evolution and Learning, Seal, Springer, 2006, pp. 497–504. [39] J. Sun, W. Xu, W. Fang, Quantum-Behaved Particle Swarm Optimization Algorithm With Controlled Diversity. International Conference on Computational Science, Springer, 2006, pp. 847–854. [40] A. Kumar, D. Viassolo, Model-Based Fault Tolerant Control, Report No.: NASA/TM-2008-215273, NASA, Washington, D.C., 2008. [41] J.A. DeCastro, Rate-based model predictive control of turbofan engine clearance, J. Propuls. Power 23 (4) (2007) 804–813, https://doi.org/10.2514/1.25846.