A controller design method which applies to time varying linear systems

MISSING IMAGE

Material Information

Title:
A controller design method which applies to time varying linear systems
Physical Description:
Book
Creator:
Koenig, Kurt Walter, 1967-
Publication Date:

Record Information

Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 20381782
oclc - 32490049
System ID:
AA00022793:00001

Table of Contents
    Title Page
        Page i
    Table of Contents
        Page ii
        Page iii
    Abstract
        Page iv
        Page v
    Chapter 1. Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
    Chapter 2. Theory behind the design method
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
    Chapter 3. A time varying second order example
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
    Chapter 4. Derivation of a model of the Emraat missile
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
    Chapter 5. The dependence of gains on flight parameters
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
    Chapter 6. Computing look-up tables
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
    Chapter 7. Gain scheduling
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
    Chapter 8. Nonlinear simulations
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
    Chapter 9. Conclusion
        Page 104
        Page 105
    Appendix A. Aerodynamic data for the Emraat airframe
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
    Appendix B. Inertial data for the Emraat airframe
        Page 117
    References
        Page 118
        Page 119
    Biographical sketch
        Page 120
        Page 121
        Page 122
Full Text










A CONTROLLER DESIGN


METHOD WHICH APPLIES TO TIME VARYING
LINEAR SYSTEMS


By

KURT WALTER KOENIG


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


1994













TABLE OF CONTENTS




A B ST R A CT . . . . . . . . . . . . . . . . . .

CHAPTERS

1 INTRODUCTION ...............................

1.1 Earlier W orks . . . . . . . . . . . . . . . .
1.2 P purpose . . . . . . . . . . . . . . . . .

2 THEORY BEHIND THE DESIGN METHOD ..................

2.1 A Geometric Inerpretation of Lyapunov's Linear Stability Theorem .
2.2 Adding Control to the System ......................
2.3 A Linear Feedback Set to Control xTp ................
2.4 One Linear Feedback Matrix to Control xTpk . . . . . . .
3 A TIME VARYING SECOND ORDER EXAMPLE . . . . . . .

4 DERIVATION OF A MODEL OF THE EMRAAT MISSILE ........

4.1 The Nonlinear Model ...........................
4.2 The Linear M odel .............................

5 THE DEPENDENCE OF GAINS ON FLIGHT PARAMETERS .....


Generating p and V from M and Q . .
The Flight Parameter Generator . . .
Initializing the Iterative Lyapunov Design
The Iterative Lyapunov Design Method .
Formulation of a State Tracker . . .
Comparing Gains with Flight Parameters


Method .......


6 COMPUTING LOOK-UP TABLES ......................
6.1 Determining a Grid of Points ......................
6.2 Formulation of :he Design Constraints ...................
6.3 Generating the Look-Up Table ......................









7 GAIN SCHEDULING .....................


7.1 C urve Fitting . . . . . . . . . . . . . . . 83
7.2 Testing the Fit . . . . . . . . . . . . . . .. .. 83

8 NONLINEAR SIMULATIONS . . . . . . . . . . ... .. 88

8.1 The Nonlinear Simulation . . . . . . . . . . ... .. 88
8.2 A Test of State Tracking . . . . . . . . . . . .. .. 90
8.3 Simulation of Flight Scenarios ... ..... . . . . . . . . .. 93

9 CONCLUSION . . . . . . . . . . . . . . . ... .. 104


APPENDICES

A AERODYNAMIC DATA FOR THE EMRAAT AIRFRAME ........ 106


B INERTIAL DATA FOR THE EMRAAT AIRFRAME ........... 117


REFERENCES ................................... 118

BIOGRAPHICAL SKETCH ............................. 120


. . . 83













Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
A CONTROLLER DESIGN METHOD WHICH APPLIES TO TIME VARYING
LINEAR SYSTEMS
By
Kurt Walter Koenig
August 1994
Chairman: Dr. Thomas E. Bullock
Major Department: Electrical Engineering

A feedback controller design method has been formulated which applies to linear

time varying systems. The motivation behind this technique is the design of an

autopilot for an air-to-air missile. The missile under study is the extended medium

'range air-to-air technology (EMRAAT) missile, a theoretical bank-to-turn missile

which is under study by the United States Air Force. Knowledge gained by this

missile will be applied to future bank-to-turn missiles.

Conventional autopilot design techniques use pole placement, linear quadratic

regulators, linear quadratic Gaussian/loop transfer recovery techniques, or eigen-

structure assignment methods. These techniques are applied to a linear model which

depends on time varying flight parameters. In applying these methods, the false as-

sumption is made that the linear model is slowly varying. Since the system is changing

rapidly, no theoretical basis exists to support the success of the resulting controller.

The proposed design method finds state feedback gains which cause the closed loop

linearized system to be stable with respect to a specified Lyapunov function. Even

though flight parameters change rapidly, local stability around the operating point is

achieved.









The design algorithm finds feedback gains which place the eigenvalues of the

derivative of a given Lyapunov function. Limitations on placement of these eigenval-

ues are stated. This concept is applied to a second order linear time varying system

where ordinary pole placement techniques fail. The design method is then applied

to a linearized model of the EMRAAT missile which is a function of time varying

flight parameters. Feedback gains are generated as a function of these flight param-

eters. It is discovered that the gains depend on dynamic pressure, mach number,

and angle-of-attack. A combination of interpolation and polynomial fitting is used

to create a look-up table for the gains. The resulting controller is programmed into

a nonlinear simulation which runs missile and target scenarios. Small miss distances

are achieved.













CHAPTER 1
INTRODUCTION

The original goal of this dissertation was to formulate an autopilot design method

which stabilizes the flight of a bank-to-turn air-to-air missile. The result was the de-

velopment of a controller design method which applies to time varying linear systems.

When used on a nonlinear system, local stability is achieved. The remainder of this

chapter gives a brief description of the missile used in this study, a summary of related

works preceding this study, and the purpose and outline of this dissertation.

An air-to-air missile is launched from a military aircraft with the intent of in-

tercepting an enemy aircraft or target. Bank-to-turn (BTT) missiles have wings

giving them greater maneuverability than the conventional skid-to-turn (STT) mis-

siles which have fins. Thus, the target must work much harder to evade a BTT

missile. The dynamics of BTT airframes are highly unstable, making for a difficult

controls problem.

The missile under study in this paper is the extended medium range air-to-air

technology (EMRAAT) missile as shown in Figure 1.1 [1]. The EMRAAT missile

is a theoretical missile formulated by the Air Force with the intent of studying the

feasibility of the BTT concept. The knowledge gained from the study of this missile

will be used by the Air Force if it ever decides to make a real BTT missile.

The EMRAAT missile is equipped with a seeker. If the seeker is infrared based,
then it measures the line of sight angles to the target. In the case that the seeker is

radar based, then, in addition to the line of sight angles, the range and range rate

of the target are measured. This information is passed to the guidance law which

determines the desired normal accelerations needed to intercept the target. The



















cross section, constant -


Start -
wedge
fairings


7.5 din


All dimensions are in inches


V A-A
View A-A


Figure 1.1. The EMRAAT missile [1].


ZZ










autopilot attempts to achieve these accelerations by applying the proper control to

the control surfaces.

1.1 Earlier Works

Early autopilot designs separated the missile dynamics into two or three lower

order models. These models were linearized, and classical techniques were used to

stabilize them. Other designs applied more advanced techniques to higher order lin-

earized models. Some of these techniques are pole placement, linear quadratic regula-

tors, linear quadratic Gaussian/loop transfer recovery techniques and eigenstructure

assignment methods [2] Feedback controllers would be designed as a function of the

linearized models. Gains would then be scheduled against the flight parameters on

which the models depend. Gain scheduling is a popular technique used in control

designs and has been studied in depth by Shamma, Athans, and Cloutier [3,4,5,6].

The problem with these design methods is that they assume that the system is
time invariant or slowly varying. The linearized models depend on rapidly changing

flight parameters. Although these methods often work, no theoretical basis exists

to support the success of these designs. More recently, H, and yi synthesis design

techniques have been used to formulate one dynamic controller which would stabilize

the system for all modeled uncertainties [7,8]. A disadvantage of H.. controllers is

that they have very high orders. Also, the existence of a robust H"" controller which

satisfies the performance requirements is not guaranteed. The only way to determine

the stability of a given design is to test it using a nonlinear simulation.

Desoer [9] states if a given time varying system

x= A(t)x (1.1)

is stable for each t when t is frozen, then it is not necessarily stable when t is allowed

to change. However, (1.1) is stable if A(t) changes slowly. A conservative upper












Alpha vs. Time
12, ,,
12o ---i--" --------'- --- ---------
10-

8

S6
12
4

2

0-


0 1 2 3 4 5 6 7 8
Seconds

Figure 1.2. Angle of attack from one trajectory of the EMRAAT missile.

bound on supt>o IIA(t)lI is given which, if enforced, guarantees asymptotic stability

of (1.1). As will now be shown, this limit is too restrictive for the EMRAAT missile.

Figure 1.2 shows the angle of attack of the trajectory of a missile flown using an

already existing autopilot. Desoer proposes the Lyapunov function


V(x, t) = xT(61I + P(t))x (1.2)


where ei > 0 and P(t) is chosen so that


AT(t)P(t) + P(t)A(t) = -31 (1.3)


Desoer gives a bound on V,

V < xTx[-3 + 2elaM + aM3m4 (1.4)


where the following definitions are made:


aM :=sup IA(t)|| < c (1.5)
t>0










aoo is positive and

Re Ai[A(t)] < -2c7o Vi, Vt > 0 (1.6)

m is a constant and depends only on ao and aM so that

|exp(rA(t))|l mexp(-oCor), Vr > 0 Vt > 0 (1.7)

Also,

iaM:=sup IIA(t)l. (1.8)
t>0
If c is allowed to be very small, then from (1.4), stability would result if
4o.2
M < -- (1.9)

For the linear pitch model of the given trajectory, o0 = 15.75. At .04 seconds into
the trajectory the closed loop matrix is

A -3.2533 .8997 (1.10)
-1678.2 -83.50
At this instant in time, m = 22.88 and aM = 2231. 42 = .0036. From this we see

'that the linear system is changing much faster than the limit shown in (1.9). The
system is not slowly varying. Note that m was only found for t = .04 sec. and not
for all time. If a greater m were found for the rest of the trajectory, then the limit in
1.9 would be even more restrictive and the result would be the same. This test does
not show that the system is unstable. In fact, a Lyapunov function is known which
shows that this linear system is stable but Desoer's theory can not support this fact.
Wilson, Cloutier, and Yedavalli [1] give a stability condition for a constant linear
system with time varying uncertainties. Given a system,

x = [A+ E(t)]x (1.11)

if a positive definite P can be found which solves the Lyapunov equation,


PA + ATP +21 = 0


(1.12)











Alpha vs. Time
12


10


8-

a)
6-




4-/
0



& 3



0 1 2 3 4 5 6
Seconds

Figure 1.3. Angle of attack from one trajectory of the EMRAAT missile.

then the system will be asymptotically stable if


,max[E(t)] <,, 1 Vt (1.13)
7maz.(P)

where 0,max is the maximum singular value. This idea is initially interesting because

it requires the error, E(t), to be bounded, but does not constrain the time variation.

The above result can not be easily applied to the EMRAAT missile. This will be

shown by applying the result to a trajectory of the EMRAAT missile whose angle of

attack is shown in figure 1.3. In order to check the stability of the EMRAAT missile,

A must be chosen since it may be any matrix which is stable. It is believed that the

best choice of A should come from a closed loop matrix in the trajectory. With this

in mind, A was taken at t = 2 sec into the trajectory.

