Polytope LPV estimation for non-linear flight control

Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011

Polytope LPV estimation for non-linear flight control Lejun Chen* and Ron J Patton* 

*Department of Engineering, University of Hull, Hull HU6 7RX, (e-mail: l.chen, [email protected]) Abstract: The current study is motivated by the need to implement the linear based model-based fault detection and isolation (FDI) methodology onto the nonlinear aircraft system directly. The nonlinear model is expressed in the linear parameter varying (LPV) manner and the corresponding LPV FDI estimator can be developed through the combination of the polytopic FDI estimators developed on each system vertex. The proposed design strategy is applied to the nonlinear longitudinal motion of a UAV aircraft (Machan) with different faults acting on the elevator actuator and wind turbulence affecting the vertical force. Keywords: Nonlinear Aircraft Application, Linear Parameter Varying, Fault Detection and Isolation, Linear Matrix Inequalities 

problems of fault detection (including delectability) and fault isolation are solved in a more efficient manner.

1. INTRODUCTION This paper focuses on the Linear Parameter Varying (LPV) polytope modeling methodology that has been widely adopted in control system design in recent years, especially related to vehicle and aerospace control (Wu 2001; Ganguli, Marcos and Balas 2002). Based upon the use of the LPV approach, the time-varying terms of a nonlinear system can be parameterised. The nonlinear system is modelled as a system with the linear structure, through the online measurement of the varying parameters. For control applications the LPV approach facilitates the direct application of classical control structures directly on the time-varying and non-linear system with robust results. FDD schemes based on LPV system have also been developed (Bokor and Balas 2004; Bokor and Kulcsar 2004; Henry and Zolghadri 2005; Casavola, Famularo, Famularo and Sorbara 2007; Casavola, Famularo, Franze and Patton 2008; Weng, Patton and Cui 2008). The LPV based FDD approaches can be divided into two categories (i) a Linear Fractional Representation formalism and (ii) the polytopic formalism. In this paper, the latter approach is adopted to design a robust LPV fault estimator for a nonlinear system. In contrast to the residual-based approach to fault detection and isolation (FDI), wherein the sensitivity of the fault is necessary, the fault estimation requires the achievement of the fault magnitude, error steady state error. Therefore, the requirement of maximizing the fault sensitivity in a residual can be transformed into one of minimizing the fault estimation error. In other words, the H  / H  fault residual generator can be viewed as an H  fault estimator. It is well known (see Chen and Patton 1999) that the ideal residual generator is a faithful estimator of the fault and this ideal is difficult to achieve in practice unless a deadbeat system design is made. By using a direct approach to robust fault estimation design (rather than residual design) the robustness 978-3-902661-93-7/11/$20.00 © 2011 IFAC

The literature of the development of the H  optimization based methods is well summarized in Marcos, Ganguli and Balas 2005. Recent work by (Patton, Putra and Klinkhieo 2009) used sliding mode estimation as a robust method of fault estimation, as a part of a fault-tolerant control compensation system. In Patton, Chen and Klinkhieo 2010, the friction effects acting in a two-manipulator robot system are viewed as actuator faults with time-varying characteristics to be estimated and compensated within a LPV Fault detection and diagnosis scheme. This approach is in general applicable to a wide class of non-linear systems with no unique equilibria (for linearization) and a wide range of fault conditions. The importance of this paper lies in the application of the robust H  LPV fault estimator approach to a very nonlinear aircraft system, the Machan UAV, including both longitudinal and lateral dynamics, wind turbulence and elevator faults. The non-linear and by parameterization timevarying model system is represented via an LPV polytope structure through the analysis of the linearised models. Additionally, the practical faults (actuator faults) are added onto the elevator and vertical velocity of the aircraft is assumed to be affected by the disturbance. Efficient interiorpoint algorithms and Linear Matrix Inequalities (LMI) are used within an H  framework, based on the well know formalism of Apkarian, Gahinet and Becker 1995. It is shown that the H   -suboptimal solution for the fault estimator on each vertex can be calculated through solving set of LMIs. The Machan model is linearised on predefined trim points using the Jacobian linearization approach. Based upon the analysis of the linearised models, the constant gain stabilising controller (both longitudinal and lateral dynamics) is calculated using eigenstructrue assignment. The polytopic

