Aero-engine system identification and output interval estimation method comprehensively considering performance degradation and model mismatch

Abstract

The invention belongs to the technical field of model parameter identification, and discloses an aero-engine system identification and output interval estimation method comprehensively considering performance degradation and model mismatch. For an aero-engine system, a linear variable parameter modeling method is introduced, state variable models of multiple working points of an aero-engine are combined through scheduling parameters, and a large-working-range global model is constructed; aiming at the problems of errors, model mismatch and the like in the modeling process, interval estimation is carried out on the system output uncertainty based on a quantile regression neural network, and the model precision and robustness are improved; a standard state space equation of an aero-engine system is augmented, health parameters of all parts of the system are regarded as states to be estimated, the health parameters are estimated in real time through an adaptive Kalman filtering method based on variational Bayes, and a system adaptive model is constructed. The method provides a research basis for the fields of fault diagnosis, nonlinear system control and the like of an aero-engine system.

Description

Technical Field

[0001] The present invention relates to the technical field of model parameter identification, and in particular to an aero-engine system identification and output interval estimation method that comprehensively considers performance degradation and model mismatch. Background Art

[0002] As an important part of the aircraft, the aircraft engine shoulders the important responsibility of providing power for the aircraft. The performance of the aircraft depends largely on the performance of the aircraft engine. Therefore, the method of constructing an adaptive model of an aircraft engine with strong real-time performance, high accuracy and simple expression is a key issue in the aviation industry. Due to the strong nonlinearity of the aircraft engine system, it is difficult for the traditional linearization method to reflect the system characteristics within a large working range. The linear variable parameter modeling method combines the linearized state variable models at multiple stable working points during the engine operation process by using scheduling parameters, thereby improving the modeling accuracy of the system. The adaptive Kalman filter algorithm can estimate the system model and noise characteristics online through observation data. Compared with the traditional Kalman filter, it has the ability to adapt to the dynamic changes of the system, so it works better when dealing with nonlinear systems and variable parameter systems.

[0003] Aircraft engines provide power for aircraft flight and they mainly work in harsh environments of high temperature and high pressure. Due to the influence of the working environment and the wear of components during use, the actual engine performance degrades and the uncertainty is strong. At present, in the research of linear system modeling methods of aircraft engine systems, most of them can only obtain point estimates of model parameters, which is difficult to reflect the system performance degradation process and complex uncertainties. The interval prediction method is used to replace the point estimate with the distribution estimate of the output result to achieve the measurement of the uncertainty of the aircraft engine system output estimate.

[0004] The objective of the present invention is to design an adaptive aircraft engine system model and predict system health parameters and performance degradation, which is achieved by combining a linear variable parameter model with an adaptive Kalman filter algorithm. Summary of the invention

[0005] The purpose of the present invention is to solve the problem of lack of a modeling of the system's large operating range characteristics when noise and performance degradation exist in aircraft engine system modeling, and to propose an aircraft engine system identification and output range estimation method that comprehensively considers performance degradation and model mismatch.

[0006] The technical solution of the present invention is as follows: a method for aircraft engine system identification and output interval estimation that comprehensively considers performance degradation and model mismatch, comprising the following steps:

[0007] Step 1: Based on the input, output and state variables of the aircraft engine system, a local linearization model at each operating point of the aircraft engine is established based on the perturbation method, and a global model of the aircraft engine in a large operating range is established based on the polynomial linear variable parameter method;

[0008] Step 2: Based on the adaptive Kalman filter theory, the performance degradation of each component of the aircraft engine is estimated in the presence of external noise disturbances based on the global model of the large working range in step 1;

[0009] Step 3: Based on the long short-term memory neural network architecture and combined with the quantile regression method, the residual between the model prediction results obtained in step 1 above and the actual aircraft engine model data is compensated:

[0010] The final model output is,

[0011] y c (n L ,u)=y c,st (n L ,u c,st )+△y c (n L ,△u c )+y pred (n L ,u c,st );

[0012] Among them, y c (n L ,u) is the final predicted value of the system output, y c,st (n L ,u c,st ) represents the estimate of the current equilibrium point output, △y c (n L ,△u c ) is the output increment estimate when the input changes, y pred (n L ,u c,st ) is the compensation value of the neural network to the global model.

[0013] The step 1 is specifically as follows:

[0014] Step 1.1: Set the nonlinear aerodynamic thermodynamic model of the aircraft engine system as follows:

[0015]

[0016] y(t)=g(x(t),u(t))

[0017] Where x∈R n is the state quantity of the aircraft engine system, y∈R mis the output of the aircraft engine system, u∈R r is the input of the aircraft engine system, while f(·) and g(·) represent the nonlinear functions in the state equation and output equation of the aircraft engine model, respectively, which are used to describe the dynamic process of engine operation and the output observation results;

[0018] At a certain equilibrium point (x st ,u st ) is used to perform Taylor approximation on the above nonlinear aerodynamic thermodynamic model, and its first-order expansion is as follows:

[0019]

[0020]

[0021] Right now,

[0022]

[0023] △y=C△x+D△u

[0024] Where A, B, C, D are system matrices of corresponding dimensions, △x=xx st Represents the increment of the true state value relative to the equilibrium state value, △u=uu st Represents the increment of the real state input relative to the equilibrium point input, △y=yy st It represents the increment of the actual state output relative to the equilibrium point output;

[0025] The above system matrix is ​​solved based on the local perturbation method and the step signal excitation method. First, when the working state of the system is at the equilibrium point, a disturbance signal is given to a state variable of the system, and the system matrix A={a ij},C={c oj} to solve;

[0026]

[0027] Among them, △x j Represents the jth state variable x j The disturbance applied is specified, and the disturbance value cannot exceed one percent of the state value of the selected equilibrium point; similarly, the elements of the matrix C are also obtained by the above method, that is,

[0028]

[0029] By inputting a step signal to observe the change of a state variable to calculate the values ​​of the system matrix B and D; use the step response of the first control variable to calculate the elements of the first column of the system matrix, and keep other control variables and state variables unchanged, that is,

[0030]

[0031]

[0032] Among them, △u1 is the step amplitude of the first control variable, which is stipulated not to exceed 1% of the input value of the selected equilibrium point; similarly, the above method is used to calculate the element value of each column of matrices B and D;

[0033] Step 1.2: Based on the benchmark values ​​of multiple equilibrium points and the local linear state variable model, a dynamic model of the system under large working conditions is constructed based on the linear variable parameter method;

[0034] The state variable model calculated above can only effectively approximate the characteristics of the system in a small neighborhood near the selected equilibrium operating point. When the deviation from the operating point is too much, the state variable model is difficult to apply. In order to obtain the dynamic characteristics of the aircraft engine in a large operating range; based on the polynomial linear variable parameter method, according to the state variable models corresponding to multiple equilibrium points of the system under different working conditions, the variable parameter model of the system is constructed;

[0035] A scheduling parameter θ is introduced to reflect the current working state of the system, and the state variable model at each equilibrium point is regarded as an affine function of the linear variable parameter model with respect to the scheduling parameter, that is,

[0036]

[0037] △y c =C c (θ)△x c +D c (θ)△u c +ν

[0038] Each element in the system matrix is ​​a variable value of the scheduling parameter, that is, A c (θ) = {a c,ij (θ)}, where a c,ij Represents the matrix A c The elements in ω and ν represent the process noise and observation noise of the system respectively, both of which obey the zero-mean Gaussian distribution, and the variances are Q and R respectively; the low-pressure turbine speed n of the aircraft engine is selected L As the dispatch parameter of the linear variable parameter model;

[0039] In order to determine the benchmark values ​​of model input, output and state under different working conditions Wherein, k represents the model reference value of the kth operating point, and multiple m-order polynomials are used to approximate the relationship between the scheduling parameters and each reference value, as shown below;

[0040]

[0041] In the fitting process, the value of t is continuously increased until the coefficient of determination R of the fitting is 2 Stop when it is greater than 90%, t=min{t:R 2 (x c,st,n (n L,c ,t))>=0.9}; By calculating each state variable, output and output reference value through the above formula, the steady-state output of the system for the corresponding working state is finally obtained;

[0042] In order to express the dynamic characteristics of the system, each element in the state space system matrix of multiple working points is fitted to obtain a power series polynomial about the scheduling parameters. The specific expression is as follows:

[0043]

[0044]

[0045]

[0046]

[0047] Among them, ta, tb, tc, and td represent the polynomial orders of different system matrices, and the values ​​of each order are continuously increased during the fitting process until the determination coefficient R of the fitting is 2 When it is greater than 90%, it stops; finally, the adaptive dynamic model of the system under large working conditions is obtained. For the working point n L,1 The state value and output value of the system are obtained through the interaction of the corresponding reference value and the variable parameter state variable model;