A-B '-4.477 .9526 (1.14)
A- BK = -1.684 -79.47 (











Sigma max[E(t)] vs. Time
2 5 0,,


200-


Fiso

E
150-
*a-
E
05100-


50



0 1 2 3 4 5 6
Seconds

Figure 1.4. A plot of axB[E(t)] versus time.

E(t) is the result of subtracting this matrix from the closed loop system from the

rest of the trajectory. Solving for P in (1.12) gives

[ 17.29 -.0454 (1.15)
S -.0454 .012 (115

Computing the upper bound in (1.13) gives

S= .0578 (1.16)
am)axj4P)
Figure 1.4 shows the plot of crmax[E(t)]. Obviously since marn[E(t)] is greater than

.0587 for most of the trajectory then (1.13) is not satisfied. This procedure was

repeated by taking the constant closed loop matrix from every point in the trajectory,

and the result was the same. The time varying linear system from this trajectory is

known to be stable because a Lyapunov function exists which can show this. The

upper bound in (1.13) can not be easily used if at all to support the claim of stability.

Desoer [9] gives a stability limit on the rate of variation of a given closed loop

time varying system. Wilson, Cloutier, and Yedavalli [1] give a stability limit on










the size of a time varying uncertainty. Both of these tests were found to be too

conservative when applied to an already existing stable autopilot, and the results
were inconclusive. In addition, these tests are analytic tools. As yet, no design
procedure exists which can give a closed loop linear system that will pass either
test. It would be desirable to formulate a design method that can yield a closed
loop time varying linear system which satisfies a less conservative stability condition.
Vidyasagar [10] gives two important results. The first can also be found in Hahn [11]

and says that given a time varying system


S= A(t)x (1.17)

if a positive definite matrix P can be found so that the matrix


PA(t) + AT(t)P (1.18)

is negative definite for all time, then the system is asymptotically stable.
The second result allows the first result to be applied to a nonlinear system. Given

a nonlinear system of the form

S= f(t,x)

where

f(t,O) = 0

and f is continuously differentiable, then let

A(t) = 8f(t, x)
ox x=O
and assume that
IIf(t,x) A(t)xl
lim sup = 0
llxll-o t>0 ||x










and A is bounded. If 0 is an asymptotically stable equilibrium point of i = A(t)z for
all time, then 0 is a locally stable equilibrium point for the system

5 = f(t,x)

These ideas can be applied in the following way. Given a nonlinear system

x= f(x,w,u) (1.19)

define

A(x,w) = B(x,w)= (1.20)
xw .(XW)nomna (X,W)nominai
where x is a vector containing the states and w is a vector containing additional
system parameters. At regions near the operating point, the system becomes close
to

Ak = A(x, w)Ax + B(x, w)Au (1.21)

where Ax and Au are small perturbations between the states and inputs and the

operating point. We would like to find a feedback control law,

u = -K(x,w)x (1.22)

so that

P[A(x, w) B(x, w)K(x, w)] + [A(x, w) B(x, w)K(x, w)]TP (1.23)

is negative definite for all x and w where P is positive definite. If such a control law
is found, then the perturbations from the nominal trajectory are locally stable for

the system

= f(x,w, -K(x,w)x) (1.24)

Shahruz and Behtash [12] give one control law which places some of the eigenvalues
of (1.23) for the case where P = I. However unnecessary limitations on where the










eigenvalues can be placed exist. Also, while many feedback control laws make (1.23)

negative definite, Shahruz and Behtash give only a small subset of these control laws.

1.2 Purpose

This paper gives an algorithm which can stabilize a linear time varying system by

placing the eigenvalues of (1.23) between desired bounds and has fewer limitations

than the control law in Shahruz and Behtash [12]. Limitations on the placement of

these bounds are stated. This algorithm applies to all systems for which a constant

positive definite P exists so that (1.23) can be made negative definite for all time.

First, theory is presented leading to the formulation of this algorithm. Then, to

demonstrate the usefulness of this algorithm, it is applied to a linear time varying

system where normal pole placement techniques fail. This algorithm is then applied

to the EMRAAT missile. A nonlinear model is made, and from this, a time vary-

ing linearized model is generated which is a function of several flight parameters.

Gains are computed, and their dependence on flight parameters determined. A gain

scheduling scheme is then implemented yielding local stability around the operating

point. The resulting design is tested in a nonlinear simulation. Finally, the results of

this test are given, and the usefulness of this design technique is evaluated.












CHAPTER 2
THEORY BEHIND THE DESIGN METHOD

The material in this section gives the theory leading to the proposed controller
design method. The first section presents a geometric interpretation of Lyapunov's
linear stability theorem. The effects of control on the velocity field of a system is
studied in the next section. The nonemptiness and convexity of the set of all feedback
gains which stabilize a system with respect to a given Lyapunov function are then
discussed. Finally, an iterative procedure that finds one element of this set is given.

2.1 A Geometric Interpretation of Lvapunov's Linear Stability Theorem

To understand why time invariance is a necessary assumption for eigenstructure
design methods the following second order linear time varying system was studied.
The example given now is from Vidyasagar's example 5.3,109 [10] and can also be
found in Khalil [13]. Given the following system

S= A(t)x (2.1)

where
[ -1 + a cos2(t) 1 a sin(t)cos(t)
A(t) -1 a sin(t) cost) -1 + a sin2(t) (2.2)
Vidyasagar [10] notes that the transition matrix is given by
P(t0) e(a-1)tcos(t) e-tsin(t) 1
q^(t,0)-= --e(a-1)sin(t) e-* cos(t) (2.3)

and the characteristic equation is

A2+(2-a)A+(2-a)=0 (2.4)

The roots of (2.4) have negative real parts for 1 < a < 2. The exponents in the first
column of D(t, 0) indicate, however, that the system is unstable for these values of a.
11











b
1.^ *- z Y \\


1r\
.-0.5 _\




-10 0 5 10 41 -1 -. 0 0-5 1
X- X1

Figure 2.1. The trajectory of (a) the system (2.1) and (b) a frozen system.

This system was simulated using MATLAB with a = 1.5 and xo = [1, 0]T. The

eigenvalues of A(t) are -!+ s7j and -- j Figure 2.la shows the state trajectory

of the system 2.1 where X, is assigned to the horizontal axis, and x2 is assigned to

the vertical axis. This plot demonstrates the instability of the system. If the system

were frozen, i.e. A(t) = A(0), then the stable trajectory of Figure 2.1b would result.

This trajectory is shown with the velocity vector field of the frozen system.

To explain the instability in Figure 2.la, we will now discuss the time varying

velocity field of equation (2.1). The plots in Figure 2.2 show the trajectory of the

time varying system in the state plane for eight instances in time. For each plot, the

velocity vector field for that instant in time is superimposed on to the trajectory. Each

plot in Figure 2.2 shows the boundaries of four pie-shaped regions which this paper

refers to as positive and negative regions. Two positive regions exist that are defined

as the set of all points whose velocities have an outer radial component, i.e. each of

these velocities have a component pointing directly away from the origin. Likewise,

there are two negative regions containing velocities with inner radial components.

The plots show that the boundaries separating the positive and negative regions

rotate in a clockwise direction. Also, the current position in the trajectory remains










inside one of the positive regions. Since the velocity of the system always has an

outer radial component, the magnitude of the state vector is always increasing so

that the system (2.1) is unstable.

A direct relationship exists between the positive and negative region boundaries

and the eigenvectors of the system. In this case, the eigenvectors are rotating clock-

wise at the same angular velocity as the region boundaries. It can be shown that if a

system has constant eigenvectors with time varying eigenvalues which have negative

real parts for all time then the system is stable. The stability question, however, is

not as easy to answer for the case with moving eigenvectors.

If the system has no positive regions, then it will be stable regardless of how

quickly it changes. However the converse is not true in general. Figure 2.3 shows a

state vector and its corresponding velocity for some dynamic system. The angle 9

can be computed in the following way.


XTk
cos(0) = xl (2.5)

For k to have an inner radial component then

S< 0 < (2.6)
22

This is true when

xT < 0 (2.7)

where

k = Ax (2.8)

So if xTAx is negative everywhere for all time, then the system is stable. Since the

system is linear, it is sufficient to check the sign of xTAx on the unit circle. If all

velocities on the unit spheroid point inside that boundary, then the same is true for





















10 S \


4 ) 4 ^ '\

0 F. -' >^ \




-5o
-10 \ \ 0 10 is



151
-5 -10 .5 0 5 10 1I

XI



is 7 -A 1-
107





15 -1 0 5 1
.S ,' />




10 7 Ic j / /

-1 -10 *I 0 5 10 15
XI



10 -, ^ '- N


4 -4 -4 - ,- *-i *^





s-10


-s.. " 4-^ f 4-




-1 -10 -3 0 5 10 1S


4
15 -



















-1 -0 5 4 O 1
10

--^ -4 --
5^ -^ 'p x S.

0 '- 05 $ 0 i$

,'F -~ F? 4\ 4 ^


-10
4* 4- -4 - <\ <- -' /
-15 -
-15 .10 -5 0 S 10 15
xi

is

'F 4 ^

-0' *, 4 ^i \ \ .




"' \ ^ C


S 4L
-,.~I ^. >^-v

-10. \ I\ t 4
\ \ I\ \ I s

-15 -10 .5 0 5 10 15
XI

4
15 1


10

*^ ^ 4 4 4 ,






.5-
C*y 4- 4-^


-151 is/


XI


Figure 2.2. The trajectory and velocity field of (2.1)








15



2

1.5



0.5
0/


-0.5
-I

-1.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
xl

Figure 2.3. x and x

all spheroids centered on the origin. This leads to the following stability condition

which is given without proof.


Stability Condition 1: If xTic = xTAx < 0 Vx : XTX = 1,V t > 0 for

the system, 5 = A(t)x, then the origin is a global asymptotically stable

equilibrium point.


Figure 2.4 illustrates an example of a second order system which meets Stability

Condition 1. Note, that while all figures given so far represent second order linear

systems, the above stability condition applies to linear systems of any order.

The above stability condition can be made less conservative by expanding the

class of shapes to ellipsoids defined in the following way.


xTPx = 1 (2.9)

where P is a symmetric, positive definite matrix. Figure 2.5 shows an ellipse of the

form (2.9) defined for a second order system. The state vector, x, is drawn to some
































-1 -0.5 0 0.5 1


Figure 2.4. An example of a
ponents


second order system with negative radial velocity corn-


Figure 2.5. One velocity vector on the ellipse xTpx = 1


0










point on the ellipse, and its velocity, :, is also shown. The normal to the ellipse at x
is Px. To ensure stability, it is sufficient to require that the projection of k onto the
normal of the ellipse, Px, is negative everywhere. The resulting term, xTpk is the
normal velocity component on the ellipse. The main thrust of this paper is to be able

to control xTP* and to make this normal velocity component negative everywhere
for all time in order to achieve stability. As before, the linearity of the system allows

one to only check xTpk on the unit circle. This generalizes the previous stability
condition to

Stability Condition 2: For the system k = A(t)x, ifxTpk = xTPA(t)x <
0 V x : XTX = 1, V t > 0 for some constant positive definite matrix P,
then the origin is a global asymptotically stable equilibrium point.

The above condition applies to linear systems of any order.

Lyapunov [14] stated that if any positive definite function of the system states
V(x) is always decreasing, i.e. V(x) < 0, then the system is stable asymptoticallyy
stable for linear systems) [10]. For the linear case, let

V(x) = xTPx (2.10)

V(x) is positive definite if and only if P is a positive definite matrix. Taking the

derivative gives

V(x) = xTPx + xTp (2.11)

= xTATPx + xTPAx (2.12)

= xT[ATp + PA]x (2.13)

Lyapunov's stability condition becomes


xT[ATP+PA]x<0 Vx,Vt > 0


(2.14)










Since

XTpx = XTPAx = xT[ATP + PA]x (2.15)
2
Stability Condition 2 is equivalent to Lyapunov's linear stability condition.

It would be useful to compute the lower and upper bound of the rate of decay of
a given positive definite function,


V(x) = xTPx (2.16)

of the states of a time varying linear system by evaluating,

c := min [AminI PA(t) + AT(t)P] (2.17)
-- 1

c := max [Ama [PA(t) + AT(t)P]] (2.18)

So,

c < xT[PA(t) + AT(t)P]x < V x : xTx = 1, V t (2.19)
_2
Then c and Z are the lower and upper bounds respectively of the normal velocity
components, xTP*, on the unit spheroid. If c < 0, then the system is stable.

2.2 Adding Control to the System

The result in section 2.1 is an analytical tool only. A design procedure is needed
to control stability for the system

x= A(t)x + B(t)u (2.20)

where u is the control vector and x is the state vector. This section focuses on two

questions:

1) What effect does u have on the velocity field of (2.20)?
2) Can the normal velocity component xTP on the spheroid xTx = 1
be controlled?
























Figure 2.6. A velocity vector (a) without control and (b) with control.

To answer question 1), consider a second order single input system frozen at some

instant in time.

S= Ax + Bu (2.21)

A is a 2 x 2 matrix, and B is a two dimensional column vector, u is a scalar input

and can take on any real value. If u = 0 then the velocity field of the system is

k = Ax. Figure 2.6a shows the velocity, k = Ax, of some state vector, x, in the

system. The control vector, B, is also shown. When u : 0, 5 has the extra term,

Bu. The direction of Bu is constant, but its magnitude is directly proportional to u.

Figure 2.6b shows many possible values for by sweeping u through a wide range of

values in small increments. As this diagram shows, the arrow head of each velocity

vector can be placed on the line drawn parallel to B and intersecting the arrow head

of Ax. This demonstrates that velocities can be controlled along the space spanned

by the columns of B. The following theorem answers question 2).













B
0-5
0.5



-1
1.5 -
~2 -1.5 1 -0 .5 o o .5 1 '.5 2


Figure 2.7. A case in which the normal velocity component cannot be controlled.

Theorem 1
Consider the system, k = A(t)x+B(t)u, where dim(x) = n, dim(u) = m
and A(t) and B(t) have compatible dimensions. Let P be a constant
positive definite matrix. The normal velocity component, xTpk, can be
arbitrarily set to any value with the right choice of u at a given time t,
at any point x, on the spheroid, xTx = 1, unless x is contained in the set


S(t) := P-'span(B_(t)) n {x: xTx = 1}


(2.22)


where Bj(t) is a basis of column vectors orthogonal to span(B(t)). If
x E S(t) then


xTpk = xTPA(t)x


(2.23)


and this velocity component cannot be controlled at time t.


