Computational Questions of Equilibrium Calculation with Application to Nonlinear Aircraft

LiTH-ISY-R-1990

Computational Questions of Equilibrium Calculation with Application to Nonlinear Aircraft

Dynamics

Mats Jirstrand and Torkel Glad Department of Electrical Engineering Linkoping University, S-581 83 Linkoping, Sweden

www:

email:

1997-11-14

REGLERTEKNIK

AUTOMATIC CONTROL

LINKÖPING

Technical reports from the Automatic Control group in Linkoping are available as UNIX-

compressed Postscript les by anonymous ftp at the address

.

Computational Questions of

Equilibrium Calculation with Application to Nonlinear Aircraft Dynamics

M. Jirstrand, S. T. Glad

Department of Electrical Engineering

Linkoping University, S-581 83 Linkoping, Sweden

WWW: http://www.control.isy.liu.se

E-mail:

matsj,torkel

@isy.liu.se

Abstract

. The computation of possible stationary orientations of an aircraft subjected to constraints on actuators is treated using methods from computer algebra. The same techniques are also used to investigate possible turn maneuvers of a missile subjected to actuator constraints. The paper also gives a brief introduction to the utilized algebraic methods, cylindrical algebraic decomposition and quantier elimination.

1 Introduction

Nonlinear aircraft dynamics is usually described by a system of nonlinear dierential equations

_

x

=

(

)

(1)

where

contains variables describing airstream orientation, rotational velocities and possi- bly orientation. The control surface deections are represented by the control variables

. The function

contains the rigid body dynamics and the aerodynamic forces and mo- ments. The latter quantities are often given by tables of values initially measured in wind tunnel experiments and later on updated during real ight. For computational reasons these quantities are often represented by polynomial approximations. The question of

nding an equilibrium of the system under constraints on the actuators then results in a system of polynomial equations and inequalities. If restrictions are put on the linearized dynamics at the equilibrium, further polynomial constraints are produced.

1

Solving the problem purely by numerical methods is made complicated by the fact that multiple solutions might exist. The aim of this paper is to show how two symbolic methods implemented in computer algebra packages, known as quantier elimination and cylindrical algebraic decomposition, can be used to solve the problem.

The paper is organized as follows. In Section 2 we recapitulate some theory about aircraft dynamics. Section 3 gives a brief introduction to cylindrical algebraic decompo- sition and quantier elimination through a number of examples. The application of the algebraic methods to nonlinear aircraft dynamics is presented in Section 4. Here both an aircraft and a missile example are treated. Section 5 concludes the paper.

2 Aircraft Dynamics

Nonlinear aircraft dynamics for conventional aircrafts is usually described by the so called 6-DOF (degrees of freedom) equations 13], which is a classical state space description of a nonlinear system, written on the form (1). The state vector has 12 components and comprise the velocity vector, the vector of Euler angles, the angular rate vector, and the position vector of the aircraft all states is measured in some convenient coordinate system.

In this paper we consider the computation of possible stationary orientations (

) of the aircraft w.r.t. its wind xed coordinate system, see Figure 1. An investigation of

x-axis (body) x-axis (stability) α

β Relative wind

x-axis (wind) z-axis

(body) y-axis

(body)

Figure 1: The orientation of an aircraft with respect to the airow.

the state equations reveals that this problem is equivalent to keeping the angular rate components of the state vector constant and equal to zero. The derivatives of the angular rate components of the state vector are given by the so called moment equations of the

2

aircraft model

_ = (

1

+

2

)

+

3

 +

4

_

Q

=

5

6

(

²

² ) +

7

_

R

= (

8

2

)

+

4

 +

9

(2) where

and

are the components of the angular velocity vector of the body xed coordinate system. The

are constants determined by the moments and cross-products of inertia of the aircraft and 

, and

are the rolling, pitching, and yawing moment, respectively. From equation (2) we observe that the condition for the aircraft to have stationary orientation is 