[0048] y c (n L,1 )=y c,st (n L,1 ,u c,st )+△y c (n L,1 ,△u c )

[0049] Among them, △u c It represents the increment of the input of the system at this moment relative to the previous sampling moment. The right side of the equal sign represents the sum of the baseline value of the observed output corresponding to the operating point and the input increment response. Similarly, the dynamic characteristics of the state variables of the system at different operating points are represented by a similar method.

[0050] The step 2 is specifically as follows: in order to estimate the system performance degradation, a health parameter is introduced to represent the efficiency and flow deviation of the engine rotating parts; the health parameter h, which reflects the performance degradation of each engine component, is augmented to the state vector of the original system, that is, Since engine performance degradation is a slow process, the partial derivative of the performance degradation factor is approximated to 0; the augmented model is

[0051]

[0052]

[0053] Right now,

[0054]

[0055] △y ext =C ext (θ)x ext +D ext (θ)△u c +ν

[0056] Among them A h (θ),B c (θ),C h (θ) is obtained by step 1, considering the change of estimated health parameters at the equilibrium point and ignoring △u c item;

[0057] Considering the state variable model at a certain working point, the linear time-invariant model is first discretized to obtain

[0058] x k =A k x k-1 +ω k

[0059] y k =H k x k-1 +ν k

[0060] Among them, ω k ~N(0,Q k ) is the system noise that follows Gaussian distribution, and ν k ~N(0,Σ k ) is the observation noise, the noise variance Σ k is a diagonal matrix; assuming that the initial state value obeys Gaussian distribution, that is, x0~N(m0,P0), the subscript k represents the sampling time;

[0061] The working process of the adaptive Kalman filter mainly includes the following steps:

[0062] Calculate the forecast distribution parameters:

[0063]

[0064]

[0065]

[0066]

[0067] Among them, α, β are parameters in the inverse gamma distribution, which are used to describe the variance of the observation noise, and ρ∈(0,1] is an adjustable constant;

[0068] Parameter update step: First, let Then the following steps are iteratively solved;

[0069]

[0070]

[0071]

[0072]

[0073] After N iterations, let Finally, the state distribution of the system at the working point is obtained, that is, the predicted value of the health parameter;

[0074] Use the current working state of the system and the scheduling parameter θ as the input x, h of the long short-term memory neural network n Represents the network output; first, the input is saved to the unit state c through the input gate t middle;

[0075]

[0076] Among them, W i is the input gate weight matrix, b i is the bias term, Indicates the output status at the previous moment; the unit state of the current input of the system is,

[0077]

[0078] After that, the forget gate determines whether to keep the cell state c at the previous moment. t-1 ;

[0079]

[0080] Determine the cell state of the LSTM neural network at the current moment;

[0081]

[0082] in Represents the Hadamard product; according to the current unit state, the final output of the long short-term memory neural network is obtained through the output gate;

[0083]

[0084]

[0085] Based on the quantile regression method, a fully connected layer is used to obtain the interval estimation of the system residual; the quantile level q∈Q={0.05,0.5,0.95} is taken to obtain the prediction results under different quantiles. At this time, the loss function of quantile regression is

[0086]

[0087] Where (y i ) + =max(y i ,0), the output when the quantile is 0.5 Represents the residual compensation term for the LPV system, that is, The prediction results when taking the other two quantiles represent interval predictions with 90% confidence.

[0088] Beneficial effects of the present invention: The present invention is aimed at the aero-engine system, and for the research field of aero-engine system identification, the problem of poor modeling accuracy of nonlinear systems in a large working range is proposed. A method for aero-engine system identification and output interval estimation that comprehensively considers performance degradation and model mismatch is proposed, and an adaptive model of multiple working points of the aero-engine system is obtained. At the same time, the performance degradation and noise interference of each system component are considered to estimate the health parameters of the system. On the one hand, the present invention introduces the linear variable parameter model method into the aero-engine system modeling process, and combines the state variable models of different steady-state working points of the system by constructing affine functions of multiple scheduling parameters to construct a system model in a large working range; on the other hand, the health parameters of each component of the system are augmented into the state vector, and the state equation of the aero-engine system is augmented. The health parameters of the system are estimated in real time through an adaptive Kalman filtering method based on the variational Bayesian idea, and the performance degradation of the aero-engine system is estimated according to the changes in the health parameters. BRIEF DESCRIPTION OF THE DRAWINGS