Theorem 1 gives the parts of the unit spheroid for which the velocity component
xTP* is uncontrollable and can therefore be used to determine if a linear time varying










system is stabilizable with an appropriate choice of constant positive definite P. If
there exists a constant positive definite P so that

max max xTPA(t)x < 0 (2.24)
t xES(t)
the system is stabilizable. The above condition is equivalent to requiring the following
matrix to be negative definite for all time.

IBT(t)[A(t)P-1 + P-1AT(t)]BL(t) (2.25)
21
This fact was established by Fields [15] and is true for the following reasons. Condi-
tion (2.24) is true if and only if

(P-l B (t),L)T(PA(t))(P-lB (t),t) <0VI' 0, Vt (2.26)
(P-1B (t)'U)T(P- (B (t)j')

Since the denominator of the above fraction is positive, with some simplification the
above statement is equivalent to

21T BT(t)[A(t)P- -1 <0 # 0, + t (2.27)
2 PA P- _L P^At)] y < 0 V IL k 0, V t (2.7

which is equivalent to

BT(t)[A(t)P-1 + P-lAT(t)]BI(t) <0 V t (2.28)

If a constant positive definite P exists satisfying condition (2.28), then the given time
varying system can be stabilized. This result is more general than a similar one given
by Shahruz [12], and the two are equivalent when P = I.
Before giving a proof of Theorem 1, we give an intuitive explanation of Theorem
1 for the second order case with a single input. Figure 2.7 shows the B vector for
such a system along with some ellipsoid xTpx = 1. If a state vector, x, is drawn
to some point on this ellipse whose normal is not perpendicular to B. then there is










always a control, u, which can place x inside the ellipsoid. This is true because the
velocity field can always be controlled in the B direction. This is not true, however,
for points on the ellipsoid whose normal is parallel to span(Bi.) as shown in figure
2.7. These points on the ellipsoid can be found by computing


P-'span(BI) (2.29)

where BI is a basis of column vectors orthogonal to span(B). We now give the proof
of Theorem 1.

Proof:
Case 1: x E S(t)

Substituting 2.20 into xTp*, we have

xTpk = xTPA(t)x + xTPB(t)u (2.30)

Since x E S then

x = P-'1(B.(t)) (2.31)

where l(B-(t)) is any linear combination of Bj1(t). So

xTPB(t)u = l(B_(t))TP-lPB(t)u (2.32)

= l(B(t))TB(t)u (2.33)

= 0 (2.34)

and therefore from 2.30 we have,


xTPk = xTPA(t)x


(2.35)










Case 2: x 0 S(t)
Suppose we want to set xTpik = c where c is some arbitrarily chosen
value. Then we want

c = xTPA(t)x + xTPB(t)u (2.36)

Since x S(t) then xTPB(t)u : 0 and there exists at least one u such
that 2.36 is true.

One could solve (2.36) for u to use as a control law, but this would not be practical
to implement because A and B can not easily be computed. Linear state feedback is
more desirable. The next section will develop a procedure for computing the set of
feedback gains that will implement a feedback control law keeping the normal velocity
component, xTPk, within some specified limit.

2.3 A Linear Feedback Set to Control xTP*

Recall from Section 2.1 that the stability of a time varying linear system

x = A(t)x (2.37)

could be analyzed by evaluating the bounds of xTPA(t)x on the spheroid

xTx = 1 (2.38)

Now we want to find the set of all linear state feedback gains for the control law,

u = -K(t)x (2.39)

for the system

x= A(t)x + B(t)u (2.40)

so that condition (2.19) holds for the closed loop system where c and are now
specified. This section will give conditions for the nonemptiness of such a set followed










by a discussion of its convexity. Finally this section will present an exhaustive search

method for finding the boundary of this set.

For the remainder of this chapter, A and B will be frozen at one instant in
time. A discussion follows on how to find a controller or set of controllers which

satisfies (2.19) at one given instant. If a given system is time varying then the results
which follow must be applied for all time with P being constant and positive definite.

These results apply to all systems which can be stabilized with respect to a Lyapunov
function using a constant P. The success of these methods depends on the existence
of a constant positive definite P so that

2BT(tt)[A(t)P-1 + P-AT(t)]B.(t) < 0 V t (2.41)

Since we require P to be zero, these results are conservative.

Let KA and TC be the set of all K which satisfy the left and right inequalities

respectively of the following expression.

c < xT[PA + ATP PBK KTBTPx V x : xTx = 1 (2.42)

where P is positive definite. The objective is to find


/C := k nC (2.43)

The following theorem gives conditions for the nonemptiness of _C and K?

Theorem 2 Given the 'nth order m-input system,

S= Ax+ Bu (2.44)

k and KT are nonempty if and only if there exists a positive definite P so
that the following conditions are true:

c < minxTPAx = min -xT[PA + ATP]x (2.45)
X ES xES 2










c> maxxT PAx = max 2xT[PA + ATp]x (2.46)
xES xeS 2

Comments:
Conditions (2.45) and (2.46) are respectively equivalent to

c < A,.m[I(vH/)-'B [AP-1 + P-1AT]BH(v')1- (2.47)

and

>_ Anax[(V'-H)-'BT[AP-1 + P-1AT]B(-vH)1] (2.48)

where

H = (P-1B)Tp-1B (2.49)

The matrix H is square, full rank, and has dimension n m. The above is true for
the following reasons. Equation (2.46) can be rewritten as


S> max 1xT[PA + ATp]x (2.50)
x=P-lBBL,,xTx=l 2
1 T T -1p1T1Rt.
= max -yB[APl +P-lATBI (2.51)
UTHL=I 2

Since H is positive definite then ~ exists, is square, and has an inverse. By making
the substitution, y = (v/H)-lz, (2.51) becomes

c > max zT(v/HB)-IBT[AP-1 + P-1A]B(v/H")-1z (2.52)
-zTz= 21
which is equivalent to

Am [vH)-1IT[AP-1 + P-1'A]B(vH)-] (2.53)


Equation (2.45) is equivalent to (2.47) for similar reasons. We now give a proof of
Theorem 2.










Proof:

We will show the nonemptiness of KX. KI can be shown to be nonempty in

a similar manner.

Case 1:

maxxxTPAx < (2.54)
xES
then

xTPAx < V x E S (2.55)



Let K = kBTP where k is a scalar value. Then


x = lxT[PA + ATP PBK KTBTP]x (2.56)

I XT[PA + ATP]x- k[xTPB][BTPx] (2.57)
2i 2

Since xTPB [BTPx]T then


xTPBBTPx > 0 V x S, V xTx = 1 (2.58)

If (2.55) is true and if k is made large enough, then the second term in

(2.57) will dominate V x S and x'p5c < V x : xTx = 1. Since

kBTP E E then E is nonempty.

Case 2:

maxxTPAx > (2.59)
xr=S
Then

XTPAx > for some x E S. (2.60)


By Theorem 1, xPAx can not be controlled V x E S; so TC is empty.










Remarks:
The nonemptiness of the intersection K: = 4c n is not guaranteed. If the designer
discovers that no intersection exists, then the upper and lower velocity bounds will
have to be adjusted.

Another useful property of CK: is convexity. This property is valuable in formulating
iterative search techniques to be described in the next section. K is convex if when

x and y are elements of K, then cax + (1 a)y is also an element of )C for 0
[16].

Theorem 3 Let K: be the set of all K so that

Amax [I[P(A BK) + (A BK)TP]] < (2.61)

Let IC be the set of all K so that

Amin[I[P(A BK) + (A BK)Tp]] > c (2.62)
2
Then /C, CK:, and : = K: fn kC are convex.

Proof:

If K: and K)C are convex, then the intersection K: is convex. Convexity of
kC will be proven here. In a similar manner the proof for the convexity of
KC can be written.

Let K1 and K2 be elements of T. We must show that aK1 + (1 ca)K2
is contained in KT when 0 < a < 1. We know that

IxT[PA+ATP-PBK1 -K KTBTP]x < V xTx = 1 (2.63)

and

lxT[PA + ATP- PBK2 KTBTP]x < V x = 1 (2.64)
2~ ~










Then

ax T[PA + ATP PBK1 KBTP]x < acE V x = 1 (2.65)

and

I(1-aQ)xT[PA+ATP-PBK2-KTBTp]x < (1-_a) VxTx 1 (2.66)
2

where 0 < a < 1. Adding (2.65) to (2.66) gives
xT[PA + ATP]x- _QXT[PBKl + KTBTP]x
(2.67)
-( a)xT[PBK2 + Kf2BTp]x < Z V xTx = 1
which becomes

1xT[PA+ATp-PB(aK+(1-a)K2)-(oaK+(1-a)K2)TBTP] < Vxx =1
2
(2.68)
So caK1 + (1 a)K2 E C for 0 < a c < 1. /C is convex.

Now that we have conditions on the nonemptiness and convexity of KC and k_ it
would be desirable to find the boundary of these sets. We know that if one of the
eigenvalues of the square matrix Q is c then

det[Q cI] = 0. (2.69)

This is helpful in understanding the following exhaustive search method for comput-
ing the boundary of T. A similar procedure exists for finding Q9_.

1) All eigenvalues of

I-[PA + ATP- PBK KTBTP] -I (2.70)
2
need to be less than or equal to zero. Fix all but one of the entries of K.
Let the i,j'th entry be the free entry and be represented by k. Let Koi,j










contain all the fixed entries of K and

Skl11


be defined as follows.
klj kin 1


kil 0 kin


kml k77j km,


Then


K = kQij + Koi,


1 j
0 0

0
0 1
0


0 0 0


(2.73)


2) Find all values of the free entry which make


det[PA + ATP 2cI PBK KTBTp] = 0


(2.74)


This problem reduces to finding the roots of a polynomial. After sub-

stituting equation (2.72) into the argument of the above determinant, it

becomes

k(-PBQi,j QJBTP) + (PA + ATP PBKoi,j KrjBTP) (2.75)

Let


Pl(i,j) := -(PBQi,j + QT.B TP)

Po(Koi,j) := PA + ATP PBKo,j KTojBTP


(2.76)

(2.77)


(2.71)


where


(2.72)


Koi,j :=


Qij "=










Then (2.74) becomes

det[kP1 + Po] = 0 (2.78)

Solve (2.78) for k.

3) Reject all complex roots. If all the roots are complex then skip the

next step.
4) Test the intervals between the real roots by checking to see if


Ama[PA + ATp 2-cI PBK KTBTp] < 0 (2.79)

The K's that bound the interval which satisfies (2.79) lie on the boundary

of T. Convexity of K implies that no more than two K's bound this

interval.

5) Repeat the process for all possible values for the fixed entries in K.

The result is 9K.

9KC can be computed in a similar way by reversing the inequalities in the above
procedure and by replacing max with Aiin. To find C the intersection of KC and TC

can be found in step 4).

The fixed entries in step 5) must be assigned to a finite number of grid points if

the above procedure is to be executed on a real system. The spacing of these grid

points must be smaller than the size of the feedback set. If the system has a high

order or a multiple number of inputs, the number of grid points will become too large,

and it will not be practical to implement this method. Understanding this procedure,

however, leads to the formulation of an iterative method that can be used on high

order, multiple input systems and will be presented in the following section. First,

an example is given applying the exhaustive search method to a second order, single

input system.











Example 1

Given the system


x= Ax + Bu


(2.80)


where

A '1 12 ]
0 -1 1
and
B = cos(75) ]
S sin(75o)
let

P=I

We want to find all feedback gains which satisfy the following constraint.


c

(2.81)


Before specifying c and c, we must check the normal velocity components
for x E S as Theorem 1 requires. Let


x=B = [ -sin75]

xAxo = -12.4075

1x TAxi =-12.40


(2.82)



(2.83)


From Theorem 2 we must have


c < -12.40 <


(2.84)


From this we choose


c = -14

-- -9


(2.85)

(2.86)


Then









In this example n = 2 and m = 1. It was decided to set i = j = 1 so
that

Q,j=Qi,i= [li 0] (2.87)

Koj = K 01,1 =[ 0 /i2 ] (2.88)

and
K = kQi,i + Koi,i = k K2 ] (2.89)

From equations (2.76) and (2.77) come

Pi_ = -[BQ,1 + Q,1B T] (2.90)
Pl = -[BQ1,, + QIBT] (2.91)

PoC = A + AT 2cI BKo0,1 KTBT (2.92)
Po- = A + AT 2I BKoi,1 KT' BT (2.93)

The roots of the following polynomials are computed in terms of k while
incrementing K2 through a wide range of values.

det[kP + Po] = 0 (2.94)

det[kP +Po0] = 0 (2.95)

Rejecting complex roots and checking the regions separated by the real
roots give

k!(K2), kc(K2) (2.96)
ku(K2), kj(g2) (2.97)

The intersection of these regions are found.

k(K2) = max[k,k-E] (2.98)
k(K2) = mink[,K] (2.99)












14

13 -

12

11 .[7,11]

10

9

81
4 5 6 7 8 9 10 11 12 13 14
ki

Figure 2.8. A plot of the boundary of IC which guarantees satisfaction of the design
constraints of Example 1.

Finally,


C C={[k(K2),K2] V K2 }u{[k(K2),K2] V K2} (2.100)