=

=

= 0 for all values of

3

4

7

and

9 satisfying

c

3

9

24

= 0 and

7

= 0. The moments 

, and

are proportional to the dimensionless aerodynamic coecients

, and

, respectively, which usually are measured and given by look-up tables and interpolation routines 13]. These quantities are functions of the angle of attack

, the sideslip angle

, and the aileron, elevator, and rudder deections.

Now, if the aerodynamic coecients are approximated by polynomial functions the set of

and

values corresponding to a stationary orientation of the aircraft can be described by the solution set of the following system of polynomial equations and inequalities

C

l

(

) = 0

(

) = 0

(

) = 0

ju

i

j

1

= 1

2

3

(3)

where the inequalities are normalized constraints for the aileron, elevator, and rudder deections, respectively.

Using tools from computer algebra such as cylindrical algebraic decomposition and quantier elimination we can now eliminate the

from (3) and obtain closed form ex- pressions in

and

describing the possible stationary orientations of the aircraft.

3 Algebraic Tools

Cylindrical algebraic decomposition (CAD) is a method for nding feasible solutions of a system of polynomial equations and inequalities, a so called real polynomial system.

A formula is an extension of real polynomial systems to expressions involving Boolean operators (

) and quantiers (

). There are algorithms to perform quantier elimination (QE) in formulas, i.e., to derive equivalent expressions without any quantied variables.

3

3.1 Cylindrical Algebraic Decomposition

Given a set of multivariate polynomials the CAD algorithm decomposes R

into compo- nents, called cells, over which the polynomials have constant signs. The algorithm also provides a point, called sample point, in each cell that can be used to determine the sign of the polynomials in the cell. This decomposition can then be utilized to decide if any real polynomial system, determined by the given polynomials, has a solution. It suces to evaluate the polynomials over the sample points to decide whether there are any solutions or not. We do not intend to give a complete description of the method but to give a brief introduction to the method by the following example.

Example 3.1 Consider the following system of polynomial inequalities

x

22

2

1

2 +

₄₁

0

(2431

1

3301)

2

2431

1 + 2685

0

How do we determine a feasible point of the system?

In Figure 2 a CAD induced by

1 and

2 is shown. The cells (regions where

1 and

2

have constant sign) of the CAD are all white patches, arcs, vertical lines, and dots. There are 63 cells in this CAD, hence the algorithm also provides us with 63 sample points, one from each cell. If we evaluate the polynomials over the sample points only three of the points satises the system of inequalities.

f

1

f

2

f

2

Figure 2: A CAD induced by

1 and

2 and the feasible sample point

1 = ¹ ₂

2 = ³ ₄ . The point

1 = ¹ ₂

2 = ³ ₄ , also shown in Figure 2 is one of the sample points and since

f

1

0,

2

0 at this point it is also a feasible point to the system of inequalities.

The CAD algorithm can be divided into three phases: projection, base and extension.

The projection phase consists of a number of steps, each in which new sets of polynomials

4

is constructed. The zero sets of the resulting polynomials of each step is the projection of

\signicant"points of the zero set of the preceding polynomials, i.e., self-crossings, isolated points, vertical tangent points etc. In each step the number of variables is decreased by one and hence the projection phase consists of

1 steps.

The base phase consists of isolation of real roots,

R ¹ of the monovariate polyno- mials, which are the outputs from the projection phase. Each root and one point in the each interval between two roots are chosen as sample points of a decomposition of R ¹ .

The purpose of the extension phase is to construct sample points of all cells of the CAD of R

. This is done by a lifting procedure that consists of

1 steps and results in a CAD for R

^.

Observe that usually it is easy to visualize zero sets of polynomials in two variables and there is no need for a CAD of R ² to solve the problem. However, for polynomials in more than two variables there are no easy way of visualizing the zero sets and then the method really shows its power.