[0089] Figure 1 The schematic diagram shows a method for aircraft engine system identification and output interval estimation that comprehensively considers performance degradation and model mismatch.

[0090] Figure 2 Schematic diagram of input disturbance signal.

[0091] Figure 3 Schematic diagram of residual adaptation for the neural network compensation model.

[0092] Figure 4 Output results for the low-pressure turbine speed of the system model with interval prediction.

[0093] Figure 5 Output results for the high-pressure turbine speed of the system model with interval prediction.

[0094] Figure 6 This is the turbine outlet temperature output result of the system model after model mismatch compensation.

[0095] Figure 7 It is the estimation result of the health parameters of each component by the adaptive Kalman filter. DETAILED DESCRIPTION

[0096] A method for aircraft engine system identification and output interval estimation that comprehensively considers performance degradation and model mismatch is shown in the flowchart below. Figure 1 As shown, it is implemented by the following steps:

[0097] Step 1: Based on the input, output and state variables of the aircraft engine system, a local linearization model at each operating point of the aircraft engine is established based on the perturbation method, and a global model of the aircraft engine in a large operating range is established based on the polynomial LPV method;

[0098] Step 1.1: The nonlinear aerodynamic thermodynamic model of the aero-engine system is as follows:

[0099]

[0100] y(t)=g(x(t),u(t))

[0101] Where x∈R n is the state quantity of the aircraft engine system, y∈R m is the observation quantity of the aero-engine system, u∈R r is the input of the aircraft engine system, while f(·) and g(·) represent the nonlinear functions in the state equation and output equation of the aircraft engine model, respectively, which are used to describe the dynamic process and observation results of the engine operation. At a certain equilibrium operating point (x st ,u st ) is used to perform Taylor approximation on the above nonlinear model, and its first-order expansion is as follows:

[0102]

[0103]

[0104] Right now,

[0105]

[0106] △y=C△x+D△u

[0107] Where A, B, C, D are system matrices of appropriate dimensions, △x = xx st Represents the increment of the true state value relative to the equilibrium state value, △u=uu st and △y=yy st Similarly, it represents the increment of system input and output.

[0108] The above system matrix is ​​solved based on the local perturbation method and the step signal excitation method. First, a disturbance signal is given to a state variable of the system when the system is in the equilibrium state, and the system matrix A = {a ij},C={c mj} to solve.

[0109]

[0110] Among them, △x j Represents the jth state variable x j The disturbance is applied, and the disturbance value is stipulated not to exceed one percent of the state value of the selected equilibrium point. Similarly, the elements of the matrix C can also be obtained by the above method, that is,

[0111]

[0112] By inputting a step signal, we observe the change of a certain state variable to calculate the values ​​of system matrices B and D. We use the step response of the first control variable to calculate the elements of the first column of the system matrix, and keep other control variables and state variables unchanged, that is,

[0113]

[0114]

[0115] Where △u1 is the step amplitude of the first control variable, and it is stipulated that the step value cannot exceed one percent of the input value of the selected equilibrium point. Similarly, the above method can be used to calculate the element value of each column of matrices B and D.

[0116] Step 1.2: Based on the benchmark values ​​of multiple equilibrium points and the state variable model, a dynamic model of the system under large working conditions is constructed based on the linear variable parameter method;

[0117] The state variable model calculated above can only effectively approximate the characteristics of the system in a small neighborhood near the selected equilibrium operating point. When the deviation from the operating point is too much, the state variable model is difficult to apply. In order to obtain the dynamic characteristics of the aircraft engine in a large operating range, based on the polynomial linear variable parameter method, according to the state variable models corresponding to multiple equilibrium points of the system under different working conditions, an adaptive model of the system is constructed.

[0118] A scheduling parameter θ that can reflect the current working state of the system is introduced into the system model, and the state variable model at each equilibrium point is regarded as an affine function of the linear variable parameter model with respect to the scheduling parameter, that is,

[0119]

[0120] △y c =C c (θ)△x c +D c (θ)△u c +ν

[0121] Each element in the system matrix is ​​a variable value of the scheduling parameter, that is, A c (θ) = {a c,ij (θ)}. ω and ν represent the process noise and observation noise of the system respectively, both of which obey zero-mean Gaussian distribution, and their variances are Q and R respectively. The low-pressure turbine speed n of the aircraft engine is selected L As the scheduling parameter of the linear variable parameter model.

