航天器轨迹优化数值方法综述

JOURNAL OF GUIDANCE, CONTROL, AND DYNAMICS Vol. 21, No. 2, March–April 1998 Survey of Numerical Methods for Trajectory Optimization John T. Betts Boeing Information and Support Services, Seattle, Washington 98124-2207 I. Introduction I T is not surprisingthat thedevelopmentof numericalmethods fortrajectory optimization have closely paralleled the exploration of space and the developmentof the digital computer. Space explo- ration provided the impetus by presenting scientists and engineers with challengingtechnicalproblems.The digital computerprovided the tool for solving these new problems.The goal of this paper is to review the state of the art in the Ž eld loosely referred to as trajectory optimization. Presentinga surveyof a Ž eld as diverse as trajectoryoptimization is a daunting task. Perhaps the most difŽ cult issue is restricting the scope of the survey to permit a meaningful discussion within a limited amount of space. To achieve this goal, I made a conscious decision to focus on the two types of methods most widely used today,namely, direct and indirect. I begin the discussionwith a brief review of the underlying mathematics in both direct and indirect methods. I then discuss the complications that occur when path and boundary constraints are imposed on the problem description. Finally, I describe unresolved issues that are the subject of ongoing research. A number of recurrent themes appear throughoutthe paper. First, the aforementioned direct vs indirect is introduced as a means of categorizing an approach. Unfortunately, not every technique falls neatly into one category or another. I will attempt to describe the beneŽ ts and deŽ ciencies in both approaches and then suggest that the techniquesmay ultimatelymerge. Second, I shall attempt to dis- criminate betweenmethod vs implementation.A numericalmethod is usually described by mathematical equations and/or algorithmic logic. Computational results are achieved by implementing the al- gorithm as software, e.g., Fortran code. A second level of imple- mentation may involve translating a (pre ight) scientiŽ c software implementationinto an (onboard)hardwareimplementation.In gen- eral, method and implementationare not the same, and I shall try to emphasizethat fact. Third, I shall focus the discussionon algorithms instead of physical models. The deŽ nition of a trajectory problem necessarilyentails a deŽ nition of the dynamic environment such as gravitational, propulsion, and aerodynamic forces. Thus it is com- mon to use the same algorithm, with different physical models, to solve different problems. Conversely, different algorithms may be applied to the same physicalmodels (with entirelydifferentresults). Finally, I shall attempt to focus on general rather than special pur- pose methods. A great deal of research has been directed toward solving speciŽ c problems. Carefully specialized techniques can ei- ther be very effectiveor very ad hoc.Unfortunately,whatworkswell for a launch vehicle guidance problemmay be totally inappropriate for a low-thrust orbit transfer. John T. Betts received a B.A. degree from Grinnell College, Grinnell, Iowa, in 1965 with a major in physics and a minor in mathematics. In 1967 he received anM.S. in astronautics from Purdue University, West Lafayette, Indiana, with a major in orbit mechanics, and in 1970 he received a Ph.D. from Purdue, specializing in optimal control theory. He joined The Aerospace Corporation in 1970 as a Member of the Technical Staff and from 1977– 1987 was manager of the Optimization Techniques Section of the Performance Analysis Department. He joined The Boeing Company, serving as manager of the Operations Research Group of Boeing Computer Services from 1987–1989.In his current positionasSeniorPrincipalScientist in the AppliedResearch andTechnologyDivision,he provides technical support to all areas of The Boeing Company.Dr. Betts is a Member of AIAA and the Society for Industrial and Applied Mathematics with active research in nonlinear programming and optimal control theory. E-mail: John.T.Betts@boeing.com. Received March 29, 1997; revision received Oct. 7, 1997; accepted for publicationOct. 20, 1997.Copyright c° 1997 by theAmerican Institute of Aeronautics and Astronautics, Inc. All rights reserved. II. Trajectory Optimization Problem Let us begin the discussion by deŽ ning the problem in a fairly general way. A trajectory optimization or optimal control problem can be formulated as a collection of N phases. In general, the inde- pendentvariable t for phase k is deŽ ned in the region t .k/0 · t · t .k/f . For many applications the independent variable t is time, and the phases are sequential, that is, t .k C 1/0 D t .k/f ; however, neither of those assumptions is required.Within phase k the dynamics of the system are described by a set of dynamic variables z D y .k/.t/ u.k/.t/ (1) made up of the n.k/y state variables and the n .k/ u control variables, respectively. In addition, the dynamics may incorporate the n.k/p parameters p.k/, which are not dependent on t . For clarity I drop the phase-dependentnotation from the remaining discussion in this section. However, it is important to remember that many complex problemdescriptionsrequire different dynamics and/or constraints, and a phase-dependent formulation accommodates this require- ment. Typically the dynamics of the system are deŽ ned by a set of ordinary differential equations written in explicit form, which are referred to as the state or system equations, Py D f [ y.t/;u.t/; p; t ] (2) where y is the n y dimension state vector. Initial conditions at time t0 are deŽ ned by Ã0l · Ã[ y.t0/;u.t0/; p; t0] · Ã0u (3) where Ã[ y.t0/;u.t0/; p; t0] ´ Ã0 , and terminal conditions at the Ž nal time t f are deŽ ned by à f l · Ã[ y.t f /; u.t f /; p; t f ] · à f u (4) where Ã[ y.t f /; u.t f /; p; t f ] ´ à f . In addition, the solution must satisfy algebraic path constraints of the form gl · g[ y.t/;u.t/; p; t] · gu (5) where g is a vector of size ng , as well as simple bounds on the state variables yl · y.t/ · yu (6) 193 D ow nl oa de d by Je t P ro pu lsi on L ab or at or y on O ct ob er 2 6, 2 01 2 | ht tp: //a rc. aia a.o rg | D OI : 1 0.2 514 /2. 423 1 194 BETTS and control variables ul · u.t/ · uu (7) Note that an equality constraint can be imposed if the upper and lower bounds are equal, e.g., .gl /k D .gu/k for some k. Finally, it may be convenient to evaluate expressionsof the form t f t0 q[ y.t/;u.t/; p; t ] dt (8) which involve the quadrature functions q. Collectively we refer to those functions evaluated during the phase, namely, F.t/ D f [ y.t/;u.t/; p; t ] g[ y.t/;u.t/; p; t ] q[ y.t/;u.t/; p; t ] (9) as the vectorof continuousfunctions.Similarly, functionsevaluated at a speciŽ c point,such as theboundaryconditionsÃ[ y.t0/,u.t0/; t0] and Ã[ y.t f /; u.t f /; t f ], are referred to as point functions. The basic optimal control problem is to determine the n.k/u - dimensionalcontrol vectorsu.k/.t/ and parametersp.k/ to minimize the performance index J D Á y t .1/0 ; t .1/0 ; y t .1/f ; p.1/; t .1/f ; : : : ; y t .N /0 ; t .N / 0 ; y t .N / f ; p .N /; t .N /f (10) Notice that the objective function may depend on quantities com- puted in each of the N phases. This formulation raises a number of points that deserve further explanation. The concept of a phase, also referred to as an arc by some authors, partitions the time domain. In this formalism the dif- ferential equations cannot change within a phase but may change fromone phase to another.An obviousreason to introducea phase is to accommodatechanges in the dynamics, for example,when simu- latingamultistagerocket.The boundaryof a phase is often calledan event or junction point. A boundary condition that uniquely deŽ nes the end of a phase is sometimescalledan event criterion(and a well- posedproblemcan have onlyone criterionat each event). Normally, the simulation of a complicated trajectorymay link phases together by forcing the states to be continuous, e.g., y[t .1/f ]D y[t .2/0 ]. How- ever, for multipath or branch trajectories this may not be the case.1 The differential equations (2) have been written as an explicit sys- tem of Ž rst-order equations, i.e., with the Ž rst derivative appearing explicitlyon the left-handside,which is the standardconventionfor aerospace applications.Although this simpliŽ es the presentation, it may not be eithernecessaryor desirable;e.g.,Newton’s lawFDma is not stated as an explicit Ž rst-order system!The objective function (10) has been written in terms of quantities evaluated at the ends of the phases, and this is referred to as the Mayer form.2 If the ob- jective function only involves an integral (8), it is referred to as a problem of Lagrange, and when both terms are present, it is called a problem of Bolza. It is well known that the Mayer form can be obtained from either the Lagrange or Bolza form by introducingan additional state variable.However, again this may have undesirable numerical consequences. III. Nonlinear Programming A. Newton’s Method Essentially all numericalmethods for solving the trajectory opti- mization problemincorporatesome typeof iterationwith a Ž nite set of unknowns. In fact, progress in optimal control solutionmethods closely parallels the progressmade in the underlyingnonlinearpro- gramming (NLP) methods. Space limitations preclude an in-depth presentation of constrained optimization methods. However, it is important to review some of the fundamental concepts. For more complete information the reader is encouraged to refer to the books by Fletcher,3 Gill et al.,4 and Dennis and Schnabel.5 The funda- mental approach to most iterative schemes was suggested over 300 years ago by Newton. Suppose we are trying to solve the nonlinear algebraic equations a.x/D 0 for the root x¤. Beginning with an es- timate x we can construct a new estimate Nx according to Nx D xC ®p (11) where the search direction p is computed by solving the linear system A.x/p D ¡a.x/ (12) The n £ n matrix A is deŽ ned by A D @a1 @x1 @a1 @x2 ¢ ¢ ¢ @a1 @xn @a2 @x1 @a2 @x2 ¢ ¢ ¢ @a2 @xn ::: : : : @an @x1 @an @x2 ¢ ¢ ¢ @an @xn (13) When the scalar step length ® is equal to 1, the iteration scheme is equivalent to replacing the nonlinear equation by a linear ap- proximation constructed about the point x. We expect the method to converge provided the initial guess is close to the root x¤. Of course, this simple scheme is not without pitfalls. First, to compute the search direction, the matrix A must be nonsingular (invertible), and for arbitrarynonlinear functionsa.x/ this may not be true. Sec- ond, when the initial guess is not close to the root, the iterationmay diverge.One common way to stabilize the iteration is to reduce the length of the step by choosing ® such that ka.Nx/k · ka.x/k (14) The procedure for adjusting the step length is called a line search, and the functionused to measureprogress(in this casekak) is called a merit function. In spite of the need for caution, Newton’s method enjoysbroad applicability,possiblybecause,when it works, the iter- ates exhibit quadratic convergence.Loosely speaking, this property means that the number of signiŽ cant Ž gures in x (as an estimate for x¤) doubles from one iteration to the next. B. Unconstrained Optimization Let us nowconsideran unconstrainedoptimizationproblem.Sup- pose that we must choose the n variables x to minimize the scalar objective functionF.x/. Necessary conditionsfor x¤ to be a station- ary point are g.x¤/ D rx F D @F @x1 @F @x2 ::: @F @xn D 0 (15) Now if Newton’s method is used to Ž nd a point where the gradi- ent (15) is zero, we must compute the search direction using H.x/p D ¡g.x/ (16) where the Hessian matrix H is the symmetric matrix of second derivativesof the objective function. Just as before, there are pitfalls in using this method to construct an estimate of the solution. First, we note that the conditiongD 0 is necessarybut not sufŽ cient. Thus a point with zero gradient can be either a maximum or a minimum. At a minimum point the Hessian matrix is positive deŽ nite, but this may not be true when H is evaluated at some point far from the solution. In fact, it is quite possible that the direction p computed by solving Eq. (16) will point uphill rather than downhill. Second, there is some ambiguity in the choice of a merit function if a line D ow nl oa de d by Je t P ro pu lsi on L ab or at or y on O ct ob er 2 6, 2 01 2 | ht tp: //a rc. aia a.o rg | D OI : 1 0.2 514 /2. 423 1 BETTS 195 search is used to stabilize the method. Certainly we would hope to reduce the objective function, that is, F.Nx/ · F .x/. However, this may not produce a decrease in the gradient kg.Nx/k · kg.x/k (17) In fact, what have been described are two issues that distinguish a direct and an indirect method for Ž nding a minimum. For an indi- rect method, a logical choice for the merit function is kg.x/k. In contrast, for a direct method, one probably would insist that the objective function is reduced at each iteration, and to achieve this it may be necessary to modify the calculation of the search direc- tion to ensure that it is downhill. A consequence of this is that the region of convergence for an indirect method may be considerably smaller than the region of convergence for a direct method. Stated differently,an indirectmethodmay requirea better initial guess than required by a directmethod. Second, to solve the equationsgD 0, it is necessary to compute the expressionsg.x/! Typically this implies that analytic expressions for the gradient are necessarywhen using an indirect method. In contrast, Ž nite difference approximations to the gradient are often used in a direct method. C. Equality Constraints Suppose that we must choose the n variables x to minimize the scalar objective function F.x/ and satisfy them equality constraints c.x/ D 0 (18) where m · n. We introduce the Lagrangian L.x;¸/ D F.x/ ¡ ¸>c.x/ (19) which is a scalar function of the n variables x and the m Lagrange multipliers ¸. Necessary conditions for the point .x¤;¸¤/ to be a constrained optimum require Ž nding a stationary point of the Lagrangian that is deŽ ned by rx L.x;¸/ D g.x/¡ G>.x/¸ D 0 (20) and r¸L.x;¸/ D ¡c.x/ D 0 (21) By analogywith the development in the preceding sections,we can useNewton’s method to Ž nd the .nCm/ variables.x;¸/ that satisfy the conditions (20) and (21). Proceeding formally to construct the linear system equivalent to Eq. (12), one obtains HL G> G 0 p ¡ N¸ D ¡g ¡c (22) This system requires the Hessian of the Lagrangian HL D r2x F ¡ m i D 1 ¸ir2x ci (23) The linear system (22) is referred to as the Kuhn–Tucker (KT) or Karush–Kuhn–Tucker system. It is important to observe that an equivalentway of deŽ ning the search direction p is to minimize the quadratic 1 2 p >HL pC g>p (24) subject to the linear constraints Gp D ¡c (25) This is referred to as a quadratic programming (QP) subproblem. Just as in the unconstrainedand root solving applicationsdiscussed earlier, when Newton’s method is applied to equality constrained problems, it may be necessary to modify either the magnitude of the step via a line search or the direction itself using a trust region approach. However, regardless of the type of stabilization invoked at points far from the solution, near the answer all methods try to mimic the behavior of Newton’s method. D. Inequality Constraints An important generalization of the preceding problem occurs when inequality constraints are imposed. Suppose that we must choose the n variables x to minimize the scalar objective function F.x/ and satisfy the m inequality constraints c.x/ ¸ 0 (26) In contrast to the equally constrained case, now m may be greater than n. However, at the optimal point x¤, the constraints will fall into one of two classes. Constraints that are strictly satisŽ ed, i.e., ci .x¤/ > 0, are called inactive.The remainingactive constraints are on their bounds, i.e., ci .x¤/D 0. If the active set of constraints is known, then one can simply ignore the remaining constraints and treat the problem using methods for an equality constrained prob- lem. However, algorithms to efŽ ciently determine the active set of constraints are nontrivial because they require repeated solution of the KT system (22) as constraints are added and deleted. In spite of the complications,methods for nonlinear programming based on the solution of a series of quadratic programming subproblems are widely used. A popular implementation of the so-called successive or sequential quadratic programming (SQP) approach is found in the software NPSOL.6;7 In summary, the NLP problem requires Ž nding the n vector x to minimize F.x/ (27) subject to the m constraints cL · c.x/ · cU (28) and bounds xL · x · xU (29) In this formulation, equality constraints can be imposed by setting cL D cU . E. Historical Perspective Progress in the developmentof NLP algorithms has been closely tied to the advent of the digital computer. Because NLP software is a primary piece of the trajectory optimization tool kit, it has been a pacing item in the development of sophisticated trajectory opti- mization software. In the early 1960s most implementations were based on a simple Newton method (11) and (12) with optimization done parametrically, i.e., by hand. The size of a typical application was nDm ¼ 10. In the 1970s quasi-Newton approximations3;5 be- came prevalent.One popular approach for dealing with constraints was to applyan unconstrainedminimizationalgorithmto a modiŽ ed formof theobjective,e.g.,minimize J .x; ½/D F.x/C 1 2 ½c.x/>c.x/, where½ is large.Althoughthose techniqueshavegenerallybeen su- perseded for generaloptimization,curiouslyenough they are funda- mental to the deŽ nition of the merit functionsused to stabilizestate- of-the-art algorithms. A second popular approach for constrained problems, referred to as the reduced gradient approach, identiŽ es a basic set of variablesthat are used to eliminate the activeconstraints, permitting choice of the nonbasic variables using an unconstrained technique.Careful implementationsof this method8¡10 can be quite effective, especially when the constraints are nearly linear and the number of inequalities is small. Most applications in the 1970s and early 1980s were of moderate size, i.e., nDm< 100. Current ap- plications have incorporated advances in numerical linear algebra that exploit matrix sparsity, thereby permitting applications with n;m ¼ 100;000 (Refs. 11–20). IV. Optimal Control A. Dynamic Constraints The optimal control problemmay be interpreted as an extension of the nonlinear programming problem to an inŽ nite number of variables. For fundamental background in the associated calculus of variations the reader should refer to Ref. 21. First let us con- sider a simple problemwith a single phase and no path constraints. D ow nl oa de d by Je t P ro pu lsi on L ab or at or y on O ct ob er 2 6, 2 01 2 | ht tp: //a rc. aia a.o rg | D OI : 1 0.2 514 /2. 423 1 196 BETTS SpeciŽ cally, suppose we must choose the control functions u.t/ to minimize J D Á[ y.t f /; t f ] (30) subject to the state equations Py D f [ y.t/;u.t/] (31) and the boundary conditions Ã[ y.t f /; u.t f /; t f ] D 0 (32) where the initial conditions y.t0/D y0 are given at the Ž xed initial time t0, and the Ž nal time t f is free. Note that this is a very simpliŽ ed version of the problem (2–10), and we have purposely chosen a problemwith only equality constraints.However, in contrast to the previous discussion,we now have a continuous equality constraint (31) as well as a discrete equality (32). In a manner analogousto the deŽ nition of the Lagrangian function (19), we form an augmented performance index OJ D [Á C º>Ã]t f C t f t