CAD was invented in the mid seventies by Collins 3] and since then many improve- ments of the original algorithm has been made, e.g., 8]. The algorithm has a high com- putational complexity 12, 2] and ecient implementations is therefore very important.

Introductions to CAD can be found in 9, 10, 1].

3.2 Quantier Elimination

QE is a method for eliminating variables, quantied by

(there exists) and

(for all), from formulas consisting of polynomial equations and inequalities combined by Boolean operators

(and),

(or), and

(implies). Hence, the objective is to nd an equivalent formula without quantiers in the unquantied, or free, variables only.

CAD or a partially built CAD is usually used as a subalgorithm in existing implemen- tations of QE algorithms 4]. Once a CAD is given, elimination of quantied variables can be performed by simple sign evaluation of the projection polynomials that describes the cells together with truth tests of formulas. This method is also illustrated by an example.

Example 3.2 Consider the following QE formula (

2 ) 

22

2

1

2 +

41

0

(2431

1

3301)

2

2431

1 + 2685

0 ]

This formula describes the projection of the solution set of the system of inequalities onto the

1 axis. Hence, the result after quantier elimination is one or several intervals for

1 .

5

Utilizing the CAD in Example 3.1 an equivalent formula in

1 can be computed 13 5

¹

1

In this paper we have used

4] developed by Hong et al. at RISC in Austria to perform QE in the nontrivial examples. Another program for QE in expressions where the quantied variables appears with low degree is

5], which seems very promising.

4 Applications to Aircraft Dynamics

To demonstrate the application of quantier elimination in computations regarding air- craft dynamics we will treat two examples. The rst example is the computation of equilibrium states of the F-16 ghter aircraft and the second example deals with the ma- neuvering of a missile. In both cases the dynamics is subjected to constraints on the control surface deections.

Example 4.1 The following expressions are scaled versions of the aerodynamical co- ecients of an F-16 ghter. The polynomial approximations have been obtained from look-up tables and interpolation routines in 13].

C

l

=

38

170

+ 148

+ 4

+

(

52

2

+ 114

79

+ 7

+ 14

)

+

(14

10

+ 37

48

+ 8

13

13

+ 20

+ 11

)

C

M

=

12

125

+

+ 6

+ 95

21

+ 17

;

202

+ 81

+ 139

C

N

= 139

112

388

+ 215

38

+ 185

+

(

11 + 35

22

+ 5

+ 10

17

)

+

(

44 + 3

63

+ 34

+ 142

+ 63

54

69

26

)

Here

1 and

2 are normalized angle of attack and sideslip angle, respectively and

1 ,

2 , and

3 the aileron, elevator, and rudder deections.

The angles of attack and sideslip angles corresponding to constant orientation of the aircraft can now be described by a formula using quantiers and Boolean operators to- gether with (3) as follows

9u

1

2

3

h

C

l

= 0

= 0

= 0

²

1

Eliminating the quantied variables gives



1

0

2

0]



3

0

2

0]]



4

0

5

0

6

0

7

0]

(4)

6

where

p

1

31

21

1

2

31

21

1

p

3

61

51

41

31

114876460x

21

1

and

4

7 are polynomials in

1 and

2 of total degree 7 of about 40 terms each.

However, although the output is quite complicated the elimination of the control variables reduce the dimension of the problem and the two dimensional solution set (4) can easily be visualized. The solution set is shown in Figure 3.

0 0.2 0.4 0.6 0.8 1

x1 -1

-0.5 0 0.5 1

x2

Figure 3: The possible stationary orientations if

1 (white).

In Figure 3 we observe that for zero sideslip angle there are no constraints (in the valid region of approximation) on the angle of attack, i.e., the aircraft can be kept stationary at any angle of attack between 0 and 1 (normalized values). Furthermore, the maximal constant sideslip angle can be obtained at an angle of attack of approximately 0

4.

A number of other questions related to equilibrium calculations can also be posed.

What happens if some actuators are further constrained, e.g., jammed or damaged? In which regions of the