[0122] In order to determine the benchmark values ​​of model input, output and state under different working conditions Where k represents the model reference value of the kth operating point, and multiple m-order polynomials are used to approximate the relationship between the scheduling parameters and each reference value, as shown below

[0123]

[0124] The fitting process is specified to continuously increase the value of t until the coefficient of determination R 2 Stop when it is greater than 90%, that is, t=min{t:R 2 (x c,st,n (n L,c ,t))>=0.9}. Similarly, each state variable, output and output reference value can be calculated through the above formula, and finally the steady-state output of the system for the corresponding working state can be obtained.

[0125] In order to express the dynamic characteristics of the system, each element in the state space system matrix of multiple working points is fitted to obtain a power series polynomial about the scheduling parameters. The specific expression is as follows:

[0126]

[0127]

[0128]

[0129]

[0130] Among them, ta, tb, tc, and td represent the polynomial orders of different system matrices, and the values ​​of each order are continuously increased during the fitting process until the determination coefficient R of the fitting is 2 When it is greater than 90%, it stops; finally, the adaptive dynamic model of the system under large working conditions is obtained. For the working point n L,1 The state value and output value of the system can be obtained by the corresponding reference value and the variable parameter state variable model.

[0131] y c (n L,1 )=y c,st (n L,1 ,u c,st )+△y c (n L,1 ,△u c )

[0132] Where △u c It represents the increment of the input of the system at this moment relative to the previous sampling moment, and the right side of the equal sign represents the sum of the baseline value of the observed output corresponding to the working point and the input increment response. Similarly, the dynamic characteristics of the state variables of different working points of the system can also be expressed by a similar method.

[0133] Step 2: Based on the adaptive Kalman filter theory, the performance degradation of each component of the aircraft engine is estimated, and the influence of external noise on the model built in step 1 is eliminated;

[0134] In order to estimate the system performance degradation, the health parameter is introduced to represent the efficiency and flow deviation of the engine rotating parts. The health parameter h is expanded to the state vector of the original system Since engine performance degradation is a slow process, the partial derivative of the performance degradation factor can be approximated to 0. The augmented model is

[0135]

[0136]

[0137] Right now

[0138]

[0139] △yext =C ext (θ)x ext +D ext (θ)△u c +ν

[0140] Among them A h (θ),B c (θ),C h (θ),D h (θ) can be obtained by step 1, considering the change of estimated health parameters at the equilibrium point and ignoring △u c item.

[0141] Considering the state variable model at a certain operating point (the other operating points can be obtained by different scheduling parameters), the linear time-invariant model is first discretized to obtain

[0142] x k =A k x k-1 +ω k

[0143] y k =H k x k-1 +ν k

[0144] Among them, ω k ~N(0,Q k ) is the system noise that follows Gaussian distribution, and ν k ~N(0,Σ k ) is the observation noise, the noise variance Σ k It is a diagonal matrix. Assume that the initial state value follows Gaussian distribution, that is, x0~N(m0,P0), and the subscript k represents the sampling time.

[0145] The working process of the adaptive Kalman filter mainly includes the following steps:

[0146] Calculate the forecast distribution parameters:

[0147]

[0148]

[0149]

[0150]

[0151] Among them, α and β are parameters in the inverse gamma distribution, which are used to describe the variance of the observation noise, and ρ∈(0,1] is an adjustable constant.

[0152] Parameter update step: First, let Then iterate and solve the following steps

[0153]

[0154]

[0155]

[0156]

[0157] After N iterations, let Finally, the state distribution of the system at the working point is obtained, that is, the predicted value of the health parameter.

[0158] Step 3: Based on the long short-term memory neural network architecture and combined with the quantile regression method, the residuals between the model prediction results obtained in the above steps and the actual aircraft engine model data are compensated:

[0159] Use the current working state of the system, that is, the scheduling parameter θ, as the input x of the long short-term memory neural network. First, save the input to the unit state c through the input gate. t middle

[0160] i t =σ(W i ·[h t-1 ,x t ]+b i )

[0161] Among them, W i is the input gate weight matrix, b i is the bias term, h t-1 Indicates the output status at the last moment.

[0162] The current input unit state of the system is

[0163]