6680

10.3182/20110828-6-IT-1002.00565

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

LPV fault estimators are then calculated, at vertices of the polytope separately, to satisfy the given global robustness specification. The LPV fault estimator is deduced through time-varying system parameters, involving the polytope description. It is demonstrated that with online measurement of the aircraft parameters, the elevator faults can be constructed effectively and robustly.

The design of a polytopic estimator can be written as  u (t )  x f (t )  Af ( ) x f (t )  B f ( )    y p (t ) 

Section 2 introduces the theoretical foundations of the LPV estimator design. The nonlinear Machan aircraft case study is illustrated in Section 3. The concluding discussion is given in Section 4.

Therefore, the estimated error vector e f (t )  fˆ (t )  f (t )

2. ROBUST LPV ESTIMATION

the estimation of the fault

The LPV system with faults is given by (1)

D f ( ) are matrices with appropriate dimensions, to be

x p (t )  Ap ( ) x p (t )  B p ( )u (t )  E p ( ) d (t )  Fp ( ) f (t )

(1)

and x f (t )  R is the state vector of the estimator, fˆ (t ) is n

f (t ) . Af ( ), B f ( ), C f ( ) , and

designed. Define x pf (t ) and wudf (t ) to be [ x p (t ) x f (t )] and T

[u (t ) d (t ) f (t )] respectively. The LPV estimator (4) can be rewritten as:

where, x p (t )   , u (t )   , y p (t )   , and d (t ) q r

 R g is minimized. Here u(t ) and y p (t ) are defined in (1),

T

y p (t )  C p ( ) x p (t )  D p ( )u (t )  G p ( ) d (t )  H p ( ) f (t ) n

(4)

 u (t )  fˆ (t )  C f ( ) x f (t )  D f ( )    y p (t ) 

m

are the states, control inputs, outputs, and disturbances. g f (t )   is the fault vector.    is a varying parameter s

vector, and Ap ( ) , B p ( ) , C p ( ) , D p ( ) , E p ( ) , Fp ( ) , G p ( ) and H p ( ) are matrices with appropriate dimensions.

Assumptions that apply to system (1) are (Apkarian, Gahinet and Becker 1995):

A F ( ) :  C

f

( )

f

( )

  A ( ) F    Co   C ( )  ( )  i  1, , N 

B f ( )  Df

f

i

B f ( i ) 

f

i

D f ( i )  

i

  (5)  

The estimator system structure is shown in Fig.1, wherein two figures represent the equivalent structure information. The figure on right side of Fig.1 is well posed for the later H  filter design.

(A1) System (1) is stable. (A2) Vector  (t ) varies in a polytope  with vertices: s 1 ,2 , , N ( N  2 ), i.e.:

 (t )   : Co{1 , 2 ,

, N }

    ii :  i  0,  i 1 N

N

 i 1

i

  1  (2)

Fig. 1. Robust LPV Estimation System Structure Eq. (1) can now be re-written as:

(A3) As the state-space matrices depend affinely on  (t ) , system (1) is assumed to be polytopic, i.e.

x p  Ax p  B1 wudf  B2 fˆ

 Ap ( ) B p ( ) E p ( ) Fp ( )   C ( ) D ( ) G ( ) H ( )   p p p  p   Ap (i ) B p (i ) E p (i ) Fp (i )     , Co  C p ( i ) D p ( i ) G p ( i ) H p ( i )   i  1, , N    (A4) C p ( ), D p ( ), G p ( ) , and H p ( ) are parameter independent, i.e. C p (i )  C p , Dp (i )  Dp , G p (i )  G p , H p (i )  H p ,

e  C1 x p  D11wudf  D12 fˆ

(6)

y p  C2 x p  D12 wudf  D22 fˆ The estimator illustrated in Fig.1 can then be expressed as:

