UFDC Home  UF Institutional Repository  UF Theses & Dissertations  Internet Archive   Help 
Material Information
Record Information

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 lookup 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 LOOKUP TABLES ...................... 6.1 Determining a Grid of Points ...................... 6.2 Formulation of :he Design Constraints ................... 6.3 Generating the LookUp 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 airtoair missile. The missile under study is the extended medium 'range airtoair technology (EMRAAT) missile, a theoretical banktoturn missile which is under study by the United States Air Force. Knowledge gained by this missile will be applied to future banktoturn 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 angleofattack. A combination of interpolation and polynomial fitting is used to create a lookup 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 banktoturn airtoair 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 airtoair missile is launched from a military aircraft with the intent of in tercepting an enemy aircraft or target. Banktoturn (BTT) missiles have wings giving them greater maneuverability than the conventional skidtoturn (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 airtoair 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 AA View AA 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. AB '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 llxllo 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(a1)tcos(t) etsin(t) 1 q^(t,0)= e(a1)sin(t) e* cos(t) (2.3) and the characteristic equation is A2+(2a)A+(2a)=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 05 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 pieshaped 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 *^ s10 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 05 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)P1 + P1AT(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 (Pl B (t),L)T(PA(t))(PlB (t),t) <0VI' 0, Vt (2.26) (P1B (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)P1 + PlAT(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))TPlPB(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)P1 + PAT(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 minput 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 [AP1 + P1AT]BH(v')1 (2.47) and >_ Anax[(V'H)'BT[AP1 + P1AT]B(vH)1] (2.48) where H = (P1B)Tp1B (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=PlBBL,,xTx=l 2 1 T T 1p1T1Rt. = max yB[APl +PlATBI (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[AP1 + P1A]B(v/H")1z (2.52) zTz= 21 which is equivalent to Am [vH)1IT[AP1 + P1'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+ATPPBK1 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(1aQ)xT[PA+ATPPBK2KTBTp]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+ATpPB(aK+(1a)K2)(oaK+(1a)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 2cI 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,kE] (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+ATBK 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)P1 + P1AT(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 ci+ > 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)P1 + PlAT(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 ci1) 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 _i1 < 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 rollyaw 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 xz plane. Sideslip, j, is the angle between x and the projection of V onto the xy 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 rollyaw 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+rpa (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 rollyaw 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= ( \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 rollyaw 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 rollyaw 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) Ixx = Qd [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 rollyaw 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 rollyaw 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 nV 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 rollyaw 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 onetoone 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 rollyaw 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 = SJS1 (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 = Six (5.25) which implies V(x) = xT[S1]TJdiagSX Our choice of P is P = S1]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 rollyaw 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(ABK) + (ABK)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 AC1 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 cli1) *.. (ni ,i1) 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 rollyaw 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 rollyaw 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 rollyaw 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(ABK)'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 rollyaw 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.2ag 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 rollyaw 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 LOOKUP 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 lookup 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 rollyaw 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 lookup 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 lookup table. a Pitch RollYaw 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 MQa 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 rollyaw 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 rollyaw 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 rollyaw 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 LookUp 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 lookup table for the rollyaw model is made in the same way. Numerical problems were encountered in parts of the lookup 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 lookup 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 lookup 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 aQ 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 rollyaw 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 7050 .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 t30, 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 rollyaw 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 (ProNay) 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 banktoturn (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 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 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. 