[0164] After that, the forget gate determines whether to keep the cell state c at the previous moment. t-1

[0165] f t =σ(W f ·[h t-1 ,x t ]+b f )

[0166] Then the unit state of the network at the current moment can be determined

[0167]

[0168] in Represents the Hadamard product. According to the current unit state, the final output of the network can be obtained through the output gate;

[0169] o t =σ(W o ·[h t-1 ,x t ]+b o )

[0170]

[0171] Based on the quantile regression method, a fully connected layer is used to obtain the interval estimation of the system residual. The quantile level q∈Q={0.05,0.5,0.95} is taken to obtain the prediction results under different quantiles. At this time, the loss function of quantile regression is

[0172]

[0173] Where (y i ) + =max(y i ,0), that is, only positive numbers are retained, and the output when the quantile is 0.5 Represents the residual compensation term for the LPV system, that is, The prediction results for the other two quantiles represent interval predictions with a 90% confidence level. The final model output is

[0174] y c (n L ,u)=y c,st (n L ,u c,st )+△y c (n L ,△u c )+y pred (n L ,u c,st )

[0175] The following is a specific embodiment:

[0176] The Simulink simulation model of the turbofan engine JT9D developed by NASA is selected to obtain the data required for the aircraft engine system identification. The state variable x is set to the low-pressure rotor speed and high-pressure rotor speed of the aircraft engine, that is, x = [n l ,n h ] T, the output variable y is set as the high-pressure turbine outlet temperature of the aircraft engine, and the input variable u is set as the input oil-air ratio FAR of the aircraft engine. A Gaussian noise is added to the output signal as the observation noise of the system, and the aircraft engine modeling is performed based on the operating data in the JT9D engine Simulink simulation model.

[0177] Near the steady-state operating point of the engine, the small disturbance method is used to excite the system, and the partial derivative method is used to obtain the approximate linear state variable model of the system at this operating point through the input-output characteristics of the system. The disturbance input curve is as follows Figure 2 .

[0178] With the scheduling parameter θ0(n l =3500rpm) as an example, the linear model of this working point is estimated by the small disturbance method as follows:

[0179]

[0180] △y=C△x+D△u

[0181] in:

[0182]

[0183] C=[0.4344,0.0064],D=[38092.33]

[0184] When the aircraft engine is in the dispatching parameter range of 80% to 100% (θ 100％ =nl 100％ =4000rpm), 17 steady-state operating points (interval of 50rpm) were selected from the working range, and the corresponding state variable models were obtained using the perturbation method. Considering that the scheduling parameters are different from the order of magnitude of the system input and output, the scheduling parameters are normalized, and then the model matrix elements corresponding to these operating points are fitted using the fitting method to obtain an affine function form about the scheduling parameters, and then the linear variable parameter model of the engine system in the working range of 80% to 100% of the scheduling parameters is obtained:

[0185]

[0186] △y c =C c (θ)△x c +D c (θ)△u c +ν

[0187] The first element a of matrix A 11 Taking as an example, the function of the scheduling parameters obtained is,

[0188] a11 (θ) = -129.78 × θ 5 +301.98×θ 4 -250.03×θ 3 +80.77×θ 2 -7.33×θ-1.09

[0189] The index determination coefficient of the fitting accuracy is R2=0.9881, and the fitting accuracy of the remaining elements of the system matrix is ​​shown in Table 1 below.

[0190] Table 1 Determination coefficient of each element in the system matrix

[0191]

[0192] It can be seen from the index data in Table 1 that the polynomial fitting accuracy of the system matrix elements with respect to the scheduling parameters meets the index requirements, and the obtained system adaptive model has high accuracy under the selected working conditions.

[0193] In order to further improve the accuracy of the model, a neural network combining long short-term memory network and quantile regression method is used to compensate for the residual of the model. The system scheduling parameters are used as the input of the network, and the residual term caused by model mismatch is used as the network output to train the network. The interval estimation result of the residual term is obtained as follows: Figure 3 After the network training with quantile q = 0.5 is completed, it is imported into Simulink to perform online compensation on the system model. The output and state change of the system adaptive model are compared with the data of the component-level model JT9D. Figures 4 to 6 .