A plot of 9aC is shown in Figure 2.8. Let K = [7 11]. We can check to see

if K E KC by evaluating


A1 = Amin,(A+AT-BK KTBT) (2.101)

= -12.91 (2.102)


A2 = Am(A + AT BK KTBT) (2.103)

= -10.52 (2.104)


As the following shows,


c < A1 (2.105)

A2 < Z (2.106)










2.4 One Linear Feedback Matrix to Control xTP*

The search procedure given in section 2.3 becomes impractical for high order
systems with multiple inputs because the number of grid points for the fixed entries
of K becomes very large. This section discusses an iterative Lyapunov design method
which saves on computations and finds one element K E K, if it exists, where KC =
KC n KC. By applying this procedure at every instant in time to a time varying linear
system, one can find a control law stabilizing the system if a constant positive definite

P exists which satisfying the following condition.

2BT(t) [A(t)P-1 + P-1AT(t)]B(t) <0 V t (2.107)

By Theorem 2, this condition guarantees the existence of a Z < 0 such that C(t)
is nonempty for every instant in time. This procedure applies to all systems which
can be stabilized with respect to a Lyapunov function given by a constant positive
definite P. The following is a discussion of the iterative Lyapunov method followed
by the algorithm itself. An example is then given applying this procedure to one

operating point of a fifth order linearized model of the EMRAAT missile.
Figure 2.9 illustrates the iterative procedure. Kf is the feedback set which satisfies
the designer's predetermined constraints for some specified P. The constraints are

Cf < xT[PA + ATP PBK KTBTp]x < Zcf V x xTx = 1 (2.108)

Let Ki be defined as the set of all K which satisfies the following constraints.

ci < xT[PA + ATp PBK TBTp]x < cZ V x : xTx = 1 (2.109)

Given Ki, we would like to find c, and ci so that if Ki 0 kf then Ki E MKi. We also
require that Kf C Ki. If Ki E Kf, then we want K, = KC. The following definitions
for c, and ci meet these requirements. Let

i := max xT[P(A BKi) + (A BK)TP]x (2.110)
Ilx||=1 2













































Figure 2.9. A geometric view of the iterative Lyapunov design method.










and

Ai:= min Ix[P(A BKi) + (A BKi)Tp]x (2.111)
||x||=i 2 'J'
Then let

S= max[cf, \i (2.112)

and
ci = min[cf,A ] (2.113)

K0 is the initial guess in the search for Kf E kfA. Co and co are computed using
(2.112) and (2.113) so that Ko is a member of 9ko and /Cf C PC0. Then a new
feedback matrix, K1 is found which lies inside of K0, but not on the boundary. New
constraining values are found in the same way as before so that the boundary of the
next feedback set contains K1. K2 is then found so that it lies inside of the present
feedback set, but not on its boundary. This process is continued until Ki = Kf e ACf.

The success of finding Kf depends on the following conditions.

1. Given that Ki E a9)C we must be able to find Ki1 such that K++i G Ki.

2. We must show that 4i+1 < c, when -, > Zf and ci+1 > ci when ci < cf.

3. We must show that )Cf C kCi+l C Ki.

4. 'Cf must be nonempty.

We now address these four points.


1. Given Ki E aki, we need to find a second feedback matrix Ki2 C dKci where
Ki2 7 Ki. Then, due to convexity, 1 Ki + 1Ki2 is a member of kC{. Figure 2.10Oa shows
a second order example of a procedure for finding Ki+A The algorithm will be given
shortly. The horizontal and vertical axis are assigned to k1l and k12 respectively. The
region PCi is enclosed by aK&L and 9Ci. Kj is known. 'i2 is found by searching along










the line that passes through K, and is parallel to the k1l axis. Point c is found by
computing Ki + -Ki2. Since K;i is convex, then c E kC. This step is repeated again
by searching along the k12 axis. Using similar arguments, point f is also in KCi. KJ+i

is set equal to point f. For higher order systems, additional boundary points are
found and interior points computed by searching along directions which are parallel
to the axis of the coordinate system.

Under normal conditions, this procedure works well. However numerical problems
do occur. These problems will now be discussed along with a cure. Figure 2.10b shows
a second order example of when the above method fails. The boundaries are shaped

in such a way that searching along line 11 and 12 yields no new boundary points. A
solution to this problem is to relax one of the constraining values by setting, for ex-
ample, zi equal to + 6. 9Ad will then move so that point b can be found as shown in
figure 2.10c. Later, the relaxed constraining value can take on its original assignment.


2. We need to show that c, > Zi+l when c, > 5f. In this case


Ti= Ai > Zf (2.114)

Since Ki+j e KCU and Ki+j OQCi, then the following condition holds with strict
inequality.

xT[P(A BKi+1) + (A BKi+I)TP]x < a Vx : xTX = 1 (2.115)


Taking (2.110) for i + 1 then

IxT[P(A BKi+1) + (A BK<+)TP x V : xT = 1 (2.116)


and there exists an x which satisfies the above condition for equality. Therefore
Zi > AX+i and Zi > maxc.f,\i+l] = -i+n. Using similar arguments, it can be shown







38




a b


OKi





S< "/ /
*// Ii0o-/ /
K. a'ttily K ^2 / / /

/ K,
a' l1
/ a1 / h

k12
c


--I-

/
f' /t9

^i+i
S b

d ki


Figure 2.10. Step 5 of the iterative Lyapunov design method for (a) a second order
example, (b) how it sometimes fails, and (c) how this problem is corrected.











that c9 < cj+1 when ci < cf.


3. We now show that Cf C kCi+i C Ki. Since c, = max[f, Aj], then ci > c. We have

already seen that ci > >Ai+. So ci > max[f, AIi+] = ci+i. Similarly c ci+. Cki+l is

the set of all K so that

c I+ < xT[P(A BK) + (A BK)TP]x < c6+1 V x xTx = 1 (2.117)

Since T, > Ei+l and ci <_ 4i+1, then for every element of Kji+ the following holds.

c < xT[P(A BK) + (A BK)TP]x < Z, V x: xTx = 1 (2.118)

Therefore, .Ci+l C K;. Since c-i+ > cf and ci+l < cf, then using a similar argument

kCf C kA+1.


4. In using this design procedure, P is chosen so that the maximum uncontrol-

lable normal velocity component is negative. Then from Theorem 2, cf can be made

negative in an attempt to achieve stability. -f must be greater than the maximum

uncontrollable normal velocity component, and cf must be less than the minimum

uncontrollable normal velocity component. From Theorem 2 this will guarantee the

nonemptiness of Cf and kf. The nonemptiness of the intersection of these two sets,

however, is unknown. If Cf is nonempty, then, as i becomes large, Ki G Kfj. If Cf

is empty then ci and -. will converge to values which do not match the desired con-

straints and Ki will yield a closed loop system that meets the constraints given by cj
and ci. The designer will either have to accept this result or try again with a different

P or different constraining values or both. Since stability is desired, one approach

would be to keep P and Cf and lower Cf until Kf becomes large enough to intersect kf.


The outline of the iterative Lyapunov design method is as follows.










1) Choose P, Cjf, and Zf. This selection must obey Theorem 2. It should

be noted that Theorem 2 guarantees the nonemptiness of Kf and )Cf, but
not their intersections. If the system is time varying, then, in order to
use this algorithm, P must be found so that

2B T(t)[A(t)P-1 + P-lAT(t)]B(t) < 0 V t (2.119)

Otherwise this algorithm cannot guarantee stability.

2) Compute co and Co so that K0 = 0 E &aC0 and )Cf C Co. K0 will be

the initial guess in the search for Kf G ACf.

3) Let i = 0

4) Let i = Z+ 1

5) Find Ki so that KA E KC1I and KJi 0 M9Ai-.

6) Compute ci and c, so that Ki E 9C)i and KCf C /C, or so that Ki E kCf.

7) Repeat steps 4) 5) and 6) until one of two events occur.

1. ci = c and i = f .

2. (ci ci-1) and (i c,_i) become very small.

Remarks:
If event 1 occurs, then kCf is nonempty and Ki E Cf. If event 2 occurs, then kf is

empty and Ki yields a closed loop system that meets the constraints corresponding

to ci and ,i.
We now explain how to perform steps 2), 5), and 6).


Step 2)










The lower and upper bound of xTP* for the open loop system is respec-
tively,

Ao = Amird (PA + ATP)] (2.120)

Ao = Am,[ 2(PA + ATp)] (2.121)

Then,

o= min[ff,Ao] (2.122)

o = max[[cf ,Ao] (2.123)


Step 5)


The algorithm is now given.


K = K,
For j = 1 to m
For k = 1 to n
Koj,k = K kJ,kQj,k
Pi = -(PBQj,k + Q7kBTP)

Po = P(A BKOj,k) + (A BKoj,k)TP 2cl
Solve det[kP1 + Po0] = 0 in terms of k for c = ci, and Zi.
Reject all values which do not meet the constraints from
steps 2) and 6). The result is two intervals whose
lower bound is k_ and k2 and whose upper bound is
k and k2.
Find the intersection of these intervals by evaluating










k max[k, k2]

k = mn[Ti, -k2.
Find the midpoint of the interval by computing

kj,k = (k +).
Let K = ki,kQj,k + Koj,k
Next k

Next j
Ki+1 = K




Step 6)



Let

AiA = min[(PA + ATP- PBKi KTBTP)] (2.124)

SAmaxI(PA + ATp PBKi KTBTp)] (2.125)
2

So,

-* = min[-Cf,A] (2.126)

"i =max[-Cf, ] (2.127)

The following example illustrates the feasibility in applying this method to the

EMRAAT missile.

Example 2
We would like to apply the iterative Lyapunov design method to the EM-

RAAT missile. The missile was flown in a simulation through a trajectory











using another autopilot design. The model was linearized, and the fol-

lowing system was taken at 4.00 seconds into the flight.


x: = Ax + Bu


39 -.1
13 -.4
36 -1
.2 -6.
06 8;



B=


612 .0086
100 .1079
557 -2.078 -
097 .0173 -
L3.29 -.0261 -

0.0 -.1296
-.0149 0.0
-1150 -31.06
-4.494 -121.7
1.234 -.2802


(2.128)


.9997
0.0
.1763
.6403
.1597


0
.10-
-12:
-4.7-
-108


0.0
-.9997
.0442
.1611
-.5701

1.0
14
22
47
.6


(2.129)






(2.130)


x = [a,0,p,q,r]T


(2.131)


Step 1)

Let P = I. The minimum and maximum uncontrollable normal velocity

components were found by computing the limiting values in equations

(2.47) and (2.48). Since P = I, these terms simplify and are evaluated as

follows:


minlxT[PA+ATp]x= A'[mr[BT(A + AT)B.]
XS .7263
= -.7263


(2.132)

(2.133)


max xT [PA + ATp]x
xEs 2


= A [axc)I(A + A)B]

= -.3017


where


-.99;
.16
-64.:
-252
-.584


and


and


(2.134)

(2.135)










Theorem 2 guarantees the nonemptyness of T and K if f > -.3017 and

cf < -.7263. With this in mind, we let Zf = -.3 and cf = -6000. Since
the intersection of T and K is not guaranteed, cf was chosen to be very

negative to increase the probability of getting an answer.

Step 2)



A = A/,,[(A+ AT)] (2.136)
2
= -781.4 (2.137)

(2.138)


A = Ama,,[(A + AT)] (2.139)
2
= 778.9 (2.140)

(2.141)

So


co = min(cf,A) (2.142)

= -6000 (2.143)


Co = max(Zf,A) (2.144)

= 778.9 (2.145)

Steps 3) 7)
A program was written for MATLAB to carry out the iterations in steps

3) 7). The program would terminate if Kf was found or when

