AlAA 90-2829 Optimal Multiple Vehicle Trajectories C. R. Hargraves and S. W. Paris ABSTRACT Three approaches for solving multiple vehicle trajectory optimization problems are Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 described. All three approaches are options that are available in a trajectory optimization computer code called OTIS (Optimal Trajectories by Implicit Simulation). Three examples of multiple vehicle problems solved by OTIS are described. INTRODUCTION Many contemporary trajectory optimization problems involve more than one vehicle. Examples are: c - Determine the minimum time intercept of a target flying a specified path. - Study orbit injections where one or more stages are to return to the launch site. Optimal control of an aircraft to evade a missile with proportional navigation g u i d a n.c. ~ W - The "standard" approach for solving these problems is to treat each one separately and to develop a specialized procedure for solving it. In this paper we are interested in demonstrating techniques that allow solution of a wide class of multiple vehicle problems as well as the more usual problems involving a single vehicle. The methods described below have been incorporated into a general purpose trajectory optimization code called OTIS. The OTIS code is based on a direct optimization approach using cubic splines, collocation and nonlinear programming. BACKGROUND A direct optimization method using collocation and nonlinear programming was described by the authors in Reference 1. These procedures were developed into a computer code called OTIS which is documented in Reference 2. Traiectory Modeling The concept of decomposing a trajectory into a linked sequence of arcs, called phases, is fundamental to providing a flexible simulation. The initial time and the final time for each phase are called events. An event may be thought of as an interruption of the current phase. Any time AlAA 90-2829 Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 the data describing a phase is to be changed, an event must occur. The desired simulation characteristics for the next phase must then be specified. The information that must be specified for each phase is: starting conditions; aerodynamic data; propulsion data; number of vehicles; jettison weight; environmental model; stopping condition; value of stopping condition; control angle functional form; numerical integration type; guidance type; integration coordinates (Cartesian, flight path or body); physical constraints; launch options; integration step size or accuracy; print step size; output options. The required number of phases is a function of the trajectory. In most cases phases are sequential. This is not required, e.g., two phases may start at the same time and end at different times. This provides the mechanism to model multiple vehicle problems A trajectory description by phase is illustrated in Figure 1. Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 phasei 1 I Phase 2 I 1 Phase 3 I Phase I Phase I I I I I I I I Alpha 0 I Thrust 0 Event E 1 E2 E3 E4 E5 E6 Figure 1 Trajectory Description by Phase Mathematical Problem Definition We want to determine the trajectories for one or more aerospace vehicles. The start and end times of each phase are defied as Ej and Ej+l, where j is equal to the phase number and NP is equal to the total number of phases. The parameters Ej, j = 1, NP + 1 are called events. Solutions are sought to the set of differential state constraints (the equations of motion, mass-flow equation, heat-flux etc.), where xl = (xi: i = 1, NS(j)) and NS(j) equals the number of states for phase j. In addition to the state variables, x, there are auxiliary variables, y, (dynamic pressure, heat flux etc.) defined as AlAA 90-2829 where yj = (xk: k = 1, NAG)) and NAG) is equal to the number of auxiliary parameters for phase j. Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 This allows path constraints of the form where NPCG) equals the number of path constraints for phase j. The control vector, u, (angle of attack, throttle factor, etc.) is defined as, where NCCj) equals the number of control variables for phase j. The discrete design parameters are, D = ( q j : p = l,NDP(j), j = 1,NP) where NDPG) equals the number discrete or design variables in phase j. The intervals (phases) are established as, [Ej, Q+l], j=l,NP At each event, boundary conditions of the form tj E can be imposed. In addition, discontinuities in the states can be imposed at the events Denote the collection of events by E = (El, E2,...,E ~ p + l )the , collection of states by x = (xl, x2,..., xNP+l), and the collection of controls by u = (ul, u2, ..., uNP+l). The optimal solution to Equations (1) to (5) is defined to be the set of vector time histories x(t), u(t) and vector parameters D and E which satisfy (1) to (5) and minimize the payoff function, J, which without loss of generality is assumed to depend only on the boundary values (i.e., conditions at events). Integral payoff functions are easily handled by introducing new states. We choose to always minimize the payoff. To maximize $, minimize J equal - $. Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 Outimal Control Solution Method The method we use to solve an optimal control problem is outlined below. States and controls are represented by piecewise polynomials. Differential equations (1) are always satisfied at nodes (Hermite Interpolation). Implicit integration is used to enforce the differential equations at the segment centers. Nonlinear programming is used to: minimize the objective; integrate the differential equations; satisfy the boundary conditions; and impose physical constraints. The free parameters are: the states and controls at nodes (X,U); design variables (D); events (or phase times) (E). When using implicit integration the phases are subdivided into a sequence of segments. The end points of segments are called nodes. All events are nodes but only a few (in general) nodes are events. This situation is illustrated in Figure 2. Within each segment, all states and controls are represented by polynomials in time. The time length of each segment is a user specified fraction of its corresponding phase time length. Segment lengths should be selected so that shorter segments occur where accelerations are high and longer segments are placed where they are low. A consequence of implicit integration is that an accurate integrated trajectory is not obtained until the optimization code has converged. AlAA 90-2829 Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 Stages t X Events E E2 Time --b E3 Figure 2 Piecewise Polynomial Representation We may integrate the equations of motion by quadratures of the form, [ T x= f(x, t) dt We impose the following conditions: x is a cubic polynomial in time on each segment, f is to be evaluated at nodes and segment centers only, any conditions to be satisfied must be local, i.e. involve only adjacent nodes. If x is known at node 1, then at node 2 we must have, x2 = xl + I' f(x, t)dt where T (=Tpt ) is the time length of the segment. 6 E4 If x is given at the first node, the above equation applied to each segment defines a set of conditions which when satisfied, define x as the integral o f f over one phase. Linear conditions are Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 used to hook the phases together. Since x is to be cubic in time, it is sufficient that p be quadratic in time. Function values a two nodes and the segment center are sufficient to define a quadratic in normalized time E =t /T on the segment, This integrates to, A little algebra shows that we must have, where where F1, F2 and FC are the values o f f at the left and right nodes and the segment center. Thus we have, Moving x2 to the right in the above gives a nonlinear equation which is to be driven to zero. We can multiply this by a constant 312T. Doing this, rearranging the F terms and calling the result A gives We call A the "defect" at the segment center. If the cubic polynomial is capable of representing the solution on the given segment, then selecting x l , x2 and Tsi to drive A to zero will produce an accurate solution to equation (I). Note that A can also be defined as the difference between ~ ( X C ) and XC'. AlAA 90-2829 The above procedure is illustrated in Figure 3. The defects for each state are evaluated at the center of each segment and equated to zero. This constitutes a set of nonlinear algebraic equations which are functions of the states and controls at each node, the events and the design parameters . The values and slopes of the controls at each node are free. Constraints are imposed to make the second derivative of the controls continuous at interior nodes and zero at the boundary nodes of each phase. Cubic interpolation is used to obtain the controls at the segment centers. The Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 boundary conditions (4) and the constraints (3) evaluated at both the nodes and the centers of the segments provide additional equations (constraints) to be satisfied. Equation ( 5 ) constitutes a set of linear Select XI, X 2 and Ts so that A=O T= TTS Node Segment Center Node Figure 3 Implicit Integration equations which connects all the phases. All of the independent variables are collected into a single vector P defined by Collecting all of the nonlinear equations into a single vector equation yields where Dij is the defect for ith statement at jth segment center. CN is the collection of all nonlinear constraints from (3). AlAA 90-2829 N is the collection of all nonlinear boundary conditions from (5). Note that the payoff function J is a function of only P. The trajectory optimization problem stated above can be expressed as minimize J = c#P) Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 subject to Np is the dimension of P. AP is comprised of all linear relationships from equations (3) and (5). 1 and u are lower and upper bounds. To handle equality constraints, li is set to ui. Single sided inequality constraints are handled by setting li= - or ui= +-. The partial derivatives of the constraints (defects, physical constraints and boundary conditions) with respect to the independent variables (states, controls, phase times and design parameters) are required. This is accomplished by using finite differences. We note that the task of computing the derivatives is greatly simplified by the sparsity of the Jacobian which results from the fact that the defects depend only on the states at adjacent nodes. DISCUSSION We consider the following three methods for describing multiple vehicle problems: (1) dumb targets, (2) linked phases, (3) extended state vector. Dumb Targets Dumb targets are vehicles that move as a specified function of time or some primary state variables. The OTIS code provides three alternative ways of propagating the position vectors for dumb targets: AlAA 90-2829 (1) motion along great circle paths, (2) tabular functions of position and velocity. (3) satellite orbits (two body propagation), Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 Linked Phases The idea of linked phases is that all trajectory problems be described by one or more phases that can be totally independent of each other. Each has its own set of gravitational, aerodynamic and propulsive forces and each may have independent boundary conditions. Information is then provided as to how these phases should be linked together in both time and state. In the simplest situation, the phases may be different stages of a single vehicle (e.g., a 3 stage booster for putting a satellite in orbit). However, much more complex situations are possible. Each phase can be totally different vehicle with position and velocity evolving at the same time or at different times. They can also represent a single vehicle which splits into two or more vehicles (e.g., a flyback booster). The linking conditions specify that the initial time and state for each phase may be: (1) free, (2) fixed, (3) equal to final time or state of another phase, (4) be a function of the final time and state of another phase. (Reference 3) Option (3) can specify, for example, complete continuity in time and state as would occur if all phases represent the same vehicle. Discontinuities in mass alone can represent spent stage jettisons. Option (4) can allow, for example, analytic propagation between phases. This can make certain orbit propagation problems much more efficient since only the thrusting phases need be integrated. Extended State Vector Dumb targets and linked phases have the disadvantage of limited influence of one vehicle on another. For example, neither of these approaches allow simulation of an airplane optimally evading a missile with proportional guidance. For these and other more general multi-vehicle problems, we extend the state vector to represent the position, velocity and mass for two or more vehicles. OTIS allows simulation of two vehicles using this approach. Controls for one vehicle AlAA 90-2829 may be used for optimization while the other must follow some fixed guidance law (e.g., proportional navigation, tabular etc.). EXAMPLES (1) Glider Minimum Time to Intercept This problem is involved with generating the minimum time intercept of a target flying at Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 constant speed, altitude, and heading. The interceptor was an unpowered aircraft with the following initial conditions: V= 20,000 ft/sec h = 90 degrees (due east) y = 0 degrees h= 200,000 ft longitude = 0 degrees latitude = 0 degrees The target state vector at time 0 was V= 1,000 fdsec h = 180 degrees (due south) y = 0 degrees h= 50,000 ft longitude = 30 degrees latitude = 0 degrees The final interceptor position was constrained to be within 1000 feet of the target. This boundary condition was imposed using one of two methods. The first method constrained the slant range between the interceptor and the target to be 1000 feet. This method did not converge well. The second method involves placing constraints on the differences between the interceptor and target's altitude, longitude, and latitude. This approach was superior in terms of convergence. The solution for this problem is shown in figures 4 and 5. The trajectory for this type of problem is rather simple. The interceptor turns and flies to a point in space. AlAA 90-2829 153000- 103oCO. sbo Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 Time, t, sec 120% 2b 4'00 6b 8'00 Time, t, sec - A 5. - 10- 15203 400 sb0 Time, t, sec 6 a, U)* 25 - 0) 2 a -U Angle of Attack. o 2b C 0 4bo Time, t, sec Figure 4 Glider Minimum Time Intercept Trajectory sbo AlAA 90-2829 2000 1500 Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 1000 500 Longitude, 8, degrees -2.3 3 Ib 13 2b 25 3b Longitude, 9, degrees Figure 5 Dumb Target Intercept (2) Flyback Booster This is a booster problem where the final weight injected into a specified orbit is maximized. The vehicle is composed of two stages (booster and orbiter) which are similar in exterior shape and size. The rocket engines are initially crossfed from the booster stage until its fuel and oxidizer are consumed. At that point, the booster is separated and the orbiter (second stage) continues to orbit injection. Two optimal trajectory problems were solved. The first trajectory was involved in maximizing the weight injected into a polar circular orbit whose altitude is 100. n.mi. No considerations were made for the recovery of the booster. The second trajectory was also concerned with the maximization of the final weight into a 100 n.mi. circular orbit, with the constraint that the booster glide back to the launch point. The specified (fixed) boundary conditions for these problems are as follows: Initial condition2 Velocity = 0. feedsecond AlAA 90-2829 Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 Flight Path Angle = 90 degrees Altitude = 0. feet Longitude = - 120. degrees Latitude = 34.5 degrees Weight = 1,249,318 pounds Stage One Separation Weight of Stage One and Stage Two = 67 1,590 pounds Staoe Two Ignition Weight = 532,931 pounds Stage One Start of Glide Weight = 138,659 pounds b g . e Two Bum Out Altitude > 300,OO feet Apogee Altitude > 100. n.mi. Perigee Altitude > 50. n.rni. Orbital Inclination = 90. degrees End of Stage One Glide 200 feetkecond < Velocity < 700. feedsecond Altitude > 0 feet Longitude = - 120. degrees Latitude = 34.5 degrees The trajectory was modeled using four phases. The first phase covers a short vertical rise where the vehicle is pitched over at a rate of -.001 radians per second. The second phase models the the two mated stages. The third phase models the flight of the orbiter. The fourth and final phase models the booster's glide from the staging point to the launch site. The payoff was the weight injected into a circular orbit whose altitude is equal to the the apogee at the end of third phase. The controls for this problem were chosen to be angle of attack and bank angle. The bank angle was constrained to be zero for phases 1,2 And 3. The initial conditions for phase 4 are specified to be the same as the final conditions for phase 2. The final conditions for the glide were constrained using nonlinear constraints on velocity, altitude and flight path angle. Equality constraints were use to specify the values for longitude and latitude. AlAA 90-2829 The final optimal trajectory and several plots displaying the flyback and non-flyback cases are given as Figures 6, 7 and 8. The maximum injected weight for the flyback case was 117,445 lb, while the value for the non-flyback case was 122,356 lb. 4030 3030 Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 2030 1coo 00 Time, t, sec Time, t, sec Figure 6 Booster Trajectories Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 AlAA 90-2829 -ldO ioo 200 rn 500 400 Figure 7 Booster Controls 500 200 Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 AlAA 90-2829 Downrange, d, n.mi Downrange, d, n.mi. Figure 8 Booster Ground Track (3) Optimal Evasion of Missile with Proportional Guidance This problem involves an aircraft evading an incoming missile. The initial conditions for the aircraft and the missile are given below. IMTIAL CONDITIONS Aircraft Missile Altitude = 5,000 feet Altitude = 2,000 feet Velocity = 500 feet/sec Velocity = 2,000 feet/sec Flight Path Angle = 0.0 degrees Weight = 22,000 pounds Heading = -90.0 degrees (W) Latitude 30.0 degrees (N) AlAA 90-2829 Longitude = 50.0 degrees (E) The following optimization problems was solved: Maximize the miss distance by flying an optimal evasive maneuver. The missile chases the aircraft. The aircraft is modeled as the first vehicle and OTIS will determine optimal angle of attack Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 and bank angle to maximize the miss distance at the point of closest approach, which is forced to be the final time as described below. The second vehicle is modeled as the missile and its angle of attack and bank angle are determined by a proportional navigation algorithm. The state vector for this problem contains the states for first vehicle followed by the states for second vehicle followed by first vehicle controls and control rates (18 total - no second vehicle control information). A final constraint on the relative position and velocity ensures that the final point is the point of closest approach. Plots for this problem are shown as Figures 9, and 10. For the initial conditions chosen the missile flies through the Earth. From a practical point of view, the simulation ended when the missile hit the ground. Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 4 0 8 . 8 9 49'. 91 &.93 49'. 95 49.97 Longitude, 8,degrees Longitude, 8, degrees Figure 9 Optimal Missile Evasion 49.99 50.a 4 AIAA 90-2829 P ; I W m Downloaded by UNIVERSITY OF NEW SOUTH WALES (UNSW) on October 27, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.1990-2829 MISlLE a Time, t, sec Time, t, sec Figure 10 Missile Evasion Controls CONCLUSIONS Efficient simulation of multi-vehicle trajectory optimization problems is greatly facilitated by use of a variety of simulation methods. A computer code (OTIS) which uses collocation and nonlinear programming has incorporated many of the procedures described. Three different examples of OTIS solutions are described to demonstrate the different approaches. REFERENCES (1) Hargraves, C.R. and Paris, S.W., "Direct Trajectory Optimization using Nonlinear Programming and Collocation", Journal of Guidance, Control and Dynamics, July-August 1987. (2) Paris, S.W., Ilgenfritz, D. H., and Hargraves, C.R., "Trajectory Design for Hypemelocity Research Vehicles", AIAA paper 88-4344. (3) Conway, B.A., Presentation on Trajectory Optimization Procedures, Kent, Wa. January 1990.