[0194] Compared with the traditional Kalman filtering method, the adaptive Kalman filtering method combined with the variational Bayesian theory has also improved the identification accuracy to a certain extent. The health parameters of each component are augmented into the state vector to obtain a new system model. The prediction of the health parameters of each component of the system based on the adaptive Kalman filtering algorithm is shown in Figure 7 The estimation accuracy comparison between the adaptive Kalman filter algorithm and the ordinary Kalman filter algorithm is shown in Table 2.

[0195] Table 2 Method comparison diagram

[0196]

[0197] It can be seen that compared with the Kalman filter method, the accuracy of output estimation is also improved.

Claims

1. A method for aircraft engine system identification and output interval estimation that comprehensively considers performance degradation and model mismatch, characterized in that: The steps include: Step 1: Based on the input, output and state variables of the aircraft engine system, a local linearization model at each operating point of the aircraft engine is established based on the perturbation method, and a global model of the aircraft engine over a large operating range is established based on the polynomial linear variable parameter method; Step 2: Based on the adaptive Kalman filter theory, the performance degradation of each component of the aircraft engine is estimated in the presence of external noise disturbances based on the global model of the large working range in step 1; Step 3: Based on the long short-term memory neural network architecture and combined with the quantile regression method, the residual between the model prediction results obtained in step 1 above and the actual aircraft engine model data is compensated: The final model output is, y c (n L ,u)=y c,st (n L ,u c,st )+△y c (n L ,△u c )+y pred (n L ,u c,st ); Among them, y c (n L ,u) is the final predicted value of the system output, y c,st (n L ,u c,st ) represents the estimate of the current equilibrium point output, △y c (n L ,△u c ) is the output increment estimate when the input changes, y pred (n L ,u c,st ) is the compensation value of the neural network to the global model.

2. The method for aircraft engine system identification and output interval estimation taking into account performance degradation and model mismatch according to claim 1, characterized in that: The step 1 is specifically as follows: Step 1.1: Set the nonlinear aerodynamic thermodynamic model of the aircraft engine system as follows: y(t)=g(x(t),u(t)) Where x∈R n is the state quantity of the aircraft engine system, y∈R m is the output of the aircraft engine system, u∈R r is the input of the aircraft engine system, while f(·) and g(·) represent the nonlinear functions in the state equation and output equation of the aircraft engine model, respectively, which are used to describe the dynamic process of engine operation and the output observation results; At a certain equilibrium point (x st ,u st ) is used to perform Taylor approximation on the above nonlinear aerodynamic thermodynamic model, and its first-order expansion is as follows: Right now, △y=C△x+D△u Where A, B, C, D are system matrices of corresponding dimensions, △x=xx st Represents the increment of the true state value relative to the equilibrium state value, △u=uu st Represents the increment of the real state input relative to the equilibrium point input, △y=yy st It represents the increment of the actual state output relative to the equilibrium point output; The above system matrix is ​​solved based on the local perturbation method and the step signal excitation method. First, when the working state of the system is at the equilibrium point, a disturbance signal is given to a state variable of the system, and the system matrix A={a ij },C={c oj } to solve; Among them, △x j Represents the jth state variable x j The disturbance applied is specified, and the disturbance value cannot exceed one percent of the state value of the selected equilibrium point; similarly, the elements of the matrix C are also obtained by the above method, that is, By inputting a step signal to observe the change of a state variable to calculate the values ​​of the system matrix B and D; use the step response of the first control variable to calculate the elements of the first column of the system matrix, and keep other control variables and state variables unchanged, that is, Among them, △u1 is the step amplitude of the first control variable, which is stipulated not to exceed 1% of the input value of the selected equilibrium point; similarly, the above method is used to calculate the element value of each column of matrices B and D; Step 1.2: Based on the benchmark values ​​of multiple equilibrium points and the local linear state variable model, a dynamic model of the system under large working conditions is constructed based on the linear variable parameter method; Based on the polynomial linear variable parameter method, a variable parameter model of the system is constructed according to the state variable models corresponding to multiple equilibrium points of the system under different working conditions; A scheduling parameter θ is introduced to reflect the current working state of the system, and the state variable model at each equilibrium point is regarded as an affine function of the linear variable parameter model with respect to the scheduling parameter, that is, △y c =C c (θ)△x c +D c (θ)△u c +n Each element in the system matrix is ​​a variable value of the scheduling parameter, that is, A c (θ) = {a c,ij (θ)}, where a c,ij Represents the matrix A c The elements in ω and ν represent the process noise and observation noise of the system respectively, both of which obey the zero-mean Gaussian distribution, and the variances are Q and R respectively; the low-pressure turbine speed n of the aircraft engine is selected L As the dispatch parameter of the linear variable parameter model; In order to determine the benchmark values ​​of model input, output and state under different working conditions Wherein, k represents the model reference value of the kth operating point, and multiple m-order polynomials are used to approximate the relationship between the scheduling parameters and each reference value, as shown below; In the fitting process, the value of t is continuously increased until the coefficient of determination R of the fitting is 2 Stop when it is greater than 90%, t=min{t:R 2 (x c,st,n (n L,c ,t))>=0.9}; By calculating each state variable, output and output reference value through the above formula, the steady-state output of the system for the corresponding working state is finally obtained; In order to express the dynamic characteristics of the system, each element in the state space system matrix of multiple working points is fitted to obtain a power series polynomial about the scheduling parameters. The specific expression is as follows: Among them, ta, tb, tc, and td represent the polynomial orders of different system matrices, and the values ​​of each order are continuously increased during the fitting process until the determination coefficient R of the fitting is 2 When it is greater than 90%, it stops; finally, the adaptive dynamic model of the system under large working conditions is obtained. For the working point n L,1 The state value and output value of the system are obtained through the interaction of the corresponding reference value and the variable parameter state variable model; y c (n L,1 )=y c,st (n L,1 ,u c,st )+△y c (n L,1 ,△u c ) Among them, △u c It represents the increment of the input of the system at this moment relative to the previous sampling moment. The right side of the equal sign represents the sum of the baseline value of the observed output corresponding to the operating point and the input increment response. Similarly, the dynamic characteristics of the state variables of the system at different operating points are represented by a similar method.

