Method for fault diagnosis of aero-engine sensor and actuator based on LFT

11203446 · 2021-12-21

Assignee

Inventors

Cpc classification

International classification

Abstract

The present invention discloses a method for fault diagnosis of the sensors and actuators of an aero-engine based on LFT, and belongs to the field of fault diagnosis of aero-engines. The method comprises: establishing an aero-engine state space model using a combination of a small perturbation method and a linear fitting method; establishing an affine parameter-dependent linear-parameter-varying (LPV) model of the aero-engine based on the model; converting the LPV model of the aero-engine having perturbation signals and sensor and actuator fault signals into a linear fractional transformation (LFT) structure to obtain an synthesis framework of an LPV fault estimator; solving a set of linear matrix inequalities (LMIs) to obtain the solution conditions of the fault estimator; and designing the fault estimator in combination with the LFT structure to realize fault diagnosis of the sensors and actuators of an aero-engine.

Claims

1. A method for fault diagnosis of the sensor and actuator of an aero-engine based on LFT, comprising the following steps: step 1: establishing an aero-engine state space model using a combination of a small perturbation method and a linear fitting method; step 1.1: inputting fuel pressure p.sub.f under steady operating points into an aero-engine; inputting a fuel pressure step signal U.sub.pf1 with an amplitude of 0.01 p.sub.f into the aero-engine after relative conversion speed n.sub.h of a high pressure rotor of the aero-engine reaches a corresponding steady state; and respectively collecting relative conversion speed response Y.sub.nh1 of the high pressure rotor and relative conversion speed response Y.sub.nl1 of a low pressure rotor outputted by the aero-engine; step 1.2: repeating the process of the step 1 for N times; and respectively collecting the relative conversion speed response Y.sub.nhi of the high pressure rotor and the relative conversion speed response Y.sub.nli of the low pressure rotor under given fuel pressure p.sub.fi, wherein i=1,2,3, . . . ,N; step 1.3: constructing a discrete small perturbation state space model of the aero-engine under the steady operating points according to the linear fitting method by taking the fuel pressure step signal U.sub.pfi as an input variable and taking the relative conversion speed response Y.sub.nhi of the high pressure rotor and the relative conversion speed response Y.sub.nli of the low pressure rotor as state variables; step 1.4: converting the discrete small perturbation state space model of the aero-engine under the steady operating points into a continuous small perturbation state space model according to a sampling period T to obtain the state space model of the aero-engine; { x . p = A pi x p + B pi u y p = C pi x p + D pi u ( 1 ) wherein the state variable is x.sub.p=[Y.sub.nl.sup.T Y.sub.nh.sup.T].sup.T∈R.sup.n; {dot over (x)}.sub.p represents a first derivative of x.sub.p; an input variable is u=U.sub.pf∈R.sup.t; an output variable is y.sub.p=Y.sub.nh∈R.sup.m; A.sub.pi, B.sub.pi, C.sub.pi and D.sub.pi are system state space matrices; C.sub.pi=C.sub.p=[0 1]; D.sub.pi=D.sub.p=0; R.sup.n, R.sup.t and R.sup.m respectively represent sets of real numbers with dimensions of n, t and m; T represents transposing for the matrices; step 2: establishing an affine parameter-dependent linear-parameter-varying LPV model of the aero-engine; step 2.1: setting the relative conversion speed n.sub.hi of the high pressure rotor of the aero-engine as a scheduling parameter θ(i), i=1,2,3, . . . ,N; step 2.2: expressing a system matrix A.sub.p(θ) and a control matrix B.sub.p(θ) of the continuous small perturbation state space model of the aero-engine as affine parameter-dependent forms, as follows:
A.sub.p(θ)=A.sub.0+θA.sub.1, B.sub.p(θ)=B.sub.0+θB.sub.1   (2) wherein A.sub.0, A.sub.1, B.sub.0 and B.sub.1 respectively represent coefficient matrices to be solved; rewriting the formula (2) into A p ( θ ) = [ I θ I ] [ A 0 A 1 ] , B p ( θ ) = [ I θ I ] [ B 0 B 1 ] ( 3 ) wherein I is a unit matrix; then [ A 0 A 1 ] = [ I θ I ] + A p ( θ ) , [ B 0 B 1 ] = [ I θ I ] + B p ( θ ) ( 4 ) wherein [I θI].sup.+ is Moore-Penrose pseudo- inverse of [I θI], i.e., the system matrix A.sub.p(θ) and the control matrix B.sub.p(θ) of the solved affine parameter-dependent LPV model of the aero-engine; step 2.3: establishing the affine parameter-dependent LPV model of the aero-engine
{dot over (x)}.sub.p=A.sub.p(θ)x.sub.p+B.sub.p(θ)u
y.sub.p=C.sub.px.sub.p+D.sub.pu   (5); step 3: converting the affine parameter-dependent LPV model of the aero-engine with perturbation and sensor and actuator fault into a linear fractional transformation (LFT) structure, and establishing an H.sub.∞ synthesis framework of an LPV fault estimator of the aero-engine; step 3.1: expressing the affine parameter-dependent LPV model P(s,θ) of the aero-engine having perturbation and sensor and actuator fault into
{dot over (x)}.sub.p=A.sub.p(θ)x.sub.p+B.sub.p(θ)u+E.sub.pd+F.sub.pf
y.sub.p=C.sub.px.sub.p+D.sub.pu+G.sub.pd+H.sub.pf   (6) wherein d∈R.sup.q is a perturbation signal; f∈R.sup.l is a fault signal comprising sensor fault and actuator fault; R.sup.q and R.sup.l respectively represent sets of real numbers with dimensions of q and l; E.sub.p, F.sub.p, G.sub.p and H.sub.p are system state space matrices; an upper LFT structure of P(s,θ) is expressed into { [ x . p z θ y p ] = [ A p B p θ B pw C p θ D p θθ D p θ w C pw D pw θ D pww ] [ x p w θ w ] , w = [ u d f ] w θ = Δ ( θ ) z θ ( 7 ) wherein an external input variable is w=[u.sup.T d.sup.T f.sup.T].sup.T∈R.sup.p1; w.sub.θ∈R.sup.r is an output variable of a time varying part Δ(θ)=θI; z.sub.θ∈R.sup.r is an input variable of the time varying part Δ(θ)=θI; A.sub.p, B.sub.pθ, B.sub.pw, C.sub.pθ, C.sub.pw, D.sub.pθθ, D.sub.pθw, D.sub.pwθ and D.sub.pww are system state space matrices; R.sup.p1 and R.sup.r respectively represent sets of real numbers with dimensions of p1 and r; p1=t+q+l, i.e., the dimension p1 of the external input variable w is equal to the sum of the dimension t of the input variable u of the aero-engine, the dimension q of the perturbation signal d and the dimension l of the fault signal f; step 3.2: setting the form of the fault estimator K(s,θ) as follows { x . K = A K ( θ ) x K + B K ( θ ) u K f ^ = C K ( θ ) x K + D K ( θ ) u K ( 8 ) wherein x.sub.K∈R.sup.k is a state variable of the fault estimator K(s,θ); {dot over (x)}.sub.K represents a first derivative of x.sub.K; R .sup.k represents a set of real numbers with a dimension of k; u.sub.K=[u.sup.T y.sub.p.sup.T].sup.T∈R.sup.p2 is an input variable of K(s,θ); p2=t+m, i.e., the dimension p2 of the input variable u.sub.K of K(s,θ) is equal to the sum of the dimension t of the input variable u of the aero-engine and the dimension m of the output variable y.sub.p of the aero-engine; {circumflex over (f)}∈R.sup.l is an output variable of K(s,θ), i.e., an estimated value of the fault signal f; A.sub.K(θ), B.sub.K(θ), C.sub.K(θ) and D.sub.K(θ) are system state space matrices; K(s,θ) is express into a lower LFT structure as follows: { [ x . K f ^ z K ] = [ A K B K 1 B K θ C K 1 D K 11 D K 1 θ C K θ D K θ1 D K θθ ] [ x K u K w K ] , u K = [ u y p ] w K = Δ K ( θ ) z K ( 9 ) wherein w.sub.K∈R.sup.r is an output variable of the time varying part Δ.sub.K(θ)=θI; z.sub.K∈R.sup.r is an input variable of the time varying part Δ.sub.K(θ)=A.sub.K, B.sub.K1, B.sub.Kθ, C.sub.K1, C.sub.Kθ, D.sub.K11, D.sub.K1θ, D.sub.Kθ1 and D.sub.Kθθ are system state space matrices; step 3.3: according to the time varying part Δ(θ) in the LPV model P(s,θ) of the aero-engine and the time varying part Δ.sub.K(θ) in the fault estimator K(s,θ), expressing the H.sub.∞ synthesis framework of the LPV fault estimator as: [ x . p x K z K z θ e f ] = [ A _ B _ θ B _ w C _ θ D _ θθ D _ θ w C _ w D _ w θ D _ ww ] [ x p x K w K w θ w ] ( 10 ) wherein e.sub.f={circumflex over (f)}−f is a fault estimation error; system matrix Ā=A.sub.0+T.sub.1ΩT.sub.2; system matrix B.sub.θ=B.sub.01+T.sub.1ΩT.sub.3; system matrix B.sub.w=B.sub.02+T.sub.1ΩT.sub.4; system matrix C.sub.θ=C.sub.01+T.sub.5ΩT.sub.2; system matrix D.sub.θθ=D.sub.01+T.sub.5ΩT.sub.3; system matrix D.sub.θw=D.sub.02+T.sub.5ΩT.sub.4; system matrix C.sub.w=C.sub.02+T.sub.6ΩT.sub.2; system matrix D.sub.wθ=D.sub.03+T.sub.6ΩT.sub.3; system matrix D.sub.ww=D.sub.04+T.sub.6ΩT.sub.4; fault estimator matrix Ω = [ A K B K 1 B K θ C K 1 D K 11 D K 1 θ C K θ D K θ1 D K θθ ] ; matrix T 1 = [ 0 B 2 0 n × r I k 0 0 ] ; matrix T 2 = [ 0 I k C 2 0 0 r × n 0 ] ; matrix T 3 = [ 0 k × r 0 0 D 2 θ I r 0 ] ; matrix T 4 = [ 0 k × p 1 D 21 0 r × p 1 ] ; matrix T 5 = [ 0 r × k 0 I r 0 D θ 2 0 ] ; matrix T.sub.6=[0.sub.p1×k D.sub.12 0.sub.p1×r]; matrix A 0 = [ A 0 0 0 k ] ; matrix B 0 1 = [ 0 B θ 0 k × r 0 ] ; matrix B 0 2 = [ B 1 0 k × p 1 ] ; matrix C 0 1 = [ 0 0 r × k C θ 0 ] ; matrix D 01 = [ 0 r 0 0 D θθ ] ; matrix D 02 = [ 0 r × p 1 D θ 1 ] ; matrix C.sub.02=[C.sub.1 0.sub.p1×k]; matrix D.sub.03=[0.sub.p1×r D.sub.1θ]; matrix D.sub.04=D.sub.11; matrix A=A.sub.p; matrix B.sub.θ=B.sub.pθ; matrix B.sub.1=B.sub.pw; matrix B.sub.2=0.sub.n×l; matrix C.sub.θ=C.sub.pθ; matrix D.sub.θθ=D.sub.pθθ; matrix D.sub.θ1=D.sub.pθw; matrix D.sub.θ2=0.sub.r×l; matrix C.sub.1=0.sub.p1×n; matrix D.sub.1θ=0.sub.p1×r; matrix D 11 = [ 0 l × ( p 1 - l ) - I l 0 ( p 1 - l ) × ( p 1 - l ) 0 ( p 1 - l ) × l ] ; matrix D 12 = [ I l 0 ( p 1 - l ) × l ] ; matrix C 2 = [ 0 t × n C pw ] ; matrix D 2 θ = [ 0 t × r C pw θ ] ; matrix D 21 = [ I t 0 t × q 0 t × l D pww ] ; D.sub.22=0.sub.p2×l; n represents the dimension of the state variable x.sub.p of the aero-engine; r represents the dimension of the output variable w.sub.θ of the time varying part Δ(θ) and the output variable w.sub.K of the time varying part Δ.sub.K(θ); k represents the dimension of the state variable x.sub.K of the fault estimator K(s,θ); step 4: solving a set of linear matrix inequalities (LMIs) to obtain the solution conditions of the fault estimator; step 4.1: obtaining the solution conditions of the fault estimator K(s,θ), i.e., [ I 0 0 0 I 0 0 0 I A _ B _ θ B _ w C _ θ D _ θθ D _ θ w C _ w D _ w θ D _ ww ] T [ 0 0 0 X 0 0 0 Q 0 0 S 0 0 0 - γ I 0 0 0 X 0 0 0 0 0 0 S T 0 0 R 0 0 0 0 0 0 1 γ I ] [ I 0 0 0 I 0 0 0 I A _ B _ θ B _ w C _ θ D _ θθ D _ θ w C _ w D _ w θ D _ ww ] < 0 ( 11 ) [ Δ ( θ ) 0 0 Δ K ( θ ) I 0 0 I ] T P [ Δ ( θ ) 0 0 Δ K ( θ ) I 0 0 I ] > 0 ( 12 ) wherein X is a symmetric positive-definite matrix; a full block scaling matrix P = [ Q S S T R ] is a symmetric matrix; γ>0 is a performance index; Q, S and R respectively represent sub-block matrices of P; step 4.2: partitioning the symmetric positive-definite matrix X and an inverse matrix X.sup.−1 thereof; X = [ L M M T E ] , X - 1 = [ J N N T F ] ( 13 ) wherein L, M and E respectively represent block matrices of X; J, N and F respectively represent sub-block matrices of X.sup.−1; partitioning the full block scaling matrix P and the inverse matrix {tilde over (P)} thereof P = [ Q S S T R ] = [ Q 1 Q 2 Q 2 T Q 3 S 1 S 2 S 3 S 4 S 1 T S 3 T S 2 T S 4 T R 1 R 2 R 2 T R 3 ] , P ~ = [ Q ~ S ~ S ~ T R ~ ] = [ Q ~ 1 Q ~ 2 Q ~ 2 T Q ~ 3 S ~ 1 S ~ 2 S ~ 3 S ~ 4 S ~ 1 T S ~ 3 T S ~ 2 T S ~ 4 T R ~ 1 R ~ 2 R ~ 2 T R ~ 3 ] ( 14 ) wherein Q.sub.1, Q.sub.2 and Q.sub.3 respectively represent sub-block matrices of Q; S.sub.1, S.sub.2, S.sub.3 and S.sub.4 respectively represent sub-block matrices of S; R.sub.1, R.sub.2 and R.sub.3 respectively represent sub-block matrices of R; {tilde over (Q)}, {tilde over (S)} and {tilde over (R)} respectively represent sub-block matrices of {tilde over (P)}; {tilde over (Q)}.sub.1, {tilde over (Q)}.sub.2 and {tilde over (Q)}.sub.3 respectively represent sub-block matrices of {tilde over (Q)}; {tilde over (S)}.sub.1, {tilde over (S)}.sub.2, {tilde over (S)}.sub.3 and {tilde over (S)}.sub.4 respectively represent sub-block matrices of {tilde over (S)}; {tilde over (R)}.sub.1, {tilde over (R)}.sub.2 and {tilde over (R)}.sub.3 respectively represent sub-block matrices of {tilde over (R)}; simplifying the solution conditions of the fault estimator K(s,θ), i.e., N L T [ I 0 0 0 I 0 0 0 I A B θ B 1 C θ D θθ D θ 1 C 1 D 1 θ D 11 ] T [ 0 0 0 L 0 0 0 Q 3 0 0 S 4 0 0 0 - γ I 0 0 0 L 0 0 0 0 0 0 S 4 T 0 0 R 3 0 0 0 0 0 0 1 γ I ] [ I 0 0 0 I 0 0 0 I A B θ B 1 C θ D θθ D θ 1 C 1 D 1 θ D 11 ] N L < 0 ( 15 ) N J T [ - A T - C θ T - C 1 T - B θ T - D θθ T - D 1 θ T - B 1 T - D θ 1 T - D 11 T I 0 0 0 I 0 0 0 I ] T [ 0 0 0 J 0 0 0 Q ~ 3 0 0 S ~ 4 0 0 0 - 1 γ I 0 0 0 J 0 0 0 0 0 0 S ~ 4 T 0 0 R ~ 3 0 0 0 0 0 0 γ I ] [ - A T - C θ T - C 1 T - B θ T - D θθ T - D 1 θ T - B 1 T - D θ 1 T - D 11 T I 0 0 0 I 0 0 0 I ] N J > 0 ( 16 ) [ L I I J ] ( 17 ) R > 0 , Q = - R , S + S T = 0 ( 18 ) wherein N.sub.L and N.sub.J respectively represent the bases of the nuclear spaces of [C.sub.2 D.sub.2θ D.sub.21] and [B.sub.2.sup.T D.sub.θ2.sup.T D.sub.12.sup.T]; step 4.3: solving the LMIs (15)-(18) to obtain matrix solutions L, J, Q.sub.3, {tilde over (R)}.sub.3, S.sub.4 and {tilde over (S)}.sub.4; step 5: designing the fault estimator in combination with the LFT structure to realize fault diagnosis of the sensor and actuator of the aero-engine; step 5.1: solving the symmetric positive-definite matrix X, the full block scaling matrix P and the inverse matrix {tilde over (P)} thereof from the formulas (13) and (14) according to the solved matrix solutions L, J, Q.sub.3, {tilde over (R)}.sub.3, S.sub.4 and {tilde over (S)}.sub.4; step 5.2: according to Schur complement Lemma, expressing the LMI (11) as [ A _ T X + X A _ X B _ θ + C _ θ T S T X B _ w C _ θ T C _ w T B _ θ T X + S C _ θ Q + D _ θθ T S T + S D _ θθ S D _ θ w D _ θθ T D _ θ w T B _ w T X D _ θ w T S T - γ I D _ θ w T D _ w w T C _ θ D _ θθ D _ θ w - R 0 C _ w D _ w θ D _ ww 0 - γ I ] < 0 ( 19 ) solving the LMI (19) to obtain a fault estimator matrix Ω; step 5.3: obtaining a state space matrix of the fault estimator K(s,θ) [ A K ( θ ) B K ( θ ) C K ( θ ) D K ( θ ) ] = [ A K B K 1 C K 1 D K 11 ] + [ B K θ D K 1 θ ] Δ K ( θ ) ( I - D K θ θ Δ K ( θ ) ) - 1 [ C K θ D K θ I ] . ( 20 )

Description

DESCRIPTION OF DRAWINGS

(1) FIG. 1 is a contrast curve of relative conversion speed response Y.sub.nh of a high pressure rotor of a state space model of an aero-engine and test data under H=0,Ma=0,n.sub.2=90% operating state.

(2) FIG. 2 is a contrast curve of relative conversion speed response Y.sub.nh of a high pressure rotor of an LPV model of an aero-engine and test data under H=0,Ma=0,n.sub.2=90% operating state.

(3) FIG. 3 is a structural diagram of an upper LFT of an LPV model P(s,θ) of an aero-engine.

(4) FIG. 4 is a system structural diagram under an LFT framework.

(5) FIG. 5 is an H.sub.∞ synthesis framework of an LPV fault estimator.

(6) FIG. 6(a) and FIG. 6(b) are simulation results of sudden fault estimation.

(7) FIG. 7(a) and FIG. 7(b) are simulation results of slow fault estimation.

(8) FIG. 8(a) and FIG. 8(b) are simulation results of intermittent fault estimation.

(9) FIG. 9 is a flow chart of the present invention.

DETAILED DESCRIPTION

(10) The embodiments of the present invention will be further described in detail below in combination with the drawings and the technical solution.

(11) The flow chart of the present invention is shown in FIG. 9, and comprises the following specific steps:

(12) step 1.1: inputting fuel pressure p.sub.f under steady operating points into an aero-engine; inputting a fuel pressure step signal U.sub.pf1 with an amplitude of 0.01 p.sub.f into the aero-engine after relative conversion speed n.sub.h of a high pressure rotor of the aero-engine reaches a corresponding steady state; and respectively collecting relative conversion speed response Y.sub.nh1 of the high pressure rotor and relative conversion speed response Y.sub.nl1 of a low pressure rotor outputted by the aero-engine.

(13) Step 1.2: repeating the above process for 13 times, i.e., respectively collecting the relative conversion speed response Y.sub.nhi of the high pressure rotor and the relative conversion speed response Y.sub.nli of the low pressure rotor under given fuel pressure p.sub.fi at balance points of 13 working conditions of (H=0,Ma=0,n.sub.h=88%, 89%, . . . , 100%), wherein i=1,2,3, . . . ,13.

(14) Step 1.3: by taking the fuel pressure step signal U.sub.pfi as an input variable and taking the relative conversion speed response Y.sub.nhi of the high pressure rotor and the relative conversion speed response Y.sub.nli of the low pressure rotor as state variables, expressing a discrete small perturbation state space model of the aero-engine as

(15) { x p k + 1 = E i x p k + F i u k y p k = G i x p k + H i u k ( 21 )

(16) wherein the state variable x.sub.p=[Y.sub.nl Y.sub.nh].sup.T∈R.sup.n; the input variable u=U.sub.pf∈R.sup.t; the output variable y.sub.p=Y.sub.nh∈R.sup.m, i=1,2,3, . . . ,13; the subscript k, k+1 are corresponding sampling moments; E.sub.i, F.sub.i, G.sub.i and H.sub.i are system state space matrices with appropriate dimensions; R.sup.n, R.sup.t and R.sup.m respectively represent sets of real numbers with dimensions of n, t and m; T represents transposing for the matrices. According to the basic idea of the fitting method, a linear least square problem is established for formula (21), and the system matrices E.sub.i, F.sub.i, G.sub.i, H.sub.i are solved by using the lsqnonlin function in MATLAB.

(17) Step 1.4: converting the discrete small perturbation state space model of the aero-engine under the steady operating points into a continuous small perturbation state space model according to a sampling period T=25 ms to obtain the state space model of the aero-engine;

(18) { x . p = A pi x p + B pi u y p = C pi x p + D pi u ( 22 )

(19) wherein A.sub.pi, B.sub.pi, C.sub.pi and D.sub.pi are system state space matrices with appropriate dimensions; C.sub.pi=C.sub.p=[0 1], and D.sub.pi=D.sub.p=0; and a relative conversion speed response Y.sub.nh curve of the high pressure rotor of the state space model at the operating pint H=0,Ma=0,n.sub.2=90% is provided, as shown in FIG. 1, and has an average relative error of 0.26% relative to the test data.

(20) Step 2.1: setting the relative conversion speed n.sub.hi of the high pressure rotor of the aero-engine as a scheduling parameter θ(i),i=1,2,3, . . . ,13.

(21) Step 2.2: expressing a system matrix A.sub.p(θ) and a control matrix B.sub.p(θ) of the continuous small perturbation state space model of the aero-engine as affine parameter-dependent forms, as follows:
A.sub.p(θ)=A.sub.0+θA.sub.1
B.sub.p(θ)=B.sub.0+θB.sub.1   (23)

(22) wherein A.sub.0, A.sub.l, B.sub.0 and B.sub.1 respectively represent coefficient matrices to be solved.

(23) Rewriting the formula (23) into

(24) [ I θ ( 1 ) I I θ ( 2 ) I .Math. .Math. I θ ( 13 ) I ] [ A 0 A 1 ] = [ A p ( θ ( 1 ) ) A p ( θ ( 2 ) ) .Math. A p ( θ ( 13 ) ) ] , [ I θ ( 1 ) I I θ ( 2 ) I .Math. .Math. I θ ( 13 ) I ] [ B 0 B 1 ] = [ B p ( θ ( 1 ) ) B p ( θ ( 2 ) ) .Math. B p ( θ ( 13 ) ) ] ( 24 )

(25) wherein I is a unit matrix.

(26) Then

(27) [ A 0 A 1 ] = [ I θ ( 1 ) I I θ ( 2 ) I .Math. .Math. I θ ( 13 ) I ] + [ A p ( θ ( 1 ) ) A p ( θ ( 2 ) ) .Math. A p ( θ ( 13 ) ) ] , [ B 0 B 1 ] = [ I θ ( 1 ) I I θ ( 2 ) I .Math. .Math. I θ ( 13 ) I ] + [ B p ( θ ( 1 ) ) B p ( θ ( 2 ) ) .Math. B p ( θ ( 13 ) ) ] ( 25 )

(28) The pinv function in MATLAB is used to solve Moore-Penrose pseudo- inverse

(29) [ I θ ( 1 ) I I θ ( 2 ) I .Math. .Math. I θ ( 13 ) I ] + of [ I θ ( 1 ) I I θ ( 2 ) I .Math. .Math. I θ ( 13 ) I ] ,
and variable transformation is conducted on a variable parameter θ, so that θ∈[−1,1], to obtain:

(30) A 0 = [ - 2.5839 0.5968 1.1613 - 4.3763 ] , A 1 = [ 0.4287 - 2.3152 0.3138 - 0.7456 ] , B 0 = [ 0.0037 0.0021 ] , B 1 = [ - 0.0003 - 0.0002 ] ( 26 )

(31) Step 2.3: establishing the affine parameter-dependent LPV model of the aero-engine
{dot over (x)}.sub.p=A.sub.p(θ)x.sub.p+B.sub.p(θ)u
y.sub.p=C.sub.px.sub.p+D.sub.pu   (27)

(32) wherein a relative conversion speed response Y.sub.nh curve of the high pressure rotor of the LPV model of the aero-engine at the operating pint H=0,Ma=0,n.sub.2=90% is provided, as shown in FIG. 2, and has an average relative error of 2.51% relative to the test data.

(33) Step 3.1: expressing the affine parameter-dependent LPV model P(s,θ) of the aero-engine having perturbation and sensor and actuator fault into
x.sub.p=A.sub.p(θ)x.sub.p+B.sub.p(θ)u+E.sub.pd+F.sub.pf
y.sub.p=C.sub.px.sub.p+D.sub.pu+G.sub.pd+H.sub.pf   (28)

(34) wherein d∈R.sup.q is a perturbation signal and takes Gaussian white noise with standard deviation of 0.001; f∈R.sup.l is a fault signal comprising sensor fault and actuator fault, which respectively take sudden fault, slow fault and intermittent fault; R.sup.q and R.sup.l respectively represent sets of real numbers with dimensions of q and

(35) l ; E p = [ 0.01 0 ] , F p = [ 0 0.1 ] ,
G.sub.p=0.2, H.sub.p=1.

(36) The upper LFT structure of P (s,θ) can be expressed as the following formula, as shown in FIG. 3:

(37) y p = F u ( P , Δ ( θ ) ) [ u d f ] ( 29 )

(38) wherein F.sub.u represents the upper LFT structure; P′ represents a time-invariable part in P(s,θ); Δ(θ)=θI represents a time varying part in P(s,θ), i.e.,

(39) 0 { [ x . p z θ y p ] = [ A p B p θ B pw C p θ D p θθ D p θ w C pw D pw θ D pww ] [ x p w θ w ] , w = [ u d f ] w θ = Δ ( θ ) z θ ( 30 )

(40) wherein an external input variable is w=[u.sup.T d.sup.T f.sup.T].sup.T∈R.sup.p1; w.sub.θ∈R.sup.r is an output variable of a time varying part Δ(θ)=θI; z.sub.θ∈R.sup.r is an input variable of the time varying part Δ(θ)=θI; R.sup.p1 and R.sup.r respectively represent sets of real numbers with dimensions of p1 and r; p1=t+q+l, i.e., the dimension p1 of the external input variable w is equal to the sum of the dimension t of the input variable u of the aero-engine, the dimension q of the perturbation signal d and the dimension l of the fault signal f; the system state space matrices are

(41) A p = [ - 2.5839 0.5968 1.1613 - 4.3763 ] B p = [ 0.4287 0 - 2.3152 0 - 0.0003 0 0 0.3138 0 - 0.7456 0 - 0.0002 ] , B pw = [ 0.0037 0.01 0 0.0021 0 1 ] C p θ = [ 1 1 0 0 0 0 0 0 1 1 0 0 ] T , C pw = [ 0 1 ] D p θθ = 0 6 × 6 , D p θ w = [ 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 ] T , D pw θ = 0 1 × 6 , D pww = [ 0 0.2 1 ] ( 31 )

(42) Step 3.2: setting the form of the fault estimator K(s,θ) as follows

(43) { x . K = A K ( θ ) x K + B K ( θ ) u K f ^ = C K ( θ ) x K + D K ( θ ) u K ( 32 )

(44) wherein X.sub.K∈R.sup.k is a state variable of the fault estimator K (s,θ); u.sub.K=[u.sup.T y.sub.p.sup.T].sup.T∈R.sup.p2 is an input variable of K(s,θ); p2=t+m, i.e., the dimension p2 of the input variable u.sub.K of K (s,θ) is equal to the sum of the dimension t of the input variable u of the aero-engine and the dimension m of the output variable y.sub.p of the aero-engine; {circumflex over (f)}∈R.sup.l is an output variable of K (s,θ), i.e., the estimated value of the fault signal f; A.sub.K(θ), B.sub.K(θ), C.sub.K(θ) and D.sub.K(θ) are system state space matrices. K(s,θ) is expressed as a lower LFT structure as follows:

(45) f ^ = F l ( K , Δ K ( θ ) ) [ u y p ] ( 33 )

(46) wherein F.sub.l represents the lower LFT structure; K′ represents a time-invariable part in K(s,θ); Δ.sub.K(θ)=θI represents a time varying part in K(s,θ), i.e.,

(47) { [ x . K f ^ z K ] = [ A K B K 1 B K θ C K 1 D K 11 D K 1 θ C K θ D K θ1 D K θθ ] [ x K u K w K ] , u k = [ u y p ] w K = Δ K ( θ ) z K ( 34 )

(48) wherein w.sub.K∈R.sup.r is an output variable of the time varying part Δ.sub.K(θ)=θI; z.sub.K∈R.sup.r is an input variable of the time varying part Δ.sub.K(θ)=θI; A.sub.K, B.sub.K1, B.sub.Kθ, C.sub.K1, C.sub.Kθ, D.sub.K11, D.sub.K1θ, D.sub.Kθ1 and D.sub.Kθθ are system state space matrices with appropriate dimensions.

(49) Step 3.3: showing a system connection diagram under the LFT framework in FIG. 4, and giving a state space expression of the system P.sub.1 in FIG. 4 as

(50) [ x . z θ e f u k ] = [ A B θ B 1 B 2 C θ D θθ D θ 1 D θ 2 C 1 D 1 θ D 11 D 12 C 2 D 2 θ D 21 D 22 ] [ x w θ w f ^ ] ( 35 )

(51) wherein system matrix A=A.sub.p; system matrix B.sub.θ=B.sub.pθ; system matrix B.sub.1=B.sub.pw; system matrix B.sub.2=0.sub.n×l; system matrix C.sub.θ=C.sub.pθ; system matrix D.sub.θθ=D.sub.pθθ; system matrix D.sub.θ1=D.sub.pθw; system matrix D.sub.θ2=0.sub.r×l; system matrix C.sub.1=0.sub.p1×n; system matrix D.sub.1θ=0.sub.p1×r; system matrix

(52) D 11 = [ 0 l × ( p 1 - l ) - I l 0 ( p 1 - l ) × ( p 1 - l ) 0 ( p 1 - l ) × l ] ;
system matrix

(53) D 12 = [ I l 0 ( p 1 - l ) × l ] ;
system matrix

(54) C 2 = [ 0 t × n C pw ] ;
system matrix

(55) D 2 θ = [ 0 t × r D pw θ ] ;
system matrix

(56) 0 D 2 1 = [ I t 0 t × q 0 t × l D pww ] ;
D.sub.22=0.sub.p2×l; n represents the dimension of the state variable x.sub.p of the aero-engine; r represents the dimension of the output variable w.sub.θ of the time varying part Δ(θ) and the output variable w.sub.K of the time varying part Δ.sub.K(θ).

(57) According to the time varying part Δ(θ) in the LPV model P(s,θ) of the aero-engine and the time varying part Δ.sub.K(θ) in the fault estimator K (s,θ), expressing the H.sub.∞ synthesis framework of the LPV fault estimator as follows, as shown in FIG. 5

(58) [ x . x . K z K z θ e f ] = [ A _ B _ θ B _ w C _ θ D _ θθ D _ θ w C _ w D _ w θ D _ ww ] [ x x K w K w θ w ] ( 36 )

(59) wherein e.sub.f={circumflex over (f)}−f is a fault estimation error, i.e., an output variable of the H.sub.∞ synthesis framework of the LPV fault estimator; system matrix Ā=A.sub.0+T.sub.1ΩT.sub.2; system matrix B.sub.θ=B.sub.01+T.sub.1ΩT.sub.3; system matrix B.sub.w=B.sub.02+T.sub.1ΩT.sub.4; system matrix C.sub.θ=C.sub.01+T.sub.5ΩT.sub.2; system matrix D.sub.θθ=D.sub.01+T.sub.5ΩT.sub.3; system matrix D.sub.θw=D.sub.02+T.sub.5ΩT.sub.4; system matrix C.sub.w=C.sub.02+T.sub.6ΩT.sub.2; system matrix D.sub.wθ=D.sub.03+T.sub.6ΩT.sub.3; system matrix D.sub.ww=D.sub.04+T.sub.6ΩT.sub.4; fault estimator matrix

(60) Ω = [ A K B K 1 B K θ C K 1 D K 11 D K 1 θ C K θ D K θ 1 D K θθ ] ;
matrix

(61) T 1 = [ 0 B 2 0 n × r I k 0 0 ] ;
matrix

(62) T 2 = [ 0 I k C 2 0 0 r × n 0 ] ;
matrix

(63) T 3 = [ 0 k × r 0 0 D 2 θ I r 0 ] ;
matrix

(64) T 4 = [ 0 k × p 1 D 21 0 r × p 1 ] ;
matrix

(65) T 5 = [ 0 r × k 0 I r 0 D θ 2 0 ] ;
matrix T.sub.6=[0.sub.p1×k D.sub.12 0.sub.p1×r]; matrix

(66) A 0 = [ A 0 0 0 k ] ;
matrix

(67) B 01 = [ 0 B θ 0 k × r 0 ] ;
matrix

(68) 0 B 02 = [ B 1 0 k × p 1 ] ;
matrix

(69) C 01 = [ 0 0 r × k C θ 0 ] ;
matrix

(70) D 01 = [ 0 r 0 0 D θθ ] ;
matrix

(71) D 02 = [ 0 r × p 1 D θ 1 ] ;
matrix C.sub.02=[C.sub.1 0.sub.p1×k]; matrix D.sub.03=[0.sub.p1×r D.sub.1θ]; and matrix D.sub.04=D.sub.11.

(72) Step 4.1: if there is a symmetric positive-definite matrix X, the symmetrical matrix

(73) P = [ Q S S T R ]
making the formula (37) and the formula (38) valid;

(74) [ I 0 0 0 I 0 0 0 I A _ B _ θ B _ w C _ θ D _ θθ D _ θ w C _ w D _ w θ D _ ww ] T [ 0 0 0 X 0 0 0 Q 0 0 S 0 0 0 - γ I 0 0 0 X 0 0 0 0 0 0 S T 0 0 R 0 0 0 0 0 0 1 γ I ] [ I 0 0 0 I 0 0 0 I A _ B _ θ B _ w C _ θ D _ θθ D _ θ w C _ w D _ w θ D _ ww ] < 0 ( 37 ) [ Δ ( θ ) 0 0 Δ K ( θ ) I 0 0 I ] T P [ Δ ( θ ) 0 0 Δ K ( θ ) I 0 0 I ] > 0 ( 38 )

(75) A closed-loop system (36) is asymptotically stable, and the L.sub.2 induced norm of the closed-loop transfer function from the external input w to the fault estimation error e.sub.f is smaller than a performance index γ(γ>0). Namely, the solution conditions of the fault estimator K(s,θ) are formula (37) and formula (38),wherein Q, S and R respectively represent sub-block matrices of P.

(76) Step 4.2: partitioning the symmetric positive-definite matrix X and an inverse matrix X.sup.−1 thereof;

(77) X = [ L M M T E ] , X - 1 = [ J N N T F ] ( 39 )

(78) wherein L, M and E respectively represent block matrices of X; J, N and F respectively represent sub-block matrices of X.sup.−1.

(79) Because x is the symmetric positive-definite matrix, then

(80) [ L I I J ] > 0 ( 40 )

(81) Partitioning the full block scaling matrix P and the inverse matrix {tilde over (P)} thereof

(82) P = [ Q S S T R ] = [ Q 1 Q 2 S 1 S 2 Q 2 T Q 3 S 3 S 4 S 1 T S 3 T R 1 R 2 S 2 T S 4 T R 2 T R 3 ] , P ~ = [ Q ~ S ~ S ~ T R ~ ] = [ Q ~ 1 Q ~ 2 S ~ 1 S ~ 2 Q ~ 2 T Q ~ 3 S ~ 3 S ~ 4 S ~ 1 T S ~ 3 T R ~ 1 R ~ 2 S ~ 2 T S ~ 4 T R ~ 2 T R ~ 3 ] ( 41 )

(83) wherein Q.sub.1, Q.sub.2 and Q.sub.3 respectively represent sub-block matrices of Q; S.sub.1, S.sub.2, S.sub.3 and S.sub.4 respectively represent sub-block matrices of S; R.sub.1, R.sub.2 and R.sub.3 respectively represent sub-block matrices of R; {tilde over (Q)}, {tilde over (S)} and {tilde over (R)} respectively represent sub-block matrices of {tilde over (P)}; {tilde over (Q)}.sub.1, {tilde over (Q)}.sub.2 and {tilde over (Q)}.sub.3 respectively represent sub-block matrices of {tilde over (Q)}; {tilde over (S)}.sub.1, {tilde over (S)}.sub.2, {tilde over (S)}.sub.3 and {tilde over (S)}.sub.4 respectively represent sub-block matrices of {tilde over (S)}; {tilde over (R)}.sub.1, {tilde over (R)}.sub.2 and {tilde over (R)}.sub.3 respectively represent sub-block matrices of {tilde over (R)}.

(84) Arranging the LMI (37) as

(85) [ I U Γ V + Z ] T [ 0 0 0 0 0 L M 0 0 0 0 0 0 0 0 M T E 0 0 0 0 0 Q 1 Q 2 0 0 0 S 1 S 2 0 0 0 Q 2 T Q 3 0 0 0 S 3 S 4 0 0 0 0 0 - γ I 0 0 0 0 0 L M 0 0 0 0 0 0 0 0 M T E 0 0 0 0 0 0 0 0 0 0 S 1 T S 3 T 0 0 0 R 1 R 2 0 0 0 S 2 T S 4 T 0 0 0 R 2 T R 3 0 0 0 0 0 0 0 0 0 0 1 γ I ] [ I U Γ V + Z ] < 0 ( 42 )

(86) wherein matrix

(87) 0 Z = [ A 0 B 01 B 02 C 01 D 01 D 02 C 02 D 03 D 04 ] ;
matrix

(88) U = [ T 1 T 5 T 6 ] ;
matrix V=[T.sub.2 T.sub.3 T .sub.4]; and matrix Γ=Ω.

(89) To satisfy the formula (38), it is required to verify that the formula (38) is valid on all possible tracks of the variable parameter θ, which is impossible. Therefore, the structure of the full block scaling matrix

(90) P = [ Q S S T R ]
is limited to make it valid. For each variable parameter θ, when R≥0, the following formula is valid

(91) [ θ I I ] T [ Q S S T R ] [ θ I I ] = θ 2 Q + R + θ ( S T + S ) θ 2 ( Q + R ) + θ ( S T + S ) 0 ( 43 )

(92) Therefore, Q=−R, and S+S.sup.T=0. Namely, the formula (38) is arranged as
R>0, Q=−R, S+S.sup.T=0   (44)

(93) To sum up, the solution conditions of the fault estimator K(s,θ) are converted into formula (40), formula (42) and formula (44).

(94) Step 4.3: arranging the LMI (42) as

(95) V T [ I Z ] T W [ I Z ] V < 0 ( 45 ) U T [ - Z T I ] T W - 1 [ - Z T I ] U > 0 ( 46 )

(96) wherein U.sub.⊥ and V.sub.⊥ are respectively the bases of the nuclear spaces of U.sup.T and V.

(97) W = [ 0 0 0 0 0 L M 0 0 0 0 0 0 0 0 M T E 0 0 0 0 0 Q 1 Q 2 0 0 0 S 1 S 2 0 0 0 Q 2 T Q 3 0 0 0 S 3 S 4 0 0 0 0 0 - γ I 0 0 0 0 0 L M 0 0 0 0 0 0 0 0 M T E 0 0 0 0 0 0 0 0 0 0 S 1 T S 3 T 0 0 0 R 1 R 2 0 0 0 S 2 T S 4 T 0 0 0 R 2 T R 3 0 0 0 0 0 0 0 0 0 0 1 γ I ]

(98) Simplifying the LMIs (45) and (46) through simple calculation as

(99) N L T [ I 0 0 0 I 0 0 0 I A B θ B 1 C θ D θθ D θ 1 C 1 D 1 θ D 11 ] T [ 0 0 0 L 0 0 0 Q 3 0 0 S 4 0 0 0 - γ I 0 0 0 L 0 0 0 0 0 0 S 4 T 0 0 R 3 0 0 0 0 0 0 1 γ I ] [ I 0 0 0 I 0 0 0 I A B θ B 1 C θ D θθ D θ 1 C 1 D 1 θ D 11 ] N L < 0 ( 47 ) N J T [ - A T - C θ T - C 1 T - B θ T - D θθ T - D 1 θ T - B 1 T - D θ1 T - D 11 T I 0 0 0 I 0 0 0 I ] T [ 0 0 0 J 0 0 0 Q ~ 3 0 0 S ~ 4 0 0 0 - 1 γ I 0 0 0 J 0 0 0 0 0 0 S ~ 4 T 0 0 R ~ 3 0 0 0 0 0 0 γ I ] [ - A T - C θ T - C 1 T - B θ T - D θθ T - D 1 θ T - B 1 T - D θ1 T - D 11 T I 0 0 0 I 0 0 0 I ] N J > 0 ( 48 )

(100) wherein N.sub.L and N.sub.J respectively represent the bases of the nuclear spaces of [C.sub.2 D.sub.2θ D.sub.21] and [B.sub.2.sup.T D.sub.θ2.sup.T D.sub.12.sup.T].

(101) Step 4.4: solving the LMIs (40), (44), (47) and (48) by using an LMI toolkit in MATLAB to obtain an optimal γ value of 0.21 and corresponding matrix solutions L, J, Q.sub.3, {tilde over (R)}.sub.3, S.sub.4 and {tilde over (S)}.sub.4.

(102) Step 5.1: solving the symmetric positive-definite matrix X, the full block scaling matrix P and the inverse matrix {tilde over (P)} thereof from the formulas (39) and (41) according to the solved matrix solutions L, J, Q.sub.3, {tilde over (R)}.sub.3, S.sub.4 and {tilde over (S)}.sub.4.

(103) Step 5.2: according to Schur complement Lemma, expressing the LMI (37) as

(104) [ A _ T X + X A _ X B _ θ + C _ θ T S T X B _ w C _ θ T C _ w T B _ θ T X + S C _ θ Q + D _ θθ T S T + S D _ θθ S D _ θ w D _ θθ T D _ w θ T B _ w T X D _ θ w T S T - γ I D _ θ w T D _ ww T C _ θ D _ θθ D _ θ w - R ~ 0 C _ w D _ w θ D _ ww 0 - γ I ] < 0 ( 49 )

(105) Substituting the value in the closed-loop system (36) to obtain
Ψ+P.sup.TΩ.sup.TQ.sub.x+Q.sub.X.sup.TΩP<0   (50)

(106) wherein

(107) Ψ = [ A 0 T X + XA 0 XB 01 + C 01 T S T XB 02 C 01 T C 02 T B 01 T X + SC 01 Q + D 01 T S T + SD 01 SD 02 D 01 T D 03 T B 02 T X D 02 T S T - γ I D 02 T D 04 T C 01 D 01 D 02 - R ~ 0 C 02 D 03 D 04 0 - γ I ] ,
P=[T.sub.2 T.sub.3 T.sub.4 0 0] and Q.sub.x=[T.sub.1.sup.TX T.sub.5.sup.TS.sup.T 0 T.sub.5.sup.T T.sub.6.sup.T].

(108) Solving the LMI (50) to obtain a fault estimator matrix Ω.

(109) Step 5.3: obtaining a state space matrix of the fault estimator K(s,θ)

(110) [ A K ( θ ) B K ( θ ) C K ( θ ) D K ( θ ) ] = [ A K B K 1 C K 1 D K 1 ] + [ B K θ D K 1 θ ] Δ K ( θ ) ( I - D K θθ Δ K ( θ ) ) - 1 [ C K θ D K θ 1 ] ( 51 )
Simulation results at the operating points H=0 km,Ma=0,n.sub.2=90% are shown in FIG. 6(a), FIG. 6(b), FIG. 7(a), FIG. 7(b), FIG. 8(a) and FIG. 8(b), and are compared with the standard H.sub.∞ method. The simulation results show that the fixed parameter fault estimator designed by the standard H.sub.∞ method cannot well cope with the change of the variable parameters. The LPV fault estimator designed by the present invention can rapidly detect the fault in the system and accurately reconstruct the fault signal, and has obvious performance advantages.