6 := I[i _i-1 < 6 = .0001 (2.146)











kf was found before 6 < .0001. The program ran 16 iterations. The final

result is


Kf =


.00046
2.057
-.00029


2.125
-.0012
-.7259


-.7974
2.339
.0263


.0805
-6.577
-.1346


.7535
.2591
-.7985


(2.147)


Now to check the result.


1
Af = A.n[(A+AT BKf KjBT)]
2f
= -818.0


Af = Amax[(A + AT BKf KfTBT)]

= -.3001

This meets the desired constraint.


af < Af < Af < Cf


(2.148)

(2.149)


(2.150)

(2.151)











CHAPTER 3
A TIME VARYING SECOND ORDER EXAMPLE

The following problem gives a case when pole placement succeeds in giving eigen-
values with negative real parts, but fails to stabilize the system. The Lyapunov
design method is then employed, and the resulting closed loop system is shown to be
asymptotically stable for all time.
We would like to find a feedback control law that stabilizes the system

x= A(t)x + B(t)u (3.1)

where

A (t) -1 + 1.Scos(t)[cos(t) + cos(t + Ir/18)] 1 1.5sin(t)[cos(t) + cos(t + r/18)] 1
A -1 1.5cos(t)[sin(t) + sin(t + 7r/18)] -1 + 1.5sin(t)[sin(t) + sin(t + 7r/18)]

(3.2)
and
B(t)= [ cos(t + r/18) (33)
sin(t + 7r/18) (
The eigenvalues of A(t) are -.4889 and 1.4661 for all time. The following control law
is proposed.

u = [ 1.5cos(t) -1.5sin(t) ]x (3.4)

The resulting closed loop system can be found in example 5.3,109 by Vidyasagar
[10] and also in Khalil [13]. The eigenvalues for the resulting closed loop system are
A = -.25 j.6614. Since the eigenvalues have negative real parts, one would expect
the closed loop system to be stable. However, Vidyasagar shows that the transition
matrix is
(t,, 0) = e5tcos) esin(t) (3.5)
1 _e.51sin^) e t cos(t)I












Trajectory of a pole placement time varying system


2501-


200

150

V 100


-50 -


-n00 -100 0 100 200 300 400 500
xl


Figure 3.1. The trajectory of the above closed loop system based
The initial conditions are x0 = [1 O]T. The system is unstable.


on pole placement.


If the initial conditions of the system are xo = [1 0]T, the resulting trajectory is

unstable as figure 3.1 shows.

We now turn to the Lyapunov design method. We choose P to be the identity

matrix. Before giving the design constraints, we need to check the the value of the

uncontrollable normal velocity components for all time. At t = 0


B = [.9848 .1736]T


(3.6)


Since P is the identity matrix, we are interested in the normal velocity component

which is on the part of the unit circle whose tangent is parallel to B. So we let


x = B = [.1736 .9848]T


The uncontrollable normal velocity component at t = 0 is


xTp = =xT[PA(O) + AT(0)P]x = -.9548
xT 2x


(3.7)


(3.8)










The uncontrollable velocity for this example is constant for all time. In selecting the

constraining values, we must have


c < -.9548 < Z (3.9)

The following assignments are made


c=-2 (3.10)

= -.5 (3.11)

Since the system is second order and has only one input it is possible to plot the set
of all feedback gains which satisfy the following.

C < lxT[P(A(t)-B(t)K(t))+(A(t)- B(t)K(t))TP]x < V x:xTx= 1, Vt (3.12)

This time varying feedback set is shown in figure 3.2. Since the time varying nature
of the system is periodic, and from inspection of figure 3.2, the following control law
is chosen.

u = 3[cos(t) sin(t)]x (3.13)

The feedback matrix in (3.13) is shown to be inside the moving feedback set in figure

3.2. The eigenvalues of the derivative of the resulting Lyapunov function is

A[A(t) B(t)K(t) + AT(t) KT(t)BT(t)] = -1.1193, -0.8579 Vt (3.14)

Therefore, the closed loop system is stable. Figure 3.3 shows a trajectory of the
Lyapunov based closed loop system where the initial condition is xO = [1 0]T.



























0


4 -2 2 4 -

____________ I) I p


4-
62 ---I -__________


6 -2 a 2 4
k11


-2


6 -4 -2 2 2 6
611
6 ____ d) t -3pir ___

4

2

0-


-6
-6 -O2 2 4
-8B 5 0 S ------ 4 ------


-26 7 2 4
k11 t.7p4


Figure 3.2. The time varying feedback set which satisfies the design constraints.








































-0.05


-0.1


-0.15


-0.2


-0.25


Trajectory of the Lyapunov based closed loop time varying system


: r /


-0.2 0 0.2 0.4 0.6 0.8 1
xl

Figure 3.3. A trajectory of the closed loop system using the Lyapunov design method.
xo = [1 O]J. The system is stable.













CHAPTER 4
DERIVATION OF A MODEL OF THE EMRAAT MISSILE

The next section derives a nonlinear model of the EMRAAT missile. A linearized

model is then generated in the following section and will be used for the design of the

autopilot. Aerodynamic and inertial cross coupling are assumed negligible in order

to reduce the order of the model. A linearized pitch model and a linearized roll-yaw

model results. All assumptions are clearly stated.

4.1 The Nonlinear Model

Before deriving the nonlinear model some variables and terms are defined. Figure

4.1 shows the missile body coordinate frame of the EMRAAT missile [2]. The three

axes, x, y, and z, are fixed to the missile as shown. The velocity of the missile is

represented by V. Angle of attack, a, is defined as the angle between x and the

projection of V onto the x-z plane. Sideslip, j, is the angle between x and the

projection of V onto the x-y plane. The velocity components, u, v, and w, are the

projections of V onto the x, y, and z axis respectively. The angle rates, p, q, and r,

are named roll rate, pitch rate, and yaw rate respectively and are defined as the rate

of rotation around the x, y, and z axis. Their directions obey the right hand rule.

These definitions are similar to those applied to aircraft as given by Etkin [17].

The following derivation of the nonlinear model is based on a similar model found
in Smith [2] and in Koenig [18]. Assumptions are made here to separate the system

into two lower order models: the pitch model and the roll-yaw model. We begin the

derivation by starting with the force equations. The forces acting on the missile are

thrust, gravity, and aerodynamic forces. The autopilot design in this paper applies

to the second phase of the flight when the engine is no longer burning so thrust is

















p



















r^ "- .- - ...q.


z




Figure 4.1. The missile body coordinate frame of the EMRAAT missile [2].










zero. Also, the weight of the missile is small in comparison to the aerodynamic force

so gravity is neglected. The aerodynamic force in the x direction is much smaller

than the aerodynamic force in the y and z directions, and therefore, will be ignored

for the rest of the paper. Newton's second law of motion implies the following:


Fy = m[i + ru pw] (4.1)

Fz = m[b + pv qu] (4.2)

Fy and Fz are the aerodynamic forces in the y and z directions. The quantities

inside the brackets are the total accelerations in the y and z directions. We know

that

V = [u2 + v2 + w2]1 (4.3)

Since v and w are much smaller than u then


V u (4.4)

Angle of attack and sideslip are given by


a = arctan (-) (4.5)

=3 arctan (-) (4.6)

If we assume that a and are small, then

w
a w- (4.7)
u
U

and
V
3-- (4.8)
U
Equations (4.1) and (4.2) can be rewritten as

Fy = mu + r p (4.9)











Fz = mu [+ p q]

which simplifies to