x pf (t )  A ( ) x pf (t )  B ( ) w udf (t ) e f (t )  C ( ) x pf (t )  D ( ) w udf (t )

(3) where:

i  1, , N

6681

(7)

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

A ( )  A0  BF ( )C B ( )  B0  BF ( )D 21 C ( )  D12F ( )C D ( )  D11  D12F ( )D 21 B 0 B2 A 0 A0  B0  1 B (8) 0 0 0 I 0 0 I 0 C  D12   0 D12  D 21  C2 0 D21 Here, Problem 1 is defined to solve the estimation problem, Problem 1 For the LPV system (1) with assumptions (A1)-(A4), design a polytopic LPV estimator (4), such that the L2 -induced

 

   

   

 S o Ap ( i )  ATp ( i ) S o  T B1 ( i ) S o  I   0 

N S  0 

0

T

S o B1 ( i )  I D11

 N S D11    0  I  0

0

  0 (11)

T

i  1,

I ,N

 Ro I   0  I S   o  Once the matrices R0 and S0 are obtained, the LPV estimator can be constructed as following Algorithm 1

norm of the operator mapping wudf (t ) into e f (t ) is bounded

Step1. Computing the full rank matrices M o , N o using

by a scalar number  for all parameter trajectories  (t ) in the polytope Θ.

SVD such that:

Considering the structure of (7) and according to the assumptions (A2)-(A4), it can be verified that the system (7) is polytopic, and the Lemma 1 can be used as an adaptation of the results from (Apkarian, Gahinet and Becker 1995). Lemma 1 For LPV system (7), the following statements are equivalent: (1) L2 -induced norm of the operator mapping wudf (t ) into

M o N o  I  Ro So T

(12)

Step 2. Computing X as the unique solution of the linear matrix equation:

 I Ro    So I  T   T   0 M o   No 0 

X

(13)

Step 3. Compute F ( i ) by solving the matrix inequality

e f (t ) is bounded by a scalar number  .

 ( i )  U x F ( i )V  V F ( i )U x  0 i  1,

(2) Parameter trajectories  (t ) in the polytope Θ,

where:

T

There exists X  X  0 satisfying the system of LMIs:

T

T

,N

(14)

T

 XA ( )  A ( ) X  B ( ) X   C ( )  T

i

i

T

i

i

XB ( i )  I

T

   

D ( i )  0 T

D ( i )

 XA0 ( i )  A0T ( i ) X XB0 ( i ) 0  T T  ( i )   B0 ( i ) X  I D11   0 D11  I   

C ( i )   I i  1,

(9)





U x  B X 0 D12 , V  C D21 0 

,N

The main result of this Section is stated in Theorem 1 which provides the solution to Problem 1, (proof omitted). Interested readers can refer to (Apkarian, Gahinet and Becker 1995). Theorem 1 Consider the LPV system in (1) with assumptions (A1)-(A3). Let N R and N S denote bases of the null spaces of

( B2T D12T ) and (C2 D21 ) , respectively. There exists a

T

T

(15)

(16)

Step 4. Solve the polytopic LPV estimator:

F ( ) 

r

  F ( ) i

p

i

i 1

(17)

Where,  i , i  1,..., N , are any solutions of the convex decomposition problem: N

    i i i 1

polytopic LPV estimator that can determine the solution of Problem

1

0  So  So  R T

if n n

the

 Ap ( i ) Ro  Ro Ap ( i ) 0  0  I  I  T T  B1 ( i ) D11 

0

T

T

n n

,

can be found such that: T

N R  0 

matrices 0  Ro  Ro  R B1 ( i ) 

N R  0   I  D11

i  1,

0

T

  0 (10)

I

,N

3. NONLINEAR MACHAN AIRCRAFT CASE STUDY The aircraft chosen in this work is an unmanned aircraft (UAV) or remotely piloted vehicle, the Machan, used as a development vehicle by Marconi Avionics, Chatham UK for research on high incidence flight with RAE Farnborough and NASA Dryden in the 1980’s (Aplin and Lamb, 1981). The Machan Euler equations relate the forces X , Y and Z and moments L , M and N in the aircraft body axes to the angular and linear velocities in the inertial axes are shown as