1

2 -plane can we guarantee that certain moments can be generated by the control surfaces? How does the computed equilibrium region depend on a design parameter or an uncertainty?

As an example we treat the following question. Suppose that the deection of the ele- vator has been further constrained to

2

1

2 . For what angles of attack can the elevator generate pitching moments corresponding to an aerodynamic coecient of magnitude greater than 40? The QE formula becomes

9u

2 

²

1600

4

² ₂

1 ]

7

Elimination

2 in this formula gives 0

1

, where

0

739182 is the unique root between 0 and 1 of

556

³

646

² + 355

134 = 0

The answer to this question is a region (in this case an interval) where we have some guaranteed robustness against disturbances, since we can generate compensating moments of a prescribed magnitude.

Note also that it is easy to incorporate symbolic design parameters or parameters representing system uncertainties in the above formulas and investigate their eect on the computed regions, since QE is a symbolic method.

Example 4.2 Here we consider the aerodynamics for pitch-axis control of a missile, see Figure 4.

x

1

u

1

C

Z C

M

Figure 4: Normal force and moment on a missile.

The quantities that determine the motion in the vertical plane of the missile are the normal force and pitching moment, which in turn are functions of the angle of attack, sideslip angle, and elevator. The following is a polynomial expression for the normal force (neglecting gravity)

C

Z

(

1

2

1 ) =

13

1

120

₂₁ + 390

₃₁

100

1

22

1350

₅₁

1400

₃₁

₂₂ + 875

1

42

4

1

and the pitching moment

C

M

(

1

2

1 ) =

13

1

10

₂₁ + 4

₂₂ + 460

₃₁

520

1

22

2600

₅₁ + 3000

₃₁

₂₂ + 680

1

42

35

1

where

1 is the angle of attack,

2 is the sideslip angle, and

1 is the elevator angle. These polynomial approximations comes from 6].

8

An analysis of how quick the missile can turn in the vertical plane can be carried out by QE methods. The missile describes a circular motion if the normal force is constant and the pitching moment is zero, the larger force the smaller radius. This implies that limits on what normal forces that can be generated limit the possible turn radii. The question how quick the missile can turn, subjected to elevator deection constraints, can thus be formulated as follows: What normal forces can be generated, for constant speed of rotation (i.e., zero pitching moment), if the maximum elevator angle is limited? The analysis is carried out in the vertical plane so the sideslip angle is

2 = 0.

Adding some constraints on the angle of attack, elevator, and normal force gives a quantier formula in the variables

,

1 ,

1

9x

1

1 

(

1

0

1 ) =

(

1

0

1 ) = 0

0

1

3

10

1

40

1

1

20

8

0 ]

Quantier elimination via cylindrical algebraic decomposition shows that the

-axis is divided into 51 cells. An investigation of the truth value of the above formula over each cell (interval) shows that there exist

1 and

1 satisfying the formula for

in the intervals



1

2 ]



3

0]

where

1 =

7

5993,

2

3

148712 is the unique real root between

13

4 and

25

8 of 38932892876800000

+ 1718497027686400000

+24272290834718720000

+ 138926900972696832000

+326132672903058884000

+ 250374264550106371777 = 0

and

3

1

504850 is the unique real root between

13

8 and

3

2 of

38932892876800000

+ 1718497027686400000

+24272290834718720000

+ 138926900972696832000

+326132672903058884000

+ 250374264550106371777 = 0

The intervals correspond to two disjoint \turn regions"shown to the left in Figure 5.

To the right in Figure 5 the moment and force functions for

1 = 0 are shown. Here we can clearly see that the constraint on the elevator makes it impossible to compensate for large negative moments. This is the reason for the limits on the turn radius. The region of smaller turn radii in Figure 5 is due to the shape of the moment curve. However, it is probably hard to control the missile to perform turns corresponding to the region of smaller radii due to the intermediate region of non-realizable turns.