Fy = mV [- +r pal

Fz =mV [ + po3- q]

We assume that the forward velocity changes slowly so that


(4.13)


Then


w
U


V


and


So (4.11) and (4.12) become
Fy
mV
Fz
mV
The aerodynamic forces are given by


=O+r-pa (4.16)

-& +p p -q (4.17)


Fy = QS [Cyfi + Cyp + Cyjr + Cy6,bp + Cyl\6r]

Fz = -QS [CNQ + CNQ& + CNq + CNd6]q I


(4.18)

(4.19)


Q is the dynamic pressure and is defined as


Q= Ppv2


(4.20)


where p is the air density. S is the surface area of the wing. The aerodynamic

coefficients come from wind tunnel tests. They depend on mach number, and some

depend on angle of attack. The values of these coefficients have been put in tabular


(4.10)


(4.11)

(4.12)


(4.14)


(4.15)








form and are given in the first appendix along with other data pertaining to the
EMRAAT missile.
Substituting (4.18) and (4.19) into (4.16) and (4.17) and solving for & and /3 gives
QS
p= a r + -t [Cy + Cyp + Cyr + Cy,6/ + Cy,lr] (4.21)

and
S+ QSCN. q = QS [CN. + CNqq + CN,,6q] (4.22)

= [ + QSC [q p (CN.aa + CNq + CN,) (4.23)

Equations (4.21) and (4.23) are two nonlinear state equations. In order to separate
pitch dynamics from roll-yaw dynamics, /3 is assumed close to zero in (4.23). Thus,

[I [1 i QSCN]1 \ QSC \ + (I-' CAq) q +--2S N6q]

(4.24)
Equation (4.21) can be written as
S= ( y -" "-p(-1 -Cr) r-(--Cy6p) p-( P+ CY6,) 4,
\mV p/ V \ mV ) \mV p \mV
(4.25)
We now turn to the moment equations of the missile. Since thrust is zero and since
gravity does not contribute any moment to the missile, then the moments around
the x, y, and z axis are given by 1, m, and n respectively, the moments due to
aerodynamic pressure. They are given by

I = QSd [Ci3 + Cip + Cr + C,6p + Clfr\] (4.26)

m = QSd [Cma + Cm,& + Cmqq + C,6q 6,q] (4.27)
and


S= QSd [Co, + Cn,p + C1r + Cn,,p6 + Cn6r6r


(4.28)










where d is the missile diameter. Again, the aerodynamic coefficients come from
wind tunnel testing and are tabulated in the first appendix. Euler's three moment
equations can be written as follows.


m =J c + H(p,q,r) (4.29)
n r
where
Ixx -IXY -Ixz
J= -IXY IYY -IYz (4.30)
-Ixz -IYz Izz
and
-(Iyy Izz)qr + Iyz(r2 q2) Ixypq + Ixyrp
H = -(Izz- Ixx)rp + Ixz(p2 r2) Ixyqr + Iyzpq (4.31)
-(Ixx Iyy)pq + Ixy(q2 p2) Iyzrp + Ixzqr
We will assume that the inertial cross products, Ixy, Ixz, and Iyz are small. Also
due to symmetry of the airframe, Iyy and Izz are assumed to be equal. Inertial data
for the EMRAAT missile is given in the second appendix. Solving (4.29) for p, 4,
and r gives
S= --I (4.32)
'xx
= -- [m + (Izz Ixx)rp] (4.33)
IYY
and
S= -- [n + (Ixx Iyy) pq] (4.34)
'zz
Substituting (4.26), (4.27), and (4.28) into (4.32), (4.33), and (4.34) gives

QSd [ + Cp + Cr + C1 6+ C,] (4.35)

iQSd Izzxx
q QSd [Cmoa + Ca& + Cmq + Cm6qSq] + Izz Ixxrp (4.36)
IYY i 1* YY
and
____ '~ xx Iyy
S= [C/3 + Cp + Cnr + C 6,6p + C, 5] + 1- pq (4.37)
Izz Izz









If gyroscope effects are considered small, then (4.24), (4.25), (4.35), (4.36), and (4.37)
can be written as two separate systems: the pitch dynamics, and the roll-yaw dy-
namics. The pitch dynamics are as follows.

&=[ + QS\ QNa l[(-2-. Qc e + (I ( -'CN,) q+ (+Q CN6) 6q]
L rV J \mV } \ mV *} \mV *
(4.38)
QSd [C.a + Cm& + Cqq + C,,,q (4.39)

The roll-yaw dynamic equations are

QS (c + 2 Cy ))p (I + Q- (,) 2Cr+(Q2S C, ('S ) c r)
\mV */ \ mV *) \ mV mV p) \mV
(4.40)
sd[C,3 + C1pp + C r + C6pSP + C \r] (4.41)
Ix-x

= Q-d- [C,3 + Cnp + Cn1r + C,6p6P + Cn,65,] (4.42)

The states of the pitch model are a and q with 6q as its input. For the roll-yaw
model, the states are /3, p, and r, and the inputs are 6p and 64.

4.2 The Linear Model
The previous section gave a model of the pitch dynamics and the roll-yaw dynam-
ics in the following form.
[ = fq(a, q, q) (4.43)

[/ ]T fpT(),p,r,6p,6r) (4.44)

We would like to have two linearized models for use with the proposed design method.
The result will be two systems in the following form.

x = Ax + Bu (4.45)









where x contains the states and u contains the inputs. A and B are matrices which
are functions of several time varying flight parameters and are computed as follows.

A = df (4.46)
dx (X,W)nominal

dB f (4.47)
B=du
d(X,W ).ominal
where w contains additional flight parameters. Note that x and u are now pertur-
bations from the point around which the linearization is taken.
We linearize the pitch model first. From inspection of (4.38) we see that

=1 ( QSCNQ 1 (-QS CN
aqll = [1 +-'---- mI --- a
M mV M mV
aq12 +( QSCN (1- QS CN,)

bq= (I + QSCN) (- CN q)
M n-V M \mV q
To linearize (4.39), we must first substitute (4.38) in for &. Then we differentiate as
in (4.46) and (4.47).
QSd C + & QSCN (-QS \
aq22 = 7y [cm + Cm ( + ~mV~) ()]
h[U [( QSCNj1 (--QSCN ]

bq21 = QSd Cm6q + Cm (11 + QSCN) (QS0)]
Ivy Irn V M
The resulting linearized pitch model is
F 1 a=1, aq12 a I + bgil I 86 (4.48)
aq2l aq22 q bqg21
The same procedure is applied to the roll-yaw model. Assuming that a is constant
in (4.40), and from inspection of equations (4.40) to (4.42) the following results.
QQSCy QSCy, QSCyr
aprIl MV apr12 = a + ,-V apr13 = -1 + m"V







59


bpr1i QSCYP bpr12 QSCy6,
mV 'mV '
QSdCi QSdCi, QSdCh,
apZ1 ~ xI apr22 ~ X ~,apr23 -
QSdCj QSdCj ,
bpr21 XX Q pr22 Ixx
QSdC,, QSdCn, QSdCnr
apr31 = Z apr32 IZZ -I pr3 ~ IZZ '

_Q~rdC'6^ QSdCn6T
bpr3l QSdC,6 bpr32 QSdC

The linearized model is given by

apri apr12 apr13 + bpr11 bpr2 1 (p
p = apr21 apr22 apr23 p + bp21 bpr22 (449
1 apr31 apr32 apr33 i "r bpr31 bpr32












CHAPTER 5
THE DEPENDENCE OF GAINS ON FLIGHT PARAMETERS

This chapter describes the procedure that was used to determine which flight
parameters would be scheduled against the gains. Figure 5.1 gives a block diagram
of the system used for this procedure. The feedback gains are computed using the

iterative Lyapunov design method described in Chapter 2. The gains depend on
the linear model which in turn depends on seven flight parameters. Six of these
parameters are held constant while the seventh one changes. The resulting gains are
checked to see if they depend on this changing parameter. The process is repeated
for each flight parameter. It is found that the gains depend on angle of attack, mach
number, and dynamic pressure. This information will be used in the gain scheduling

process. The next section will discuss the atmospheric tables used to generate p and
V from M and Q. The flight parameter generator will then be presented. The third
section describes the initializer. The details of the iterative Lyapunov design method
are given in the fourth section. Finally, the results of the comparison between the
gains and flight parameters will be given in the fifth section.

5.1 Generating p and V from M and Q

By inspection, the linear models from Chapter 4 clearly depend on a, 0, p, q, r.
p, V, and the aerodynamic coefficients. The coefficients, however, can be eliminated
from the list because they depend on mach number and angle of attack. It is desirable

to replace p and V with M and Q since the later two can be easily measured on the
missile. We know the following.

2
Q = ^P2(5.1)

M = (5.2)
Vs0
60



















































Dependent
Flight
Parameters


Figure 5.1. A block diagram of the system used to determine the dependence of
gains on flight parameters.









where Vsos is the speed of sound. Both p and Vos are functions of altitude.

p = fp(h)

V8s0 = f(h)


where h is altitude in feet above sea level. Here,
atmospheric tables and are implemented by linear
and substituting the result into (5.1) gives


Q =


We generate a third table in the following way.


fp and fs are functions based on
interpolation. Solving (5.2) for V


(5.5)


f3(h) := f(h)f(h) = pV


= 2


(5.6)


(5.7)


The function, f3, is a one-to-one function so that its inverse can be found by reading
the table backwards. With this in mind we can solve (5.7) for h.


h f (2Q)
h= f31M


(5.8)


From (5.2) and (5.4)


V = MV,,, = Mf,(h)


(5.9)


Substituting (5.8) into (5.9) to eliminate h gives


S= Mf, (f- (2Q))

f (2Q))
-fP ( M2^


(5.10)


Also from (5.3)


(5.11)


(5.3)

(5.4)










Equations (5.10) and (5.11) are used to eliminate p and V in the linear model. As
a result, the linearized model can now be generated from the following seven flight

parameters.

[M Q ca 3 p q r] (5.12)

5.2 The Flight Parameter Generator

A series of flight conditions are made and used to generate many linear models.
Feedback gains are generated for each condition. The first of the series is called
the nominal flight condition. The values of the parameters for the nominal flight

condition are

M = 2, Q = 1250 psf, c = 8,/3 = 0,

p = 00/s,q = 100/s,r = 00/s.

Next, one of the parameters is allowed to vary while the other six are held constant.

This process is repeated six times so that each parameter can be tested. Table 5.1

shows the starting point, and the minimum and maximum values of each changing

parameter. Parameters with nonzero starting points begin at the starting point,

increment to the maximum value, return to the starting point, and then decrement

to the minimum value. Due to symmetry, the remaining variables have starting points

at zero and increment to their maximum values. M and Q sweep through a narrow

range because of restrictions of the atmospheric tables.

5.3 Initializing the Iterative Lyapunov Design Method

The iterative Lyapunov design method requires an initial guess, Kqo and Kpro and
two positive definite matrices, Pq and Pp,,. The initializer supplies these values and

will be discussed in this section.







64


Table 5.1. The initial point and range of changing flight parameters.
Initial point Minimum value Maximum value
M 2 1 2.6
Q 1250 psf 700 psf 5000 psf
a 8 _-80 200
/3 00 00 100
p 00/s 00/s 500/S
q 10/s -101/s 200/s
r 00/s 0/s 200/s

The initial feedback gains are found by using a pole placement algorithm. At the
nominal flight condition, the linear models are given by

A [ -1.1345 0.9996 [Bq -0.14631 (5.13)
S -261.4732 0.6209 ,Bq -123.1091 '

Apr -25. [ 24 1 .1350.9996 -017.5143 5. (5.13)
-0.459 0.140 -.100 0.018 0.117
Ap, = -2255.5 -2.41 .066 ,Bp = -1173.5 -1335.6 (5.14)
73.0 -0.181 -.648 L 2.01 -114.4
The desired eigenvalues for the closed loop pitch dynamics have been chosen to be
-40jlO0. For the closed loop roll-yaw model the desired pole locations are -20j5,
-80. The resulting feedback gains are

S\ 11.1633 -.6233 1 = -1.4166 -.0758 .3758 (5.15)
L 2.9337 .0086 -.3300 -



For both closed loop systems, P must be found so that xTpx is a Lyapunov
function. The following problem is stated.

Given a stable linear system : = Ax, find a positive definite function,
V = xTpx so that V1 = xT(PA + ATP)x is a negative definite function.

Let A be put into Jordan canonical form.

A = SJS-1 (5.16)










where S is an invertible matrix. The diagonals of J are the real parts of the eigen-

values of A, and imaginary parts of the eigenvalues lie in skew symmetric locations

off of the diagonals. For example, if the eigenvalues of A are a jb, c, then

'a -b 0
J= b a 0 (5.17)
0 0 c
Let

J = Jdiag + Jskew (5.18)

where Jdiag is the symmetric part of J and Jakew is the skew symmetric part of J. In

the above example

a 0 0 0 -b 0
Jdiag 0 a 0 Jskew = b 0 0 (5.19)
0 0 c 0 0
Making the following transformation on the system k = Ax, let


x = Sz (5.20)

Then

Si = ASz (5.21)

and

z = S-'ASz = Jz. (5.22)

Let Pz = -Jdiag. We wish to check the velocities of the system z = Jz on the ellipsoid

zTPzz = -zTJd"gz = 1 (5.23)

The normal to the ellipsoid at z is -Jdiagz. The projection of z = Jz onto the normal

is JdiagZ. Here the normal velocity component has the same magnitude but opposite

direction to the normal of the ellipsoid. If the system has no complex eigenvalues

then the velocities are orthogonal to the ellipsoid everywhere. For this reason, the










choice of P- = -Jdiag is the best choice for a positive definite function for the system

S= Jz.


V(z) = -zT JdiagZ


(5.24)


Making the following transformation into the x coordinate system gives


z = S-ix


(5.25)


which implies


V(x) = -xT[S-1]TJdiagS-X


Our choice of P is

P = -S-1]TJdiagS-.

For the nominal flight condition, the Jordan canonical form of Ag -
Apr- BprKpr is found and from (5.27)


S261772 60 1 ] 6792 10.7
Pq 609 14.9 Pp= 10.7 4.02
P 60 14[9 -316.4 -.509

The eigenvalues of (PqAq + ATpq PqBqKq KTBTpq) at


-316.4
-.509 (5.28)
15.7 j

the nominal point are


-1.0715 x 106, -40.002


(5.29)


Likewise, for the closed loop roll-yaw system the eigenvalues are


-1.3614 x 105, -320.00,-20.003


(5.30)


5.4 The Iterative LvaDunov Design Method


The iterative Lyapunov design method generates feedback gains so that xT PqX
and xTPp'x are Lyapunov functions for each closed loop system. The algorithm
requires the initial guesses Kgo and KprO, for the first point, and the positive definite
matrices Pq and Ppr. As a given flight condition changes, the feedback gains from


(5.26)


(5.27)

BqKq and










the previous point become the initial guess for the present point. The result of

this algorithm is a series of gains; one set for each flight condition generated. The

next portion of this section discusses some modifications made to the algorithm from

section 2.4. The material which follows describes how the design constraints are

selected.

The design method of section 2.4 finds a K so that the eigenvalues of (PA+ATp)

fall between cf and Cf where A is the closed loop system. The design algorithm

used in Figure 5.1 has been modified so that the resulting K satisfies the following

requirements.

cl < Ai ([P(A-BK) + (A-BK)TP]] < Z,
c2 < A2 [P(A- BK) + (A -BK)TP]] < 2 (5.31)
... ... ... ... ,...

c, < A,[I[P(A BK)+ (A BK)TP]] < n
where A1 is the smallest eigenvalue, A2 is the second smallest eigenvalue, and so on.

The values of the c's are supplied by the designer. This modification has been made

with the belief if more design constraints are made, then the performance will vary

less with changes in the flight conditions. The modified algorithm is as follows.

1) Choose P, e1f, cif, ... c,, and Z,,f. The selection of clf and cnf must

obey Theorem 2. However, Theorem 2 does not guarantee the nonempti-

ness of kC.

2) Compute cl0Cio, ... ,cn,,no so that the initial guess, K0 E Co0 and

JCf C /Co.

3) Let i = 0

4) Let i = i + 1

5) Compute Ki so that Ki E AC-1 and K, 9iC,-I.

6) Compute cli. ci, ... cni, Cni so that Ki E 9C, and KACf C Ci.







68

7) Repeat steps 4) 5) and 6) until one of two events occur.

1. cli= =Cl, Zli = Cll, ... cini = Cnf and Zni = Enf.
2. (Cli cli-1) *.. (ni- -,i-1) become very small.


Event 1. implies that the final answer has been found. Event 2. implies that /f
is empty and Ki satisfies the constraints corresponding to cli, ,i, ... cn, Z,. Steps

2), 5), and 6) are accomplished in the following way.


Steps 2) and 6)



Let
e = A1 [!(PA + ATp PBKi KTBTP)]

e = An [(PA + ATP- PBKi KTBTP) (5.32)
Then,
Cli = min[clf, el]
Eli = max[-cl, el]
... ...... (5.33)
fc = mzn[f _, e-]
-,i = max[nf, en]

Step 5)



Let K = Ki
For j = 1 to m
For k = 1 to n

P1 = -(PBQj,k + QBTP)
PO = P(A BKO,k) + (A BKo,k)TP 2cI

Koj,k = K kj,kQj,k










Solve det[kPi + P0] = 0 for in terms of k for c = cli, ... Ci.

Reject all values of k which do not meet the constraints from
steps 2) or 6). The result is 2n intervals whose lower bounds are

designated by k1, ... k2,n and whose upper bounds is named k, ..., k2n.

Find the intersection of these intervals by evaluating

k = max[k1, ...,1k2,]

k = mn [kl,...,k2n]
Find the midpoint of the interval by computing

kj,k = (k+k).
Let K = kj,kQj,k + Koj,k

Next k

Next j

KA+i = K



Remarks:

Let the feedback sets ),C1, k 2i &2, 2... ,, and C n be defined respectively as the

set of all K so the constraints (5.31) are met. Theorem 2 provides conditions for the

nonemptiness of K, and kAn. But conditions for the nonemptiness of the remaining

sets are still unknown. Also, )CI, ... ,_ are not convex in general. These are the

limitations of using the modified design algorithm.

We now turn to selecting constraining values for the eigenvalues of [P(A BK)+

(A BK)TP]. It is necessary to evaluate the uncontrollable normal velocity compo-

nents for each flight condition that will be generated in Figure 5.1. Tables 5.2 and 5.3

show for both models the minimum and maximum uncontrollable normal velocities

for each changing flight parameter.











Table 5.2. Uncontrollable normal velocity components for the pitch model
changing parameter minimum uncontrollable velocity maximum uncontrollable velocity
M -43.1456 -41.7253
Q -45.0144 -42.2289....
........ -43.3856 -42.2362....
__________ -42.6047 -42.6047
p -42.6047 -42.6045
q ..... -42.6047 -42.6047
r -42.6047 -42.6044


Table 5.3. Uncontrollable normal velocity components for the roll-yaw model
changing parameter minimumuncontrollable velocity maximum uncontrollable velocity
M -21.3937 -20.9751
Q -22.3016 -21.1574
_a_______ -21.7136 -21.1140
________ -21.3180 -21.3180
p -21.3181 -21.3180
q -21.3184 -21.3180
r -21.3180 "-21.3180..


For the pitch model, the uncontrollable normal velocity components range from

-45.0144 to -41.7253. Theorem 2 requires that


-41.7253 < q2
!ql < -45.0144


(5.34)


In addition, from (5.29), we want


cq1 <
_q2 <


-1.0715 x 106
-40.002


< Cql
< Cq2


(5.35)


From this the constraining values for the pitch model have been chosen to be


!21 = -1.2 X 106, qC- = -1 x 106,
Sq2 = -45, zq2 = -35


(5.36)


Likewise for the roll-yaw model, when looking at Table 5.3, Theorem 2 requires that


-20.9751
gprl


Cpr3
-22.3016


(5.37)










The eigenvalues in (5.30) suggest the following.

Cp l < -1.3614 x 105 < cprl
cpr2 < -320.00 < Cpr2 (5.38)
Cp,3 < -20.003 < Cp,3
The following constraining values have been chosen for the roll-yaw model.

cpi = -150000, Cprl = -110000
Cp,2 = -340, Zp, = -300 (5.39)
Cpr3 = -22, Cpr3 = -17
5.5 Formulation of a State Tracker

The autopilot of the EMRAAT missile will be a state tracker. That is, we want to
be able to change the location of the equilibrium point in order to control the values

of some of the states. The following shows how this will be accomplished.

Given the linear system

x= Ax + Bu (5.40)

y =Cx, (5.41)

we would like to find a control law


u = -Kx + KrefV (5.42)

so that y, the output, tracks v ,the reference input, asymptotically. We require y = v

when 5 = 0. When 5 = 0, then


0 = Ax BKx + BKfv. (5.43)

Since K is chosen so that the system is stable, then A BK is invertible and


x = -(A BK)-1BKrefV (5.44)

Also,


v = y = Cx = -C(A BK)-BKref V


(5.45)









Because (5.45) is true for any v, then

I = -C(A BK)-BKref, (5.46)

Controllability of the system implies that C(A-BK)-'B is invertible. Solving (5.46)
for Kref gives

Kref = -[C(A BK)-'B]-1 (5.47)

The EMRAAT missile has three inputs and therefore only three states can be
tracked. Controlling a in the pitch model and /3 and p in the roll yaw model is
desirable. For the pitch model y = a, implying that

Cq =[1 0] (5.48)

and, therefore,

Kef, = -[Cq(Aq BqKq)-'Bq]- (5.49)

For the roll-yaw model

Pi 0 o1 J (5.50)
Y = p 0 1 0
and, thus,

Krefr = -[Cpr(Apr BpKp1Bpr]'1 (5.51)

Krfq and Krefpr are computed for each flight condition and then compared along
with Kq and Kpr to the flight conditions.

5.6 Comparing Gains with Flight Parameters

A series of gains have been generated as a function of different flight conditions.
Each flight parameter has been swept through a range of points while the remaining
six have been held constant. In order to compare different gains on the same input
it has been decided to use the products of gains and their corresponding terms at
typical values. For example a typical value of a is 8. So we set ao equal to 8











Table 5.4. Extreme values and range of the pitch channel control terms as individual
flight parameters vary. Gains depend mostly on M, Q, and a._ __
_____M_. M] Q ... p q r
control min min min min min min min
term max max max max max max max
diff diff diff diff diff diff diff
kqlIao -2.15 -3.1 -2.1 -1.60 -1.61 -1.61 -1.60
-.72 -.16 -1.4 -1.59 -1.59 -1.59 -1.59
1.43 2.95 .69 .017 .017 .017 .017
kqi2qo -.1509 -.21 -.16 -.118 -.118 -.118 -.118
-.0529 -.031 -.10 -.117 -.117 -.117 -.117
.099 .179 .054 .0005 .0005 .0005 .0005
kref qllcO -2.53 -3.5 -2.7 -1.98 -1.98 -1.98 -1.98
-.916 -.53 -1.7 -1.97 -1.96 -1.96 -1.96
S 1.61 3.0 .94 .017 .017 .017 .017

and look at kqllao. We set qo and aco equal to 10/s and 8 respectively so that

we can look at kql2qo and kref qllao. The sum of these three terms are fed into to

the elevator. For the terms which are fed into the remaining inputs, the following

assignments are made.
/o = = 2
po = p = 100/s (5.52)
ro = 25/s
Figures 5.2a-g show kqi ao plotted against all seven flight parameters. These seven
figures show that kqnao changes with M, Q, and a but remains nearly constant when

/3, p, q, and r change. Similar figures exist for the remaining eleven gains and are

summarized in Tables 5.4 and 5.5. The minimum and maximum values of each term

is listed for each changing flight parameter along with the difference between the

minimum and maximum values. Angles are expressed in radians. From this table it
was determined that all gains will be scheduled against M, Q, and a.






















































































-1.59



-I,








.1.61 o v ) ) --------
0 20 40 6o a6 106 120 140 160 10 206



Figure 5.2. kqliao vs. a)M ; b)Q; c)a; d)3; e) p; f) q; g) r