6682

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

m(u  qw  rv)  X m(v  ur  pw)  Y m( w  vp  qu )  Z I x p  ( I z  I y )rq  L I y q  ( I x  I z ) pr  M I z r  ( I y  I x )qr  N

aircraft is trimmed based upon the trajectory of true airspeed

vT  u  v  w and altitude h inside the flight envelope separately, the corresponding eigenvalue trajectory of the linearised models are shown in Fig.2 and Fig.3. An investigation has shown that the aircraft dynamics are more affected by variation in true air speed than by changes in altitude. 2

(18)

where, I x , I y , I z are the moments of inertia about the axes through the centre of gravity but parallel to the aircraft body axes. u , v and w are the forward, side and vertical velocity of the aircraft respectively. p , q and r are the roll rate, pitch rate and yaw rate, respectively.

2

2

pole migrations vs machan velocity 6

4

start 25m/s end 50m/s intermediate velocity points

2

X  X E  D cos   ( Lw  LT ) sin   mg sin  Y  Ya  mg cos  sin  Z  ( Lw  LT ) cos   D sin   mg cos cos L  LE  La M  M   Lw (cg  0.25)c  (lt  0.25  cg )Lq N  N

imaginary part

The aerodynamic force and moment equations are:

0

-2

-4

(19)

-6 -25

-20

-15

-10 real part

-5

0

5

Fig. 2. Eigenvalue Migrations vs Machan Velocity

 is the incidence,  and  are pitch angle and roll angle respectively, Ya is the side force, cg is the position of the aircraft centre of gravity, X E is the thrust force due to the engine, lt is the tail moment, D is the force acting on the

pole migrations vs machan height 4 3 2

imaginary part

airframe, Lw , LT and Lq represent the wing lift, total tail lift and tail lift due to the pitch rate respectively, M a , N a and La are the pitching moment, yawing moment and rolling

Ke

0 -1

-3 -4 -25

The engine (thrust) dynamic is given as XE 

1

-2

moment component respectively and LE is the rolling moment due to the engine.

( Pmax TH  p  X EU 2 )

start 20m end 120m intermediate height points

-20

-15

-10 real part

-5

0

5

Fig. 3. Eigenvalue Migrations vs Machan Altitude (20)

Fig.1 shows that for this system,

u (t )  u p (t )

and

where, Pmax , TH ,  p , K e and U 2 represent the maximum

wudf (t )  [utrim (t ) d (t ) f (t )] . u p (t ) and utrim (t ) represent the

engine power, the throttle demand, the propeller efficiency, engine rise rate and the air flow rate respectively. The details of the parameters are given in Aslin 1985.

plant inputs and the stabilizer inputs utilised for trimming purpose respectively. The plant inputs are the sum of the controller outputs and the stabilizer inputs. The controller and fault estimator are designed with fully measurable system states. The time-varying system with additive actuator faults can now be expressed as:

This nonlinear Machan model contains fourteen states: u , v ,

w , yaw angle  ,  ,  , q , altitude h , X E , elevator deflection  , roll rate p , yaw rate r , rudder deflection  , aileron deflection  . The widely adopted linearization approaches for control system design include: Jacobean linearization approach, state transformation approach and the function substitution approach (Marcos 2001). Here the Jacobian linearization approach is used for the purpose of Machan linearization, which requires multiple trim points to be chosen. However, it is well posed for polytopic LPV system design, since the polytopic system is affine and the scheduling parameters can be selected as trimmed variables. Therefore, the number of the scheduling parameters will be small, which greatly decreases the computational requirements. The nonlinear

T

x(t )  A(vt (t )) x(t )  B(vt (t ))u(t )  Fa f a (t )

(21)

where Fa is fault distribution matrix and f a represents actuator faults. The actuator fault estimate fˆa (t ) is implemented using Algorithm 1. The open-loop Machan is unstable and hence a stable closedloop system must be configured to satisfy (A1). A common constant gain eigenstructure assignment control (lateral and longitudinal) is constructed to stabilise the fault-free openloop system. Eigenstructure assignment is a widely adopted approach to stability augmentation and to achieve good handling quality modal de-coupling.