9

-4 -3 -2 -1 0 1

0 0.1 0.2 0.3 0.4

-12 -10 -8 -6 -4 -2 0

0 0.1 0.2 0.3 0.4

x

1

x

1

M

C

Z

Figure 5: Left: Possible turns. Right: The normal force and pitching moment (

1 = 0).

5 Conclusions

We have given a brief introduction to two symbolic methods, cylindrical algebraic decom- position and quantier elimination, and have shown how these methods can be utilized in equilibrium calculations of nonlinear aircraft dynamics. The possible stationary orienta- tions of an aircraft w.r.t. a wind xed coordinate system were studied where constraints on control surface deections also were taken into account. A second application was the calculation of what turn radii a missile can obtain in the vertical plane if constraints on the deection of the elevator are present.

The symbolic nature of the methods makes it easy to extend the computations to include design parameters and parameters that describe uncertainties. The symbolic form of the result of a computation, in terms of parameters to be determined, often facilitates an optimal choice of such parameters.

The main limitation of the quantier elimination method is its algorithmic complexity which at present limits its usefulness for large scale applications. The development of symbolic algorithms and existing implementations have not yet been fully recognized by the control community. Many problems in control can be treated successfully with symbolic methods and the interested reader is referred to 9, 11, 7] for further applications.

Acknowledgement

This work was supported by the Swedish Research Council for Engineering Sciences (TFR), which is gratefully acknowledged. The authors are also grateful to Dr. Hong at RISC, Austria for providing them with a copy of the

program.

10

References

1] D.S. Arnon, G.E. Collins, and S. McCallum. Cylindrical algebraic decomposition I:

The basic algorithm. SIAM J. Comput., 13(4):865{877, November 1984.

2] S. Basu, R. Pollack, and M.-F. Roy. On the combinatorial and algebraic complexity of quantier elimination. In S. Goldwasser, editor, Proc. 35th Annual Symposium on Foundations of Computer Science, pages 632{641, Santa Fe, NM, USA, 1994.

3] G.E. Collins. Quantier elimination for real closed elds by cylindrical algebraic de- composition. In Second GI Conf. Automata Theory and Formal Languages, Kaiser- slauten, volume 33 of Lecture Notes Comp. Sci., pages 134{183. Springer, 1975.

4] G.E. Collins and H. Hong. Partial cylindrical algebraic decomposition for quantier elimination. J. Symbolic Comput., 12:299{328, September 1991.

5] A. Dolzmann and Th. Sturm. REDLOG computer algebra meets computer logic.

Technical report, FMI, Universitat Passau, D-94030 Passau, Germany, February 1996.

6] P.G. Eriksson and T. Finnstrom. Exact linearization and analysis of non-linear systems. Master's thesis LiTH-ISY-EX-1199, Department of Electrical Engineering, Linkoping University, 1992.

7] K. Forsman. Constructive Commutative Algebra in Nonlinear Control Theory. PhD thesis, Dept. of Electrical Engineering, Linkoping University, S-581 83 Linkoping, Sweden, 1991.

8] H. Hong. Improvements in CAD-based Quantier Elimination. PhD thesis, The Ohio State University, 1992.

9] M. Jirstrand. Algebraic Methods for Modeling and Design in Control. Licentiate thesis LIU-TEK-LIC-1996:05, Dept. of Electrical Engineering, Linkoping University, Sweden, March 1996.

10] B. Mishra. Algorithmic Algebra. Texts and Monographs in Computer Science.

Springer-Verlag, 1993.

11] D. Nesic. Dead-Beat Control for Polynomial Systems. PhD thesis, The Australian National University, 1996.

12] J. Renegar. On the computational complexity and geometry of the rst-order theory of the reals. Parts I{III. J. Symbolic Comput., 13:255{352, 1992.

13] B. L. Stevens and F. L. Lewis. Aircraft Control and Simulation. Wiley, 1992.

11