Table 5.5. Extreme values and range of the roll-yaw channel control terms as indi-
vidual flight parameters vary. Gains depend mostly on M, Q, and a.
SM Q 1 p q -r
control min min min min minm imm mmin
term max max max max max max max
diff diff diff diff diff diff diff
kpli/3o -.050 -.11 -.22 -.050 -.050 .050 .050
.022 .045 -.013 -.038 -.036 -.036 -.036
.072 .157 .21 .011 .014 .014 .013
kpTrl2po -.17 -.14 -.29 -.13 -.13 -.13 -.13
-.011 -.061 -.13 -.13 -.13 -.13 -.13
.16 .20 .15 .001 .001 .001 .001
kpr13ro .074 .047 .15 .16 .16 .16 .16
.18 .28 .25 .16 .17 .17 .17
.10 .23 .099 .0006 .001 .001 .0009
kpr2ito .026 .0062 .065 .093 .092 .092 .092
.11 .167 .15 .10 .101 .101 .101
.080 .16 .084 .008 .009 .009 .009
kpr22po -.054 -.18 .0015 .014 .014 .014 .014
.015 .015 .152 .014 .015 .015 .015
.068 .20 .15 .0007 .0008 .0005 .0007
kpr23ro -.151 -.25 -.24 -.15 -.15 -.15 -.15
-.067 -.034 -.11 -.14 -.14 -.14 .14
.084 .22 .12 .001 .002 .001 .002
kref pr11ii/3co -.14 -.21 -.16 -.14 -.14 -.14 -.15
-.064 -.050 -.13 -.13 -.13 -.13 -.13
.08 .16 .028 .011 .014 .013 .013
kref pr12PcO -.071 -.046 -.39 -.040 -.040 -.045 -.040
.029 .15 .056 -.039 -.040 .013 -.038
.10 .19 .45 .001 .001 .057 .002
kref pr21l3co .055 .030 .11 .12 .12 .12 .12
.126 .19 .17 .13 .13 .13 .13
.071 .16 .058 .008 .009 .010 .009
kref pr22PcO -.010 -.33 -.17 -.071 -.071 -.12 -.071
-.070 -.070 .24 -.070 -.070 -.066 -.070
.027 .26 .41 .001 .001 .051 .002












CHAPTER 6
COMPUTING LOOK-UP TABLES

In the last chapter we showed that the gains depend mostly on a, M, and Q.
This chapter describes the process of generating a look-up table of gains verses M,

Q, and a. First a grid of points is formulated. Design constraints are then formulated

based on the knowledge of uncontrollable velocity components of the linear models.

Finally, the gains are computed.

6.1 Determining a Grid of Points

A two dimensional grid of points for M and Q has been made and used for each
entry of a in the table for both the pitch channel and the roll-yaw channel. Q sweeps
through a wide range of values starting with 100 psf and ending at 15,000 psf. The

values of M were chosen so that each entry of M and Q lie in the atmospheric tables

used to compute p and V. As a result, the grid points are not rectangular. All

entries are restricted to values between M = .6 and M = 3.5 and must correspond

to altitudes between sea level and 50,000 ft.

Table 6.1 shows the values of a used in the look-up table for both channels. M
and Q sweep through all values of the grid previously mentioned for each value of a
in Table 6.1. The spacing between the grid points was determined in a trial and error

process. During the iterative procedure for computing feedback gains, the initial
guess for each point came from the result of an adjacent grid point. The closer the

spacing between adjacent points, the fewer iterations were needed to find the next

feedback gain. Numerical problems as described in section 2.4 and shown in figure
2.10 were encountered. When this happened some of the constraints were relaxed

so interior points in the feedback set could be found. Later these constraints were








77


Table 6.1. Values of a in the pitch controller look-up table.
a
Pitch Roll-Yaw
-3.0 -1.0
1.0 1.0
4.0 2.5
8.0 4.0
12.0 8.0
16.0 12.0
20.0 16.0
_____ 20.0

made more restrictive and returned to their original assignments. This became a

tedious process for some parts of the grid. When the number of iterations exceeded

100 it was decided to add more grid points so that the desired feedback gains could

be found in fewer iterations.

6.2 Formulation of the Design Constraints

Before using the iterative Lyapunov design algorithm, the uncontrollable normal

velocity components for the entire M-Q-a grid must be determined. This information

is needed to formulate the design constraints. This process is described by the block

diagram in Figure 6.1. Table 6.2 presents the minimum and maximum uncontrollable

normal velocity components of the pitch model throughout the M Q grid for each

value of a. Here P = Pq, the matrix computed during the initializing procedure

of the last chapter. Table 6.3 gives the same result for the roll-yaw model where

P = Ppr,. The extreme values of uncontrollable normal velocities for the pitch model

are -49.5985 and -37.3480. Also the eigenvalues of

I[Pq(Aq BqKq) + (Aq BqKq)TPq] (6.1)







P P
q'^ pr


Uncontrollable
Normal
Velocities


Formulation
of Design
Constraints


Figure 6.1. Formulation of the Design Constraints.


Table 6.2. Uncontrollable normal velocity components for the pitch model.
a(degrees) minimum uncontrollable velocity maximum uncontrollable velocity
-3 -47.1539 -37.4260
1 -47.0806 -37.3873
4 -47.1611 -37.3480
8 -47.9797 -37.4247
.. 12 -48.8796 -37.3906
16 -49.5915 -37.6144
20 -49.5985 -37.4355

Table 6.3. Uncontrollable normal velocity components for the roll-yaw model.
a(degrees) minimum uncontrollable velocity [maximum uncontrollable velocity
-1 -23.3399 -19.6307
1 -23.3774 -19.6409
2.5 -23.4445 -19.6883
4 -23.4971 -19.7450
8 -24.5852 -17.2696
12 -25.3163 -19.9600
16 -25.6480 -19.8625
20 -26.4797 -20.30717


Flight
Parameter
Generator










at the nominal flight condition from the previous chapter were found to be -1.0715 x

106 and -40.002. Theorem 2 requires that

-37.3480 < cq2 (6.2)
cqI < -49.5985
But we also want
cq1 < -1.0715 x 106 < ql (6.3)
-q2 < -40.002 < (6q.
since we desire the eigenvalues of (6.1) to be close to those at the nominal flight

condition. From this, the constraining values were chosen to be

CqI = -1.2 x 106, q1 = -106,

Cq2 =-45, cq2 =-35

Likewise, for the roll-yaw model, the uncontrollable normal velocity components range

from -26.4797 to -17.2696. From Theorem 2 we must have

-17.2696 < p3 (6.4)c
cpr1 < -26.4797 )
With the eigenvalues of

S[Pp,(Apr BprKpr) + (Ap. BprKp,)TPp.]
2
at the nominal flight condition being -1.3614 x 105, -320.00, and -20.003, the

following is desirable.

cp,1 < -1.3614 x 105 < CplI
cpr2 < -320.00 < Cp,2 (6.5)
cpr3 < -20.003 < Cp,3
With this in mind, the following selections were made.


.prl = -200000, CprIl =-90000,

pr2 = -400, Tp,2 = -250,


cpr3 = -25, Cp,3 = -15











6.3 Generating the Look-Up Table

With the design constraints set, feedback gains and feedforward gains are gen-

erated for each grid point. Figure 6.2 gives a block diagram of the system used to

accomplish this. For the pitch model, the initial guess comes from one of the gains

that was generated in the previous chapter when determining the dependence be-

tween gains and flight parameters. The operating point from which this initial guess

originates is

M = 1, Q = 1250psf, a = 8 (6.6)

and is one of the extreme values listed in Table 5.1. The result of the first point is

used to start adjacent points which, in turn, start new adjacent points until gains

have been computed for the entire grid. The look-up table for the roll-yaw model is

made in the same way.

Numerical problems were encountered in parts of the look-up table for the roll-

yaw model. They were similar to the problems that were predicted in step 5) of the

iterative design method given in Chapter 2. To overcome these difficulties, some of the

constraining values were relaxed for a number of iterations and were later returned to

their original assignments in the algorithm. Eventually, the desired feedback matrix

was found.

Figure 6.3 shows kq ii verses mach number and dynamic pressure when a = 8. As

this figure would indicate, the gains generated from the iterative Lyapunov method

are smooth with respect to the dependent flight parameters. This fact gives hope

that the gain scheduling scheme will be easy.

















































Gain
Schedules


Figure 6.2. Formulation of the look-up table.




















































10000


5000


0 0.5


Figure 6.3. kql, verses M and Q when a = 80













CHAPTER 7
GAIN SCHEDULING

The result of the previous chapter is thirteen tables of gains in terms of mach

number, dynamic pressure, and angle of attack. This chapter discusses the process

of curve fitting used to implement the look-up table. The results of a test of this

scheme will follow.

7.1 Curve Fitting

It was decided to use a combination of polynomial fitting and interpolation to

implement the autopilot. Third order polynomials were fit to the tables as a function

of mach number. However, low order polynomials could not achieve close fits as a