6683

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

Fault test signals - (fault 1) and oscillatory faults (fault 2) are added onto the elevator (Fig.4). Fault 1 depicts a multi-step up-step down signal (rad) which is useful for exciting nonlinearity and testing for fault delectability properties. Fault 2 (rad) is constructed through a multi-sine signal.

Longitudinal System Outputs

90

f 2 (t )  0.005(sin(10t )  2sin(t ))

altitude Xe vertical velocity forward velocity

80 70 60 50 40 30 20 10 0 -10

0

50

100

150

200

250

300

350

400

Time sec 0.4

Longitudinal System Outputs

0.03

Elevator Fault

0.02

0.01

0

-0.01

-0.02

-0.03

0.2 0.1 0 -0.1 -0.2 -0.3

0

50

100

150

200

250

300

350

pitch rate pitch angle elevator deflection

-0.4 -0.5 -0.6

0

50

100

150

200

250

300

350

400

Time sec

400

Fig.6. Longitudinal Output Responses with Elevator Fault 2

Time sec

From the longitudinal system outputs as shown in Figs.5 and 6, the altitude variation does not affect the longitudinal dynamics substantially. Hence, although the height could be used as an LPV scheduling parameter a more useful scheduling parameter is the total velocity or air speed chosen in the known speed interval vT  [32 36]. From this

0.03

0.02

Elevator Fault

0.3

0.01

0

-0.01

-0.02

-0.03 300

information the estimators on each polytope vertex Fi ( ) 310

320

330

340

350

360

370

380

390

400

Time sec

Fig.4. Faults 1and 2 acting on the elevator These fault levels are of practical significance in terms of their magnitude of less than 1 actuator degree. The longitudinal system states with elevator faults 1 and 2 are shown in Figs.5 and 6, respectively. The initial values of the longitudinal states are predefined. The lateral dynamics do not change much since the full nonlinear system is trimmed based on the routine of longitudinal states (details omitted).

i  1,2 can be calculated. After 41 iterations, the sub-optimal solutions  are 8.491 and 5.345 respectively, for each vertex. The LPV estimator is then built through the combination of Fi ( ) i  1,2. Fig.7 shows the results of the elevator fault estimation for both fault 1 and fault 2, with a zero-mean Gaussian disturbance d (t ) with variance 0.005 added onto w , the aircraft vertical velocity (considered as a simple wind gust effect). A suitable shaping filter is added behind the estimated fault, to establish a trade-off between improving the rise time and decreasing the steady state estimation error.

Longitudinal System Outputs

Elevator fault estimation rad

0.03 100

altitude Xe forward velocity vertical velocity

80

60

40

20

0

0.01 0 -0.01 -0.02 -0.03

-20

0

50

100

150

200

250

300

350

fault estimated fault

0.02

200

220

240

260

280

300

320

340

360

380

400

Time sec

400

Time sec 0.05

Estimated Fault Fault

0.04 0.03

Elevator Estimation

Longitudinal System Outputs

0.4 0.3 0.2 0.1 0 -0.1

-0.3

0 -0.01 -0.02

-0.04 -0.05 250

-0.4 -0.5

0.01

-0.03

elevator deflection ptich rate pitch angle

-0.2

0.02

260

270

280

290

300

310

320

330

340

350

Time sec 0

50

100

150

200

250

300

350

400

Fig.7. Fault estimation provided by polytopic LPV estimator

Time sec

Fig.5. Longitudinal Output Responses with Elevator Fault 1

Simulation results show that the robust LPV fault estimator provides a good online estimation performance for both faults 6684

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

1 and 2. The blue curve represents the fault signal and the red curve depicts the estimated fault signal. Fig. 8 also shows the change of aircraft altitude (40 to 35 m).

45

5. ACKNOWLEDGMENTS The authors thank the European Commission for research funding in the contract FP7-233815, Advanced Fault Diagnosis for Safer Flight Guidance and Control (ADDSAFE).