3. The method for aircraft engine system identification and output interval estimation taking into account performance degradation and model mismatch according to claim 2, characterized in that: The step 2 is specifically as follows: in order to estimate the system performance degradation, a health parameter is introduced to represent the efficiency and flow deviation of the engine rotating parts; the health parameter h, which reflects the performance degradation of each engine component, is augmented to the state vector of the original system, that is, Since engine performance degradation is a slow process, the partial derivative of the performance degradation factor is approximated to 0; the augmented model is Right now, △y ext =C ext (θ)x ext +D ext (θ)△u c +n Among them A h (θ),B c (θ),C h (θ) is obtained by step 1, considering the change of estimated health parameters at the equilibrium point and ignoring △u c item; Considering the state variable model at a certain working point, the linear time-invariant model is first discretized to obtain x k =A k x k-1 +oh k y k =H k x k-1 +n k Among them, ω k ~N(0,Q k ) is the system noise that follows Gaussian distribution, and ν k ~N(0,Σ k ) is the observation noise, the noise variance Σ k is a diagonal matrix; assuming that the initial state value obeys Gaussian distribution, that is, x0~N(m0,P0), the subscript k represents the sampling time; The working process of the adaptive Kalman filter mainly includes the following steps: Calculate the forecast distribution parameters: Among them, α, β are parameters in the inverse gamma distribution, which are used to describe the variance of the observation noise, and ρ∈(0,1] is an adjustable constant; Parameter update step: First, let Then the following steps are iteratively solved; After N iterations, let Finally, the state distribution of the system at the working point is obtained, that is, the predicted value of the health parameter.

4. The method for aircraft engine system identification and output interval estimation taking into account performance degradation and model mismatch according to claim 1 or 2, characterized in that: Use the current working state of the system and the scheduling parameter θ as the input x, h of the long short-term memory neural network n Represents the network output; first, the input is saved to the unit state c through the input gate t middle; Among them, W i is the input gate weight matrix, b i is the bias term, Indicates the output status at the previous moment; the unit state of the current input of the system is, After that, the forget gate determines whether to keep the cell state c at the previous moment. t-1 ; Determine the cell state of the LSTM neural network at the current moment; in Represents the Hadamard product; according to the current unit state, the final output of the long short-term memory neural network is obtained through the output gate; Based on the quantile regression method, a fully connected layer is used to obtain the interval estimation of the system residual; the quantile level q∈Q={0.05,0.5,0.95} is taken to obtain the prediction results under different quantiles. At this time, the loss function of quantile regression is Where (y i ) + =max(y i ,0), the output when the quantile is 0.5 Represents the residual compensation term for the LPV system, that is, The prediction results when taking the other two quantiles represent interval predictions with 90% confidence.

Images