function of angle of attack or dynamic pressure, so linear interpolation was used for

these two variables. Figure 7.la shows kq11 as a function of M for various constant

values of Q at a = 8. An example of polynomial fitting of one of these curves is

shown in Figure 7.1b where Q = 150psf. A polynomial has been made for every

a-Q pair and the polynomial coefficients are interpolated as a function of these two

variables. The tables of the three pitch controller gains each have 7 entries for a

and 22 entries for Q. Since each polynomial has 4 coefficients the total number of

coefficients for each gain is 7 x 22 x 4 = 616. Similarly the tables for the roll-yaw

controller have 8 entries for a and the same number of entries for Q. Each of the ten

gains are then scheduled using 8 x 22 x 4 = 704 polynomial coefficients.

7.2 Testing the Fit

Figure 7.2 gives the system for testing the polynomial and interpolation routines.
Gains were generated from these routines at locations centered between the original

grid points. These new locations were found by taking the average of the coordinates

83












a 0
a
10
0_ .......

.20 .- -47
30,
i- .. 48.


-60
-70-50
.5 1 1.5 2 as 3 3.5 .6 0.65 0.7 0.75 0.8 0.A5 0.9 0.95 1
M M

Figure 7.1. A plot of the kqii verses M with (a) various constant values of Q and
Q = 8, and (b) a least squares third order polynomial fit for the plot where Q =
150psf and ac = 8.


of adjacent grid points. Linear models were also made at each test point where /, p,

q, and r were set to zero. The eigenvalues of


I[Pq(Aq- BqKq) + (Aq BqKq)TPq] (7.1)


and

j[Pp,(Ap,, Bp,Kp,) + (Ap. BpKp)TPP] (7.2)

were computed to see if they remained within the desired limits. Table 7.1 shows

the maximum eigenvalue of (7.2) for all of the test points at a = 1.75. These values

are plotted against their indices in Figure 7.3. Most of the eigenvalues of Table

7.1 lie within the desired limits of -25 and -15. The eigenvalue in seventeenth row,

second column, however, is -3.96, the worst value found out of all of the test points.

Although the deviation is high, this value is still negative indicating stability for the

closed loop system. For the pitch controller, the actual limiting values of A1 and A2

at all the test points are

-1.37 x 106 < Aq1 < -9.85 x 105 (73)
-47.7 < Aq2 < -38.1 ( )

















































Figure 7.2. A test of the curve fitting routines used to implement the autopilot.












Table 7.1. The maximum eigenvalues of (6.1) for all test points at a = 1.75.
Q/M indices 1 2 3 4 5 6 7[ 8
1 -20.79 -20.82 -20.84 -20.87 -20.88 -20.89 -20.90 -20.90
2 -20.68 -20.69 -20.70 -20.72 -20.74 -20.76 -20.76 -20.76
3 -20.84 -20.85 -20.84 -20.83 -20.81 -20.79 -20.78 -20.79
4 -20.59 -20.61 -20.49 -20.30 -20.18 -20.33 -20.34 -20.21
5 -20.30 -20.22 -19.93 -19.71 -19.96 -19.88 -20.01 -20.62
6 -20.41 -19.97 -19.77 -20.16 -19.88 -20.20 -20.79 -20.97
7 -19.10 -17.70 -17.86 -18.92 -20.03 -20.42 -19.75 -20.81
8 -18.62 -19.79 -20.67 -20.37 -16.16 -19.14 -20.21 -19.48
9 -19.90 -20.97 -19.33 -14.48 -19.58 -21.08 -20.93 -19.64
10 -20.23 -20.88 -11.67 -16.89 -19.96 -20.40 -19.72 -19.41
11 -20.57 -20.38 -12.40 -18.26 -20.60 -20.85 -20.28 -20.30
12 -20.77 -17.82 -17.66 -20.97 -21.85 -21.65 -21.14 -21.32
13 -19.83 -18.51 -21.91 -20.81 -20.87 -22.17 -22.13 -22.15
14 -20.45 -22.23 -20.57 -17.42 -20.84 -22.21 -22.42 -22.39
15 -22.42 -22.13 -18.83 -18.24 -20.65 -21.23 -21.49 -21.90
16 -21.82 -19.45 -15.17 -20.54 -21.09 -20.92 -20.94 -21.73
17 -21.85 -3.96 -20.63 -22.26 -22.35 -22.13 -22.03 -22.48
18 -18.88 -19.74 -21.48 -22.29 -22.58 -22.65 -22.68 -22.87
19 -20.85 -21.06 -21.71 -22.27 -22.69 -22.98 -23.19 -23.37
20 -21.65 -21.93 -22.29 -22.63 -22.93 -23.19 -23.39 -23.52
21 -20.62 -22.15 -22.52 -22.62 -22.61 -22.53 -22.46 -22.42


-10

M~ -20,
E
t-30,

-40-

-50,,
30 -


10 04
I ,,^ F 0 0 . .


Indices of M


Figure 7.3. A plot of the eigenvalues from Table 6.1.


lll in i c -








87

Likewise, for the roll-yaw controller, the limiting values at all the test points are

-2.03 x 105 < Aprl K -9.28 x 104
-2392.7 < Apr2 < -182.6 (7.4)
-26.27 < Apr3 < -3.96

Some of these values differ significantly from the desired constraints; however, since

these values are still negative, these deviations are acceptable and indicate that the

closed loop system will be stable. It should also be noted that most eigenvalues

remain well within their desired constraints as shown in Table 7.1.













CHAPTER 8
NONLINEAR SIMULATIONS

A nonlinear simulation has been used to test the proposed autopilot for the EM-

RAAT missile. First a section follows giving an overview of the nonlinear simulation.

A test module is then made to generate state commands in order to evaluate the

autopilot's tracking ability. Finally, a series of flight scenarios are run to determine

the ability the missile has to intercept the target.

8.1 The Nonlinear Simulation

Figure 8.1 shows a block diagram of the simulation used to test the EMRAAT

missile. The program is written in FORTRAN. Initial conditions of the target and

missile are specified by the user. The simulation is then run and a trajectory of both

the target and missile results. All target and missile variables can be observed. The

target is programmed to fly in a straight line until the range between the target and

missile falls below 5,000 ft. The target then makes a 9 g turn to the right. The

simulation terminates when the closing velocity becomes positive.

The seeker measures the line of sight angles and the range rate of the target.

The simulation uses exact measurements of these values and does not assume any

noise. These values are sent to the guidance law which, in this case, implements

proportional navigation. A derivation of this guidance law can be found in Bryson

and Ho [19]. The outputs of the guidance law are two desired accelerations, a, and

a,. The BTT logic makes the conversion from the acceleration commands to the

three state commands ac, 0, and pc. Since the missile can achieve a much higher

acceleration with angle of attack than with sideslip, the BTT logic uses Pc to rotate

the desired accelerations into the pitch plane. If this roll maneuver is successful then


















Seeker
(RF)


Gain
Schedule




M Q .


K

K~f


n Guidance Law
R (Pro-Nay)
OQ
vc


ayc, azc

BTT Logic


aI O cl


PC


oi p q r


Sp
6q
Sr


Nonlinear
Missile
Dynamics




a Iaz


Exact Computation of Missile and Target Variables


Figure 8.1. A Block Diagram of the Nonlinear Missile Simulation.


Target
Position


U = -Kx + K,,fv


i[ l, i










ay, will become small and a, will become positive. ac and /30, are computed in an

attempt to match a^ and ayo respectively.

The autopilot implements the control law


u = -Kx + Krefv (8.1)

where x is a vector containing the actual states and v contains the state commands

from the bank-to-turn (BTT) logic. The states come from exact measurements in the

simulation. If this autopilot were to be implemented in an actual missile, the states

would be measured using an inertial platform. The gains K and Krf come from the

gain schedule implemented with a combination of polynomials and interpolation. In

this simulation there is no delay in the gain schedule and K and h,,f are produced

instantaneously. The output of the autopilot is the control surface angles 6p, 6q, and

86. Linear and angular accelerations are computed by the missile dynamics module
of the program. The simulation uses the output of the missile dynamics to compute

all of the flight variables including the position and velocity of the missile.

8.2 A Test of State Tracking

The model for the EMRAAT missile has five states and three inputs. The autopi-

lot is designed to track three state commands: ac, 3, and pc. Before running missile

target scenarios it was decided to test the autopilot's tracking ability. The BTT logic

was disconnected, and the following commands were applied to the reference inputs

of the autopilot.
100 for 0 s < t < .5 s
0 for .5 s < t < 2.75 s (8.2)
c= 10 for 2.75 s < t <3.75s (.
0 for 3.75 s < t
0' for 0 s /= 5 for 2 s 0 for 2.5s < t











Table 8.1. The rise times of each commanded state.
Altitude(ft) Mach t,(s) t, (s) 4.(s)
20,000 2.0 .100 .243 .0284
50,000 3.0 .108 .316 .0181

0/s for 0 s < t < 1 s
100/s for 1 s < t < 1.5 s
0/s for 1.5 s < t< 2.75 s .
Pc 100o/s for 2.75 s < t< 3.25 s (8.4)
-100/s for 3.25 s 0/s for 3.75 s < t
The experiment was performed once at an initial altitude of 20,000 ft with an initial

mach number of 2.0 and a second time at an initial altitude of 50,000 ft with an initial

mach number of 3.0. Figure 8.2 shows the results of the first test. In addition to the

commanded states, figure 8.2 also shows the remaining two states q and r. The rise

time as defined as the time needed to achieve 90% of the desired value was found

for each commanded state and is shown in table 8.1. It should be noted that mach

number and dynamic pressure do not change significantly during this test. This is a

test of accuracy in state tracking and is not a valid test of gain scheduling in terms

of M and Q. However, a changes very rapidly and does not appear to hinder the

performance.

Cross coupling is evident from this experiment. The step in the roll rate of figure

8.2c at t = 2 s is due to the sideslip command in figure 8.2b. A lesser degree of cross

coupling occurs during the roll command for 2.75 s < t < 3.75 s when a which is at

10 is rotated into sideslip. The same effects are present in the second experiment at

50,000 ft.

Figure 8.2d and 8.2e show the two states that have no commands. The spikes in

the pitch rate occur when ac changes value, because a pitching maneuver is required

to change the angle of attack. The two spikes in the yaw rate occur for the same





































b) C--) Co,,m,..c 0ta C-)


0 = a a a r a

a) p C -) ~,. Oo,,a~ 0 (-)
'Sc.


'ool


S.o


-0)


Figure 8.2. A test of state tracking at an initial altitude of 20,000 ft with an initial

mach number of 2.0. The states shown are (a) a, (b) 03, (c) p, (d) q, and (e) r.


S.---


. T


. -a ora.










reason when sideslip changes. The two steps in the yaw rate, however, are due to

cross coupling with the roll rate.

The rise times in table 8.1 are small and indicate that the missile will perform well

in flight scenarios. While cross coupling effects are noticeable, they are not believed

to be great enough to significantly hinder the performance of the missile; thus it was

decided to run the missile through a series of flight scenarios to see how well the

missile can intercept the target.

8.3 Simulation of Flight Scenarios.

A series of flight scenarios has been run to test the autopilot. Figure 8.3 shows

the trajectory of the missile and target for a flight scenario at an altitude of 20,000

ft. The miss distance is .64 ft. A hit is considered to be any miss distance under 10

ft. Figure 8.4 shows the commanded and the actual y and z accelerations. The devi-

ations from the commanded normal accelerations are mostly due to a simplification

used in implementing the BTT logic. Instead of using two aerodynamic coefficients,

proportionality constants were assumed. The errors are not large enough to prevent

the interception of the target. The state tracking of the reference inputs appears to

be working well as seen in figure 8.5 and indicates that the autopilot is performing

well.

The trajectories of a 40,000 foot altitude scenario is shown in figure 8.6. The miss

distance is 0.05 ft. A similar scenario at 10,000 ft is shown in figure 8.7. The miss

distance is 2.8 ft.

Figure 8.8 gives a case were a miss occurred. The scenario took place at 50,000

ft and the target was missed by 452 ft. Figure 8.9a and b shows that the missile was

unable to achieve the desired z acceleration and angle of attack after 1.3 seconds into

the flight. This is because the elevator reached its -40 limit as figure 8.9c indicates.























a) Top View of Missile (-) and Target (-) Trajectory


6000

5000

4000

3000


0 1000 2000 3000 400 5000 6008 7000 800O 9
Position in the x direction in ft.
xl04 b) Side View of Missile (-) and Target (-) Trajectory


-2.05 k


0 1000 2000 3000 4000 5000 6000 7000 8000
Position in thex direction in ft.


Figure 8.3. (a) Top view and (b) side view of missile and target
scenario which occurred at 20,000 ft. The miss distance is 0.64
numbers of the target and missile are 2.5 and .92 respectively.


trajectories of a
ft. Initial mach






















a) Commanded z Acceleration (-) and Acwal z Acceleration (-)


Seconds
b) Conmmanded y Acceleration (-) and Actal y Acceleration (-)


Seconds


Figure 8.4. Commanded and actual (a) z accelerations and (b) y accelerations for
the scenario in figure 8.3.