Altitude

40

REFERENCES

35

30 200

220

240

260

280

300

320

340

360

380

Apkarian, P., Gahinet, P. and Becker, G. (1995). "Self-scheduled H  control of Linear parameter-varying systems: a design example." Automatica 31: 1215-1261.

400

Time sec

Fig.7. Change of aircraft altitude at 300s Fig. 9 shows the corresponding OFC fault estimation illustrating the strong estimation recovery following the change of the aircraft altitude. 0.05

Estimated Fault Fault

0.04

Elevator Fault

0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 -0.05 250

300

350

400

Time sec

Fig.9. Fault estimation result due to the change of altitude

4. CONCLUDING DISCUSSION The principle of LPV design for robust fault estimation has been introduced through the use of a set of LMIs using efficient interior-point algorithms. The robust LPV fault estimator design is well posed for the fault detection and fault isolation of a nonlinear aircraft system. The fault magnitudes considered in this study are very low (less than 1 degree of actuator rotation). The LPV study shows that the robustness of the estimation error is very high (even during a manoeuvre), corresponding to the system control inputs, disturbances and faults, but also the structure of the robust LPV estimator can be modified online through the measurement of the varying parameters. The full force and moment nonlinear UAV Machan has been used in this study as it has very non-linear longitudinal dynamics and is correspondingly hard to control. The corresponding development of the LPV-based approach for directly estimating the faults of nonlinear aircraft system is clearly demonstrated. The simulation results show that incipent (hard to detect) faults can be estimated robustly and effectively using the polytope estimator.

Aplin, J. and Lamb, P. (1981). "The Machan Simulation", Marconi Avionics, FARL report no. 262/1332/TN/S. Aslin, P. (1985). Aircraft simulation and robust flight control system design. DPhil thesis, Department of Electronics, University of York. Bokor, J. and Balas, G. (2004). "Detection filter design for LPV systems - a geometric approach." Automatica 40: 511-518. Bokor, J. and Kulcsar, B. (2004). "Robust LPV detection filter design using geometric approach." Proc.12th Mediterranean Control Conference on Control and Automation. Casavola, A., Famularo, D., Famularo, D. and Sorbara, M. (2007). "A fault detection filter design method for linear parameter-varying systems." Proc.J.IMechE Part I: 865874. Casavola, A., Famularo, D., Franze, G. and Patton, R. J. (2008). A Fault detection filter design method for a class of linear time-varying systems. 16th Mediterranean conference on control and automation. Corsica. Chen, J. and Patton, R. J. (1999). Robust model-based fault diagnosis for dynamic systems, Kluwer. Ganguli, S. A., Marcos, A. and Balas, G. (2002). Reconfigurable LPV control design for Boeing 747-100/200 longitudinal axis. Proc. ACC2002. Anchorage: 3612-3617. Henry, D. and Zolghadri, A. (2005). "Design and analysis of robust residual generators for systems under feedback control." Automatica 41: 251-264. Marcos, A. (2001). A linear parameter varying model of the Boeing 747-100/200 longitudinal motion. Master thesis, Department of Aerospace and Engineering Mechanics, University of Minnesota. Marcos, A., Ganguli, S. and Balas, G. J. (2005). "An application of H  fault detection and isolation to a transport aircraft." Control Engineering Practice 13(1): 105-119. Patton, R. J., Chen, L. and Klinkhieo, S. (2010). LPV approach to friction estimation as a fault diagnosis problem. 15th International conference on methods and models in automation and robotics, Miedzyzdroje, Poland, Sept 2010 Patton, R. J., Putra, D. and Klinkhieo, S. (2010). Friction compensation as a fault-tolerant control Problem, Int. J. Sys. Sci., 41, (8), August 2010, pp987-1001 Weng, Z., Patton, R. J. and Cui, P. (2010). Robust fault estimation for linear parameter-varying time-delay systems. 16th Mediterranean conference on control and automation. Congress Center, Ajaccio, France. Wu, F. (2001). "A generalised LPV system analysis and control synthesis framework." Int. J. of Control 74: 745-749.

6685