{'fraction_non_alphanumeric': 0.0414202392957401, 'fraction_numerical': 0.024392902891717316, 'mean_word_length': 4.884723336006416, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Advances in Trajectory Optimization for Space Vehicle Control

Space mission design places a premium on cost and operational efficiency. The search for new science and life beyond Earth calls for spacecraft that can deliver scientific payloads to geologically rich yet hazardous landing sites. At the same time, the last four decades of optimization research have put a suite of powerful optimization tools at the fingertips of the controls engineer. As we enter the new decade, optimization theory, algorithms, and software tooling have reached a critical mass to start seeing serious application in space vehicle guidance and control systems. This survey paper provides a detailed overview of recent advances, successes, and promising directions for optimization-based space vehicle control. The considered applications include planetary landing, rendezvous and proximity operations, small body landing, constrained attitude reorientation, endo-atmospheric flight including ascent and reentry, and orbit transfer and injection. The primary focus is on the last ten years of progress, which have seen a veritable rise in the number of applications using three core technologies: lossless convexification, sequential convex programming, and model predictive control. The reader will come away with a well-rounded understanding of the state-of-the-art in each space vehicle control application, and will be well positioned to tackle important current open problems using convex optimization as a core technology.


Introduction
Improvements in computing hardware and maturing software libraries have made optimization technology become practical for space vehicle control. The term computational guidance and control (CGC) was recently coined to refer to control techniques that are iterative in nature and that rely on the onboard computation of control actions (Lu, 2017;Tsiotras and Mesbahi, 2017). This paper surveys optimization-based methods, which are a subset of CGC for space vehicles. We consider applications for launchers, planetary landers, satellites, and spacecraft. The common theme across all applications is the use of an optimization problem to achieve a control objective. Generally speaking, the goal is to solve: where J : R n → R is a cost function, C ⊆ R n is a feasible set, and x is an ndimensional vector of decision variables. Optimization is a relevant area of study for modern space vehicle control for two reasons: effectiveness of formulation, and the (emerging) existence of efficient solution methods.
To answer why an optimization formulation is effective, consider the physical and operational constraints on the tasks that recent and future space vehicles aim to perform. Future launchers and planetary landers will require advanced entry, descent, and landing (EDL) algorithms to drive down cost via reusability, or to access scientifically interesting sites (Blackmore, 2016). Instead of landing in open terrain, future landers will navigate challenging environments such as volcanic vents and jagged blades of ice (San Martín and Robertson, 2017;Europa Study Team, 2012). Meanwhile, human exploration missions will likely be preceded by cargo delivery, requiring landings to occur in close proximity (Dwyer-Cianciolo et al., 2019). Motivated by the presence of water ice, upcoming missions to the Moon will target its south pole (NASA Science, 2019), where extreme light-dark lighting conditions call for an automated sensor-based landing (Robinson, 2018). Even for robotic missions, new onboard technology such as vision-based terrain relative navigation requires the satisfaction of challenging constraints that couple motion and sensing. Regardless of whether one achieves the lowest cost (1a) or not, optimization is indeed one of the most compelling frameworks for finding feasible solutions in the presence of challenging constraints (Tsiotras and Mesbahi, 2017).
In orbit, foreseeable space missions will necessitate robotic docking for sample return, debris capture, and human load alleviation (Woffinden and Geller, 2007). Early forms of the capability have been shown on the Japanese ETS-VII, European ATV, Russian Soyuz, US XSS-11, and US DART demonstrators. Most recently, the human-rated SpaceX Crew Dragon performed autonomous docking with the ISS, and the Orion Spacecraft is set to also feature this ability (Stephens et al., 2013;D'Souza et al., 2007). Further development in autonomous spacecraft rendezvous calls for smaller and cheaper sensors as well as a reduction in the degree of cooperation by the target spacecraft. This will require more flexibility in the chaser's autonomy, which is practically achievable using onboard optimization.
The above mission objectives suggest that future space vehicle autonomy will have to adeptly operate within a multitude of operational constraints. However, optimality usually stipulates operation near the boundary of the set of feasible solutions. In other words, the vehicle must activate its constraints (i.e., touch Figure 1: Taxonomy of optimization problems. Going from inside to outside, each class becomes more difficult to solve. Roughly speaking, SOCP is currently the most general class that can be solved reliably enough to be deployed on a space system in short order. the constraint set boundary) at important or prolonged periods of its motion. By virtue of the feasible set C in (1b), optimization is one of the few suitable methods (and is perhaps the most natural one) to directly impose system constraints (Mayne et al., 2000).
The benefit of an appropriate formulation, however, is limited if no algorithm exists to solve Problem 1 efficiently, which means quickly and utilizing few computational resources. Convex optimization has been a popular approach for formulating problems since it enables efficient solution methods. Figure 1 illustrates a taxonomy of optimization problem classes or families, of which convex optimization is a part. The inner-most class in Figure 1 is the linear program (LP). Next comes the quadratic program (QP), followed by the second-order cone program (SOCP). Troves of detail on each class may be found in many excellent optimization textbooks (Rockafellar, 1970;Nocedal and Wright, 1999;Boyd and Vandenberghe, 2004). Roughly speaking, SOCP is the most general class of problems that state-of-the-art algorithms can solve with high reliability and rigorous performance guarantees Domahidi et al., 2013). Beyond SOCP, the semidefinite program (SDP) class enables optimization over the space of positive semidefinite matrices, which leads to many important robust control design algorithms (Skogestad and Postlethwaite, 2005;Boyd et al., 1994). SDP is the most general class of convex optimization for which off-the-shelf solvers are available, and many advances have been made in recent years towards more scalable and robust SDP solvers (Majumdar et al., 2020).
Although convex optimization can solve a large number of practical engineering problems, future space system requirements often surpass the flexibility of "vanilla" convex optimization. Solving nonconvex optimization problems will be required for many foreseeable space vehicles . Thus, extending beyond SDP, we introduce three nonconvex problem classes.
First, one can abandon the convexity requirement, but retain function continuity, leading to the nonlinear program (NLP) class. Here the objective and constraint functions are continuous, albeit nonconvex. Alternatively, one could retain convexity but abandon continuity. This leads to the mixed-integer convex program (MICP) class, where binary variables are introduced to emulate dis-crete switches, such as those of valves, relays, or pulsing space thrusters Achterberg, 2007;Achterberg and Wunderling, 2013). Note that the MICP and NLP classes overlap, since some constraints admit both forms of expression. For example, the mixed-integer constraint: can be equivalently formulated as a nonlinear continuous constraint: In the most general case, nonlinearity and discontinuity are combined to form the mixed-integer nonlinear program (MINLP) class. Since integer variables are nowhere continuous and the corresponding solution methods are of a quite different breed to continuous nonlinear programming, we reserve MINLP as the largest and toughest problem class. Algorithms for NLP, MICP, and MINLP typically suffer either from exponential complexity, a lack of convergence guarantees, or both . Nevertheless, the optimization community has had many successes in finding practical solution methods even for these most challenging problems (Achterberg and Wunderling, 2013;Szmuk et al., 2019b). This paper stands in good company of numerous surveys on aerospace optimization. (Betts, 1998) presents an eloquent, albeit somewhat dated, treatise on trajectory optimization methods. (Trélat, 2012) provides a comprehensive survey of modern optimal control theory and indirect methods for aerospace problems, covering geometric optimal control, homotopy methods, and favorable properies of orbital mechanics that can be leveraged for trajectory optimization. (Tsiotras and Mesbahi, 2017) corroborate the importance of optimization in forthcoming space missions.  survey the various appearances of lossless convexification and sequential convex programming in aerospace guidance methods.  cover extensively the topic of model predictive control for aerospace applications, where Problem 1 is solved recursively to compute control actions. (Mao et al., 2018a) survey three particular topics: lossless convexification, sequential convex programming, and solver customization for real-time computation. (Shirazi et al., 2018) provide a thorough discussion on the general philosophy and specific methods and solutions for in-space trajectory optimization. Recently, (Song et al., 2020) surveyed optimization methods in rocket powered descent guidance with a focus on feasibility, dynamic accuracy, and real-time performance.
This paper contributes the most recent broad survey of convex optimizationbased space vehicle control methods. We consider rockets for payload launch, rocket-powered planetary and small body landers, satellites, interplanetary spacecraft, and atmospheric entry vehicles. However, we do not cover some related topics like guidance of purely atmospheric vehicles (e.g., missiles and hypersonic aircraft), and control of satellite swarms, due to sufficiently unique distinctions. For a start in these areas, we refer the reader to (Palumbo et al., 2010;Murillo and Lu, 2010;Zarchan, 2019;Tewari, 2011) for hypersonic vehicle control, and (Rahmani et al., 2019;Morgan et al., 2012;Scharf et al., 2003;Tillerson et al., 2002) for swarm control.
From the algorithmic perspective, our focus is on convex optimization-based methods for solving the full spectrum of optimization classes in Figure 1. The motivation for focusing on convex methods comes from the great leaps in the reliability of convex solvers and the availability of flight heritage, which gives convex optimization a technology infusion advantage for future onboard and ground-based algorithms Blackmore, 2016). We nevertheless make side references to other important, but not convex optimization-based, algorithms throughout the text. Lastly, this paper discusses algorithms at a high level, and chooses to cover a large number of applications and methods in favor of providing deep technical detail for each algorithm. The goal, in the end, is to expose the reader to dominant recent and future directions in convex optimization-based space vehicle control research.
The paper is organized as follows. Section 2 covers general theory of important optimization methods used throughout spaceflight applications. Section 3 then surveys each space vehicle control application individually. Section 3.1 surveys powered descent guidance for planetary rocket landing. Section 3.2 discusses spacecraft rendezvous and proximity operations, followed by a discussion in Section 3.3 of its close cousin, small body landing. Constrained attitude reorientation is covered in Section 3.4. Section 3.5 surveys endo-atmospheric flight, including ascent and entry. Last but not least, orbit insertion and transfer are surveyed in Section 3.6. We conclude the paper with a perspective on what lies ahead for computational guidance and control. As such, Section 4 highlights some recent applications of machine learning to select problems. This final section also tabulates some of the optimization software tooling now available for getting started in optimization methods for spaceflight applications.
Notation. Binary numbers belong to the set I {0, 1}. Vectors are written in bold, such as x ∈ R n versus y ∈ R. The identity matrix is generally written as I, and sometimes as I n ∈ R n×n in order to be explicit about size. The zero scalar, vector, or matrix is always written as 0, with its size derived from context. The vector of ones is written as 1, with size again derived from context. Starred quantities denote optimal values, for example x * is the optimal value of x. We use (a; b; c) to concatenate elements into a column vector, like in MATLAB. The symbol ⊗ denotes the Kronecker matrix product or quaternion multiplication, depending on context. The positive-part function [x] + max{0, x} saturates negative elements of x to zero. Given a function f (x(t), y(t), t), we simplify the argument list via the shorthand f [t]. Throughout the paper, we interchangeably use the terms "optimization" and "programming", courtesy of linear optimization historically being used for planning military operations (Wright, 2011). When we talk about "nonlinear programming", we mean more precisely "nonconvex programming". Convexity is now known to be the true separator of efficient algorithms, however this discovery came after linear programming already established itself as the dominant class that can be efficiently solved via the Simplex method (Rockafellar, 1993). Finally, "guidance" means "trajectory generation", while "navigation" means "state estimation".

Background on Optimization Methods
This section provides a broad overview of key algorithms for space vehicle trajectory optimization. The main focus is on methods that exploit convexity, since convex optimization is where state-of-the-art solvers provide the strongest convergence guarantees at the smallest computational cost (Nocedal and Wright, 1999;Boyd and Vandenberghe, 2004).
Our algorithm overview proceeds as follows. First, Section 2.1 introduces the general continuous-time optimal control problem. Then, Section 2.2 describes how the problem is discretized to yield a finite-dimensional problem that can be solved on a computer. Following this introduction, Sections 2.3-2.6 overview important algorithms for space vehicle trajectory optimization.

Optimal Control Theory
Optimal control theory is the bedrock of every trajectory optimization problem (Pontryagin et al., 1986;Berkovitz, 1974). The goal is to find an optimal input trajectory for the following optimal control problem (OCP): In Problem 4, x : [0, t f ] → R nx is the state trajectory and u : [0, t f ] → R nu is the input trajectory, while t f ∈ R is the final time (i.e., the trajectory duration). The state evolves according to the dynamics f : R nx × R nu × R → R nx , and satisfies at all times a set of constraints defined by g : R nx × R nu × R → R nc . At the temporal boundaries, the state satisfies conditions provided by a boundary constraint b : R nx × R nx × R → R n b . The quality of an input trajectory is measured by a cost function consisting of a running cost L : R nx × R nu × R → R and a terminal cost L f : R nx × R → R.
Two aspects differentiate Problem 4 from a typical parameter optimization problem. First, the constraints encode a physical process governed by ordinary differential equations (ODEs) (4b). Second, due to the continuity of time, the input trajectory has an infinite number of design parameters. This makes Problem 4 a semi-infinite optimization problem that cannot be directly implemented on a computer. In the following subsections, we provide a brief overview of two approaches for solving this problem, called the direct and indirect methods. Roughly speaking, direct methods discretize Problem 4 and solve it as a parameter optimization problem, while indirect methods attempt to satisfy the necessary conditions of optimality.

Indirect Methods
The maximum principle, developed since the 1960s, extends the classical calculus of variations and provides a set of necessary conditions of optimality for Problem 4 (Pontryagin et al., 1986;Hartl et al., 1995). The maximum principle has found numerous aerospace applications (Longuski et al., 2014).
The indirect family of optimization methods solves the necessary conditions of optimality, which involves a two-point boundary value problem (TPBVP) corresponding to the state and costate dynamics and their boundary conditions. Traditionally, this is solved by a single-or multiple-shooting method. One limitation of these methods is the requirement to specify in advance the time intervals over which constraint (4c) is active (Betts, 1998). Other issues that hinder onboard implementation include poor convergence stemming from a sensitivity to the initial guess, and long computation time.
Despite these challenges, the indirect approach is often the only practical solution method when aspects like numerical sensitivity and trajectory duration rule out direct methods. Low-thrust trajectory optimization, discussed in Section 3.6, is a frequent candidate for the indirect approach since the low thrust-to-weight ratios and long trajectory durations (from weeks to years) create extreme numerical challenges when formulated as a parameter optimization problem.
Most indirect methods in aerospace literature solve only the necessary conditions of optimality for Problem 4. However, nonlinear optimization problems can have stationary points that are not local minima, such as saddle points and local maxima. This has prompted interest in using second-order conditions of optimality to ensure that the solution is indeed a local minimum (Cesari, 1983). At the turn of the century, researchers showed how second-order information can be incorporated in orbit transfer applications (Jo and Prussing, 2000). In the last decade, further work used second-order optimality conditions for orbit transfer and constrained attitude reorientation problems (Caillau et al., 2012b;Picot, 2012).
A promising modern indirect method family relies on homotopy in order to solve the TPBVP (Pan et al., 2016Pan and Pan, 2020;Taheri and Junkins, 2019;Trélat, 2012;Rasotto et al., 2015;Lizia et al., 2014). Homotopy aims to address the aforementioned challenges of slow convergence, initial guess quality, and active constraint specification. The core idea is to describe the problem as a family of problems parametrized by a homotopy parameter κ ∈ [0, 1], such that the original problem is recovered for κ = 1, and the problem for κ = 0 is trivially solved. For example, consider solving a non-trivial root-finding problem: where y ∈ R n and F : R n → R n is a smooth mapping. A (linear) homotopy method will have us define the following homotopy function: where G : R n → R n is a smooth function that has a known or easily computable root y 0 ∈ R n . Popular choices are G(y) = F (y) − F (y 0 ), called Newton homotopy, and G(y) = y − y 0 , called fixed-point homotopy. In nonlinear homotopy, the function Γ(y, κ) is a nonlinear function of κ, but otherwise similar relationships continue to hold. The locus of points (y, κ) where (6) holds is called a zero curve of the rootfinding problem. Success of the homotopy approach relies on the zero curve being continuous in κ on the interval κ ∈ [0, 1], albeit possibly discontinuous in y. Unfortunately, the existence of such a curve is not guaranteed except for a few restricted problems (Watson, 2002;Pan et al., 2018). In general, the loci of points satisfying (6) may include bifurcations, escapes to infinity, and limit points. Furthermore, the solution at κ = 1 may not be unique.
Nevertheless, homotopy methods have been developed to successfully traverse the κ ∈ [0, 1] interval when a zero curve does exist. The essence of the homotopy approach is to judiciously update an intermediate solution (y k , κ k ) so as to follow a κ-continuous zero curve from y 0 to y K , where Γ(y K , 1) = F (y K ) = 0 and K is the final iteration counter. At each iteration, some methods use a Newton-based root finding approach (Pan et al., 2016), while others rely solely on numerical integration (Caillau et al., 2012a). For further details on the homotopy approach, we refer the reader to (Pan et al., 2016).

Direct Methods
Direct methods offer a compelling alternative where one discretizes Problem 4 and solves it as a parameter optimization problem via numerical optimization. The resulting solution in the convex case is usually very close to the optimal continuous-time one. As discussed in the next section, the solution can satisfy (4b) exactly if an exact discretization method is used . The optimization step is most often performed by a primal-dual interior point method (IPM), for which a considerable software ecosystem now exists thanks to 40 years of active development (Nocedal and Wright, 1999;Forsgren et al., 2002;Wright, 2005). Some of this software is listed in Section 4.1.
Thanks to this expansive software ecosystem, and the large research community actively working on numerical optimization algorithms, direct methods may be considered as the most popular approach today. Their ability to "effortlessly" handle constraints like (4c) makes them particularly attractive (Betts, 1998). In the remainder of this paper, our main focus is on direct methods that use convex optimization.
Nevertheless, as mentioned in the previous section, indirect methods are still relevant for problems that exhibit peculiarities such as extreme numerical sensitivity. It must further be emphasized that some of the best modern algorithms have resulted from the combined use of an indirect and a direct approach. Typically an indirect method, and in particular the necessary conditions of optimality, can be used to discover the solution structure, which informs more efficient customized algorithm design for a direct solution method. We will discuss how this "fusion" approach was taken for powered descent guidance and atmospheric entry applications in Section 2.3 and Section 3.5.2 respectively. Last but not least, the maximum principle can sometimes be used to find the analytic globally optimal solution for problems where no direct method can do so (e.g., nonlinear problems). In this case, an indirect method can provide a   Figure 2: Matrix sparsity patterns for the linear dynamics equation F X = GU for the massspring-damper system in Section 2.2.1, using each of the discretization methods discussed in Section 2.2. A salient feature of the ZOH and RK4 methods is their relative sparsity compared to pseudospectral methods. The adaptive collocation in (d) increases sparsity by subdividing time into multiple intervals within which separate collocations are applied, and which are coupled only through continuity at their interface (in this figure, at x 2 ). The nonzero element colors merely serve to visually separate the elements by their value (larger values correspond to warmer colors).
reference solution against which one can benchmark a direct method's performance Sundström and Guzzella, 2009). In summary, indirect and direct methods play complementary roles: the former is a good ground-truth and analysis tool, while the latter is preferred for real-time onboard implementation.

Discretization
To be solvable on a computer, the semi-infinite Problem 4 must generally be reduced to a finite-dimensional problem. This is done by the process of discretization, where the goal is to convert the differential constraint (4b) into a finite-dimensional algebraic constraint. This is especially important for the family of direct methods discussed in Section 2.1.2, which rely on discretization to solve Problem 4 as a parameter optimization problem.
Generally, discretization is achieved by partitioning time into a grid of N nodes and fixing a basis for the state signal, the control signal, or both (Malyuta et al., 2019). The following subsections discuss three popular approaches: an exact discretization based on zeroth-order hold (Section 2.2.2), an approximate discretization based on the classic Runge-Kutta method (Section 2.2.3), and a pseudospectral discretization that is either global or adaptive (Section 2.2.4). We frame the discussion in terms of three salient features: 1) sparsity of the discrete representation of the dynamics (4b), 2) the mechanics of obtaining a continuous-time trajectory from the discrete representation, and 3) the connection, if any, between the discrete solution and the optimal costates of the original Problem 4 derived via the maximum principle. Our goal is to give the reader enough insight into discretization to appreciate the algorithmic choices for spaceflight applications in Section 3. For a more thorough discussion, we defer to the specialized papers (Betts, 1998(Betts, , 2010Kelly, 2017;Conway, 2011;Ross and Karpenko, 2012;Rao, 2010;Agamawi and Rao, 2020;Malyuta et al., 2019;Phogat et al., 2018b).

Example Dynamical System
To ground our coverage of discretization in a concrete application, let us restrict (4b) to a linear time-invariant (LTI) system of the form: Discretization for (7) is easily generalized to handle linearized dynamics of a nonlinear system (like (4b)) about a reference state-input trajectory (x,ū) : R → R nx × R nu . To do so, replaceÃ andB with: and add a residual term r(t) =f [t] −Ã(t)x(t) −B(t)ū(t) to the right-hand side of (7). Note that in this case, (7) generally becomes a linear time-varying (LTV) system. While the zeroth-order hold method as presented in Section 2.2.2 requires linear dynamics, the Runge-Kutta and pseudospectral methods in their full generality can in fact handle the unmodified nonlinear dynamics (4b).
As a particular example of (7), we will consider a simple mass-spring-damper system with n x = 2 states and n u = 1 input: where r is the position, m is the mass, ζ is the damping ratio, ω n is the natural frequency, and f is the force (input). We set m = 1 kg, ζ = 0.2, and ω n = 2 rad s −1 . Furthermore, consider a staircase input signal where f (t) = 1 for t ∈ [0, t step ) and f (t) = 0 for t ≥ t step . We shall use t step = 1 s. The initial condition is r(0) =ṙ(0) = 0. The simulation source code for this example is publically available 1 .
The dynamics (9) can be written in the form (7) by using the state x = (r;ṙ) ∈ R 2 , the input u = f ∈ R, and the Jacobians:

Zeroth-order Hold
Zeroth-order hold (ZOH) is a discretization method that assumes the input to be a staircase signal on the temporal grid. ZOH is called an exact discretization method because, if the input satisfies this staircase property, then the discrete-time system state will exactly match the continuous-time system state at the temporal grid nodes. In practice, ZOH is a highly relevant discretization type because off-the-shelf actuators in most engineering domains, including spaceflight, output staircase commands .
Optimization routines that use ZOH typically consider a uniform temporal grid, although the method generally allows for arbitrarily distributed grid nodes: The input trajectory is then reduced to a finite number of inputs u k ∈ R nu , k = 1, . . . , N − 1, that define the aforementioned staircase signal: It then becomes possible to find the explicit update equation for the state across any [t k , t k+1 ] time interval, using standard linear systems theory (Antsaklis and Michel, 2006): where ∆t k = t k+1 − t k , x k ≡ x(t k ) and Φ(·, t k ) : R → R nx×nx is the state transition matrix. Since we assumed the system to be LTI, we have Φ(t, t k ) = eÃ (t−t k ) where e is the matrix exponential. If the system is LTV, the state transition matrix can be computed by integrating the following ODE: As it was said before, ZOH is an exact discretization method if the input behaves according to (12). The reason behind this becomes clear by inspecting (13), which exactly propagates the state from one time step to the next. This is different from forward Euler discretization, where the update is: For a general non-staircase input signal, however, there is a subtle connection between ZOH and forward Euler discretization. The former does a zeroth-order hold on the input signal, and integrates the state exactly, while the latter does a zeroth-order hold on the output signal (i.e., the time derivative of the system state). Thus, unlike in forward Euler discretization, state propagation for ZOH discretization cannot diverge for a stable system.
The incremental update equation (13a) can be written in "stacked form" to expose how the discrete dynamics are in fact a system of linear equations. To begin, note that writing (13a) for k = 1, . . . , N − 1 is mathematically equivalent to: where X = (x 1 ; x 2 ; . . . ; x N ) ∈ R N nx is the stacked state, U = (u 1 ; u 2 ; . . . ; u N −1 ) ∈ R (N −1)nu is the stacked input, and ⊗ denotes the Kronecker product. Zeros in (16) denote blocks of commensurate dimensions. Clearly, we can then write the discrete dynamics as: The sparsity pattern for (17) using the mass-spring-damper system (9) with N = 4 and t f = 0.6 s is shown in Figure 2a. Both F and G consist largely of zeros, which has important consequences for customizing optimization routines that exploit this sparsity to speed up computation (Malyuta et al., 2019;Dueri et al., 2017).
An initial value problem (IVP) using a fixed x 1 can be solved either by recursively applying (13a), or by solving (17): where F 1 represents the first n x columns of F , F 2: represents the remaining columns, and X 2: = (x 2 ; . . . ; x N ). We use the left pseudoinverse of F 2: , and note that the solution is unique since F 2: has a trivial nullspace. Figure 3 shows an example of applying (13a) to the mass-spring-damper system (9) using t f = 10 s and N = 51. Because the input step at t step = 1 s falls exactly at a time node, the discretization is exact.
To connect ZOH discretization back to Problem 4, it is now possible to write the problem as a finite-dimensional nonlinear parameter optimization:  Figure 3: Simulation of a mass-spring-damper system using four discretization types, whose matrix sparsity patterns are shown in Figure 2. In each case, the simulation final time is t f = 10 s and a temporal grid of N = 51 is used. The pseudospectral simulations are drawn as continuous curves by using (24) to recover a continuous-time trajectory from a finite set of states at the temporal grid nodes.
the optimal solution is only approximately optimal for Problem 4 due to an inexact running cost integration in (19a) and a finite (N − 1)-dimensional basis with which the ZOH input is constructed. Third, the path constraints (4c) are satisfied pointwise in time via (19c), in other words there is a possibility of inter-sample constraint violation, which can nevertheless sometimes be avoided (Açıkmeşe et al., 2008). Finally, in many important special cases Problem 19 is convex, and a multitude of algorithms exploit this fact, as shall be seen throughout this paper.

Runge-Kutta Discretization
The classic Runge-Kutta (RK4) discretization method can be viewed as an advanced version of the forward Euler method (15). Unlike ZOH, which explicitly assumes a staircase input signal, RK4 is a general numerical integration method that can be applied to any state and control signals. As such, it is an inexact discretization method like forward Euler, albeit a much more accurate one. In particular, if we use the uniform temporal grid (11), then the accumulated RK4 integration error shrinks as O(N −4 ), whereas forward Euler integration error shrinks at the much slower rate O(N −1 ) (Betts, 2010;Butcher, 2016).
Directly from the definition of the RK4 method, we can write the following state update equation: where u k+1/2 = 0.5(u k + u k+1 ). By reshuffling terms in (20), we can write the state update in a similar form to (13a): where {A, B − , B + } are constructed from {I,Ã,B} according to (20). Taking inspiration from (16), (21) can be written in stacked form: where the matrix E combines columns in order to share the input values u k for 2 ≤ k ≤ N − 1 (i.e., the internal time grid nodes): Unlike ZOH, RK4 actually uses the input value at the N th time node, hence there is one extra degree-of-freedom (DoF) leading to a slightly larger stacked input, U = (u 1 ; u 2 ; . . . ; u N ) ∈ R N nu . By defining A and B according to (22), the discrete dynamics take the same form as (17). In this case, the sparsity pattern for the same mass-spring-damper example, using N = 4 and t f = 0.6 s, is shown in Figure 2b. Like ZOH, RK4 yields a sparse representation of the dynamics.
Like ZOH, an IVP using the discretized dynamics can be solved either by recursively applying (21), or via a pseudoinverse like (18). An example is shown in Figure 3. Clearly, RK4 is an inexact discretization method. In this case, the interpolated input u k+1/2 in (20c)-(20d) is erroneous just after the input step at t step = 1 s. Increasing the grid resolution will quickly decrease the integration error, at the expense of a larger linear system (17).
When discretized with RK4, Problem 4 looks much like Problem 19, except u N is an extra decision variable and the user may also choose RK4 to integrate the running cost in (4a). However, a subtle but very important difference is that the solution to the discretized problem generally no longer exactly satisfies the original continuous-time dynamics. Thus, although RK4 may be computationally slightly cheaper than ZOH (especially for LTV dynamics, since it does not require computing integrals like (13b)), it is used less often than ZOH or pseudospectral methods discussed next.

Pseudospectral Discretization
A key property of ZOH discretization from Section 2.2.2 is that it does not parametrize the state signal. As a result, numerical integration is required to recover the continuous-time state trajectory from the solution of Problem 19. In trajectory optimization literature, this is known as explicit simulation or time marching (Rao, 2010). An alternative to this approach is to approximate the state trajectory upfront by a function that is generated from a finite-dimensional basis of polynomials: where τ = 2t −1 f t − 1 and φ i : [−1, 1] → R are known as Lagrange interpolating polynomials of degree N −1. Note that the polynomials are defined on a normalized time interval. Since the temporal mapping is bijective, we can equivalently talk about either t or τ .
Given a temporal grid, the Lagrange interpolating polynomials satisfy an isolation property: φ i (τ i ) = 1 and φ i (τ j ) = 0 for all j = i. Hence, the basis coefficients x i correspond exactly to trajectory values at the temporal grid nodes. Moreover, the trajectory at all other times is known automatically thanks to (24). This is known as implicit simulation or collocation. In effect, solving for the N trajectory values at the temporal grid nodes is enough to recover the complete (approximate) continuous-time trajectory. Because (24) approximates the state trajectory over the entire [0, t f ] interval, this technique is known as a global collocation.
Collocation conditions are used in order to make the polynomial obtained from (24) behave according to the system dynamics (7): where the prime operator denotes differentiation with respect to τ (i.e., dφ i /dτ ) and C ⊆ {1, . . . , N } is the set of collocation points (Malyuta et al., 2019). Note that we have already seen a disguised form of (25) earlier for the RK4 method.
In particular, the well-known (20b)-(20e) are essentially collocation conditions. According to the Stone-Weierstrass theorem (Boyd, 1989), (24) approximates a smooth signal with arbitrary accuracy as N is increased. To avoid the so-called Runge's divergence phenomenon, time is discretized according to one of several special non-uniform distributions, known as orthogonal collocations (de Boor and Swartz, 1973). In this scheme, the grid nodes τ k are chosen to be the roots of a polynomial that is a member of a family of orthogonal polynomials. For example, Chebyshev-Gauss-Lobatto (CGL) orthogonal collocation places the scaled temporal grid nodes at the roots of (1 − τ ) 2 c N −1 (τ ) = 0, where c N (τ ) = cos(N arccos(τ )) is a Chebyshev polynomial of degree N . This particular collocation admits an explicit solution: A pseudospectral method is a discretization scheme that approximates the state trajectory using (24), and selects a particular orthogonal collocation for the collocation points (Rao, 2010;Ross and Karpenko, 2012;Kelly, 2017). In fact, the choice of collocation points is so crucial that flavors of pseudospectral methods are named after them (e.g., the CGL pseudospectral method). Given this choice, if the dynamics and control are smooth, the approximation (24) will converge spectrally (i.e., at an exponential rate in N ) to the exact state trajectory (Rao, 2010).
Associated with any C is a differentiation matrix D ∈ R |C|×N such that D ji = φ i (τ j ). Some collocations (e.g., CGL) admit an explicit differentiation matrix, while for others the matrix can be efficiently computed to within machine rounding error via barycentric Lagrange interpolation (Berrut and Trefethen, 2004). Having D available allows us to write the collocation conditions (25) in stacked form: where the stacked state and input have the same dimensions as in RK4. We may thus write the discrete dynamics in the form of (19b) by defining: The sparsity pattern for the mass-spring-damper example, using N = 4 and t f = 10 s, is shown in Figure 2c. This time, due to F the dynamics constraint is not sparse. This has historically been a source of computational difficulty and a performance bottleneck for pseudospectral discretization-based optimal control (Malyuta et al., 2019;Sagliano, 2019).
Unlike for ZOH and RK4, where an IVP can be solved by recursively applying an update equation, pseudospectral methods require solving (27) simultaneously, which yields the state values at the temporal grid nodes all at once. In general, the solution is once again obtained via the pseudoinverse (18). However, some pseudospectral methods such as Legendre-Gauss (LG) and Legendre-Gauss-Radau (LGR) produce a square and invertible F 2: (furthermore, F −1 2: F 1 = 1 ⊗ I nx ). This can be used to write (19b) in an "integral form" that has certain computation advantages (Françolin et al., 2014): Returning to our example of the mass-spring-damper system, Figure 3 shows a simulation using CGL collocation. Like RK4, pseudospectral discretization is an inexact method, and only approaches exactness as N grows large. In this case, the method struggles in particular due to the discontinuous nature of the input signal, which steps from one to zero at t step = 1 s. The control trajectory is not smooth due to this discontinuity, hence the aforementioned spectral convergence guarantee does not apply. Indeed, it takes disproportionately more grid nodes to deal with this discontinuity, than if we were to subdivide the simulation into two segments t ∈ [0, t step ) and t ∈ [t step , t f ], where the pre-discontinuity input applies over the first interval and the post-discontinuity input applies over the second interval (Darby et al., 2010).
This idea is at the core of so-called adaptive, or local, collocation methods (Darby et al., 2010;Sagliano, 2019;Koeppen et al., 2019;Zhao and Shang, 2018). These methods use schemes such as hp-adaptation (h and p stand for segment width and polynomial degree, respectively) in order to search for points like t step , and to subdivide the [0, t f ] interval into multiple segments according to an error criterion. We defer to the above papers for the description of these adaptation schemes. For our purposes, suppose that a partition of [0, t f ] into S segments is available. The th segment has a basis of N polynomials, a set of collocation points C , and is of duration t s, such that S =1 t s, = t f . Each segment has its own version of (28): where the newly defined F is not to be confused with the earlier F 1 and F 2: matrices. The new matrices in (30) are now combined to write a monolithic dynamics constraint (19b). Doing so is straightforward for the input, which can be discontinuous across segment interfaces: The state trajectory, however, must remain continuous across segment interfaces. To this end, the same coefficient x i is used in (24) for both the final node of segment , and the start node of segment + 1. Understanding this, we can write: where E combines the columns of F in a similar way to (23). The net result is that the final n x columns of F sit above the first n x columns of F +1 . An example of the sparsity pattern for an adaptive collocation scheme is shown in Figure 2d using N 1 = N 2 = 2, t s,1 = t step = 1 s, and t f = 10 s. One can observe that a second benefit of adaptive collocation is that it results in a more sparse representation of the dynamics. This helps to improve optimization algorithm performance (Darby et al., 2010;Sagliano, 2019).
Solving an IVP with adaptive collocation works in the same way as for global collocation. An example is shown in Figure 3, where two segments are used with the split occuring exactly at t step . In this manner, two instances of (24) are used to approximate the state trajectory, which is smooth in the interior of both temporal segments. As such, the approximation is extremely accurate and, for practical purposes, may be considered exact in this case.
A final, and sometimes very important, aspect of pseudospectral discretization is that certain collocation schemes yield direct correspondence to the maximum principle costates of the original optimal control problem (Problem 4). This is known as the covector mapping theorem . One example is the integral form (29) for LG and LGR collocation (Françolin et al., 2014). Roughly speaking, the Lagrange multipliers of the corresponding parameter optimization Problem 19 can be mapped to the costates of Problem 4. Note that this requires approximating the running cost integral in (4a) using quadrature weights {w k } N k=1 defined by the collocation scheme: This unique aspect of pseudospectral methods is why some of the optimal control problem solvers in Table 1 at the end of this article, such as DIDO, GPOPS-II, and SPARTAN, are listed as both direct and indirect solvers. In fact, they all solve a discretized version of Problem 4. Nevertheless, they are able to recover the maximum principle costate trajectories from the optimal solution (Ross and Karpenko, 2012;Patterson and Rao, 2014;Sagliano, 2019).

Convex Optimization
We now come to the question of how to actually solve a finite-dimensional optimization problem such as Problem 19. As mentioned in the introduction, this can be done relatively reliably using well established tools if the problem is convex. Convexity has pervaded optimization algorithm design due to the following property. If a function is convex, global statements can be made from local function evaluations. The ramifications of this property cannot be understated, ranging from the guarantee of finding a global optimum (Rockafellar, 1970) to precise statements on the maximum iteration count (Peng et al., 2002). For the purposes of this review, it is sufficient to keep in mind that a set C ⊆ R n is convex if it contains the line segment connecting any two of its points: where Similarly, a function f : R n → R is convex if its domain is convex and it lies below the line segment connecting any two of its points: Countless resources cover the theory of convex optimization, among which are the notable books by (Boyd and Vandenberghe, 2004;Rockafellar, 1970;Nocedal and Wright, 1999). After applying a discretization technique akin to those in Section 2.2, a trajectory design convex optimization problem takes the following form (which is just another way of writing Problem 19): (36a) Figure 4: Illustration of the convex relaxation technique used throughout much of lossless convexification literature for powered descent guidance. Using the maximum principle, lossless convexification proves that the optimal solution (T * (t); σ * (t)) lies on the green boundary of the setŪ .
is an affine function defining the trajectory boundary conditions. If t f is a decision variable, a sequence of Problem 36 instances can be solved using a line search that computes min t f J * (t f ) (Blackmore et al., 2010). The sequence {A k , B k , d k } N −1 k=1 of matrices of commensurate dimensions represents the linear time-varying discretized dynamics (4b). In numerous aerospace applications, including rocket landing and spacecraft rendezvous, the dynamics are at least approximately of this form (Açıkmeşe and Ploen, 2007;de Ruiter et al., 2013).
The path constraints (36c) are where nonconvexity appears most often for a space vehicle trajectory optimization problems. Sometimes the nonconvexity can be removed by a clever reformulation of the problem, a process called convexification. If the reformulation is exact, in other words the solution set is neither reduced nor expanded, the convexification is lossless. One example of lossless convexification that has pervaded rocket landing literature is a thrust lower-bound constraint. Let T (t) ∈ R 3 be a thrust vector, then combustion stability and engine performance dictate the following constraint: The lower-bound is nonconvex, but it was shown that it admits the following lossless convexification (Açıkmeşe and Ploen, 2007): The convexification "lifts" the input space into an extra dimension, as illustrated in Figure 4. Clearly the lifted feasible setŪ is convex, and its projection onto the original coordinates contains the feasible set defined by (37). The backbone of lossless convexification is a proof via the maximum principle that the optimal solution lies on the boundary ofŪ, as highlighted in Figure 4. Thus, it can be shown that the solution using (38) is optimal for the original problem which uses (37).
Another example of lossless convexification comes from the constrained reorientation problem. Let q(t) ∈ R 4 with q(t) 2 = 1 be a unit quaternion vector describing the attitude of a spacecraft. During the reorientation maneuver, it is critical that sensitive instruments on the spacecraft are not exposed to bright celestial objects. This dictates the following path constraint: where M ∈ R 4×4 is a symmetric matrix that is not positive semidefinite, making the constraint nonconvex. However, when considered together with the implicit constraint q(t) 2 = 1, (39) can be losslessly replaced with the following convex constraint (Kim and Mesbahi, 2004): where µ is chosen such that the matrix M + µI is positive semidefinite. Instead of the maximum principle, the proof of this lossless convexification hinges on the geometry of the unit quaternion.

Sequential Convex Programming
Sequential convex programming (SCP) is an umbrella name for a family of nonconvex local optimization methods. It is one of many available tools alongside nonlinear programming, dynamic programming, and genetic algorithms (Floudas and Pardalos, 2009). If lossless convexification is a surgical knife to remove acute nonconvexity, SCP is a catch-all sledgehammer for nonconvex trajectory design. Clearly many aerospace problems are nonconvex, and SCP has proven to be a competitive solution method for many of them Liu and Lu, 2014;Bonalli et al., 2017). This section provides an intuition about how SCP algorithms work as well as their advantages and limitations. The interested reader can find further information in  which provides a comprehensive tutorial on SCP.
At the core, every SCP algorithm is based on the following idea: iteratively solve a convex approximation of Problem 4, and update the approximation as new solutions are obtained. Figure 5 provides an illustration and highlights how SCP can be thought of as a predictor-corrector algorithm. In the forward predictor path, the current solution is evaluated for its quality. If the quality check fails, the reverse corrector path improves the quality by solving a subproblem that is a better convex approximation. Examples of SCP algorithms include cases where the subproblem is linear (Palacios-Gomez et al., 1982), second-order conic (Mao et al., 2018b), and semidefinite (Fares et al., 2002).
Consider Problem 36 for a simple concrete example of the SCP philosophy. Assume that g is the only nonconvex element and that t f is fixed. At the The forward path can be seen as a "predictor" step, while the reverse path that calls the convex optimization solver can be seen as a "corrector" step. Although the test criterion can be guaranteed to trigger for certain SCP flavors, the solution may not be feasible for the original problem.
location 1 in Figure 5, the SCP method provides an iterate in the form of a current trajectory guess {x k ,ū k } N −1 k=1 . In its most basic form, SCP linearizes and relaxes the g function: where α k ∈ R nc is a virtual buffer zone and we defineḡ g(x k ,ū k , t k ), ∆x k x k −x k , and ∆u k u k −ū k . The subproblem solved at location 2 in Figure 5 is then given by: Problem 42 introduces several new elements. The variables η k regulate the size of trust regions around the previous solution, and the weights w tr , w vc ∈ R are set to large positive values that encourage convergence and zero constraint violation. The best choice of p-norm · in (42d) depends on the problem structure. The stopping criterion used in 3 of Figure 5 may be, for example: where is a user-chosen convergence tolerance constant that can be interpreted as a small "numerical error". SCP denotes a family of solution methods and, as such, countless variations of Problem 42 exist. Early versions of SCP for trajectory generation focused on motion kinematics alone (Schulman et al., 2014) or included dynamics but with few convergence guarantees (Augugliaro et al., 2012). Today, a family of methods is emerging with stronger convergence guarantees, including SCvx (Mao et al., 2018b), GuSTO (Bonalli et al., 2019, and penalized trust region (PTR) . Problem 42 exemplifies the PTR method, where the trust region sizes η k are themselves optimization variables that are kept small using a penalty in the cost (42a). PTR is often the fastest method, but its theoretical convergence properties are relatively unexplored. In comparison, SCvx and GuSTO provide a guarantee that the algorithm converges to a locally optimal solution, albeit with potentially non-zero η k and α k . When these variables are zero, however, the solution is locally optimal for the original optimal control problem.
The main algorithmic differences across SCP methods lie in how the convex approximations are formulated, what methods are used to update the intermediate solutions and to measure progress towards optimality, and how all of this lends itself to theoretical analysis. For example, SCvx uses a discrete-time convergence proof while GuSTO uses the continuous-time maximum principle. The main difference with the PTR method is that both SCvx and GuSTO update η k outside of the optimization problem. Interestingly, the PTR method has been observed to yield much faster convergence in practice, and a theoretical explanation recently appeared .

Related Algorithms
In the general context of optimization, SCP belongs to the class of so-called trust region methods (Nocedal and Wright, 1999;Conn et al., 2000). However, SCP is not to be confused with another popular trust region method, sequential quadratic programming (SQP). First of all, SCP solves its subproblems to full optimality. While this increases the number of iterations in the reverse path of Figure 5, it vastly reduces the number of forward passes. Owing to the growing maturity of IPM solvers and the advent of solver customization (Domahidi et al., 2013;Dueri et al., 2014), iterations in the reverse path are relatively "cheap", making the trade-off a good one. Second, SCP requires only first-order problem information, since nonconvexities are handled by a simple linearization such as in (41). On the other hand, SQP is a second-order method that requires the factorization of a Hessian. This raises concerns about matrix positive semidefiniteness and may require computationally expensive techniques such as the BFGS update (Gill and Wong, 2011).
Differential dynamic programming (DDP) is another family of algorithms that, like SCP, is built around the idea of linearization (Mayne, 1966;Jacobson, 1968). More precisely, DDP solves a discrete-time optimal control problem with an additive cost function like the one in (19a). Although DDP falls outside the scope of this survey paper, we will mention that it has major successful applications in space vehicle trajectory optimization and provides an interesting variation of the linearize-and-solve framework of Figure 5. DDP is particularly popular for low-thrust orbit transfer trajectory optimization. For example, the NASA Mystic software used DDP for low-thrust trajectory optimization of the Dawn Discovery mission to the Vesta and Ceres protoplanets of the asteroid belt (Whiffen and Sims, 2001;Whiffen, 2006). Other appearances of DDP include include multi-revolution and multi-target orbit transfer (Lantoine and Russell, 2012a,b), Earth to Moon transfer with an exclusion zone (Pellegrini and Russell, 2020a,b), and low-thrust flyby trajectory planning to near-Earth objects (Colombo et al., 2009).
A disadvantage of the original DDP algorithm is that it is an unconstrained optimization method. This means that while SCP naturally handles state and input constraints like (36c), implementing such constraints is still an active research area for DDP. Most attempts to incorporate constraints make use of penalty, barrier, augmented Lagrangian, and active set methods (Tassa et al., 2014;Xie et al., 2017). Most recently, extensions of DDP were proposed to handle general nonconvex state and input constraints using a primal-dual interior point method (Pavlov et al., 2020;Aoyama et al., 2020).
Another disadvantage of DDP is that it is a second-order method. Like SQP, this makes DDP more computationally expensive than SCP which only requires first-order information. Nevertheless, there is numerical evidence that DDP is faster than SQP (Lantoine and Russell, 2012b). Furthermore, related algorithms have been developed that only use first-order information, such as the iterative linear quadratic regulator (iLQR). The ALTRO software is a popular modern trajectory optimization toolbox based on the iLQR and augmented Lagrangian methods (Howell et al., 2019).

Mixed-integer Programming
Mixed-integer programming (MIP) solves problems where some decision variables are binary. Consider a concrete example in the context of Problem 4. Without loss of generality, suppose that the control input is partitioned into continuous and binary variables: Binary variables naturally encode discrete events such as the opening of valves and relays, the pulsing of space thrusters, and mission phase transitions (Bemporad and Morari, 1999;Sun et al., 2019). Furthermore, binary variables can help to approximate nonlinear gravity, aerodynamic drag, and other salient features of a space vehicle trajectory optimization problem Marcucci and Tedrake, 2019).
To formally discuss how MIP might be relevant for a spacecraft trajectory optimization problem like Problem 4, consider the space vehicle to be an autonomously switched hybrid system (Saranathan and Grant, 2018). In particular, suppose that the vehicle dynamics (4b) and its constraints (4c) are continuous except for the following extra "if-then" condition: where z ∈ R nz is some mixture of inputs, states, and time. In this formulation, the constraint function c : R nz → R nc is activated if the trigger function q : R nz → R n ζ maps to the negative orthant R n ζ <0 . The conditional statement (45) can be formulated as the following set of mixed-integer constraints: where M ∈ R is a sufficiently large positive number. The function σ : I n ζ → I is called the activation function, and it imposes the if-then logic of (45) through (46a)-(46b). When the binary variable ζ i equals one, the ith trigger is activated. Thus, when σ(ζ) = 1, the left-hand side of (45) holds. In fact, q(z) = 0 is also possible, but this case is irrelevant since an optimal solution will not activate the constraint function unnecessarily. Mixed-integer programming can be used to solve Problem 4 in the presence of the constraints (46). Traditional MIP solvers are based on the branch-andbound method (Nemhauser and Wolsey, 1999;Cook et al., 1998). At their core is a divide-and-conquer logic that often, though not always, speeds up the solution process by eliminating large numbers of possible ζ combinations. Modern MIP solvers also improve runtime through methods like pre-solving (which reduces n ζ ), cutting planes (which introduce new constraints to tighten the feasible space), heuristics, parallelism, branching variable selection, symmetry detection, and so on (Achterberg, 2007;Achterberg and Wunderling, 2013). In the worst case, however, MIP runtime remains exponential in n ζ . This is a large hindrance to onboard implementation, since space vehicle hardware is often not able to support the large MIP computational demand (Malyuta et al., 2020;Malyuta and Açıkmeşe, 2020b,a).
As usual in optimization, one can trade the global optimality of MIP for solution speed by accepting local optimality or by approximating the precise statement (45) with a more efficient formulation. In the following subsections, we will introduce two popular approaches that have recently emerged in both direct and indirect solution methods for solving MIP problems without introducing integer variables.

State-triggered Constraints
State-triggered constraints (STCs) take the direct solution approach, and are under active study using the SCP framework from Section 2.4 Malyuta et al., 2020). Roughly speaking, STCs embed the discrete ifthen logic from (45) into a continuous direct formulation with minimal penalty to the solution time . In its most basic form, an STC models (45) for the scalar case n ζ = 1 and n c = 1. While there are several useful theoretical connections between STCs and the linear complementarity problem (LCP) (Cottle et al., 2009), STCs encode a much larger feasible space Figure 6: A Venn diagram of the sets of solution variables z ∈ R nz that satisfy the trigger condition Q, constraint condition C, and all other constraints Z. STCs ensure that feasible solutions satisfy z / ∈ Q ∩ C c (illustrated by the red set). The feasible set with an STC is shaded green, and the feasible set with the bidirectional constraint (47) is shaded in yellow. The sets with bold outlines are padded to help visual separation only.
than LCP constraints . Namely, STCs only encode forward implications, and they are not bidirectional statements of the following form: (47) Figure 6 illustrates the distinction between (45) and (47). The green set denotes the feasible set with the STC, while the yellow set denotes the feasible set of the more restrictive constraint (47). Clearly, the feasible space is larger when using the STC, and this can translate into a more optimal solution.
Continuing our discussion for the scalar case, it can be shown that (45) is equivalent to either one of the following two continuous constraints: where σ ∈ R is a slack variable that plays the role of the activation function from (46c). Although both constraints in (48) are nonconvex, they are readily ingested by the SCP linearization process described in Section 2.4. A notable feature of STCs is that they readily extend to the multivariable case of (45), and have been shown to handle both AND and OR combinations of triggers and constraints (Szmuk et al., 2019a,b): In the general context of optimization, STCs do for trajectory optimization what the S-procedure from linear matrix inequalities (LMIs) does for stability analysis and controller synthesis (Boyd et al., 1994), and what sum-of-squares (SOS) programming does to impose polynomial non-negativity over basic semialgebraic sets (Majumdar and Tedrake, 2017). That is, STCs embed an otherwise difficult logic constraint into a tractable continuous formulation. In particular, note that the scalar version of (45) can be written as: On the other hand, the S-procedure and SOS programming consider the following constraints, respectively: where f i , i = 0, . . . , p, are quadratic functions, while p, p eq , and p ineq are polynomials. Comparing (50) with (51) makes the connection to STCs clear.

Homotopy Methods
Homotopy methods, also known as numerical continuation schemes, were previously introduced in Section 2.1.1 in the context of solving standard optimal control problems. It turns out that homotopy can also be used to encode (45) in a continuous framework, and has been successfully embedded into recent indirect trajectory optimization algorithms. In this section, we briefly introduce the relaxed autonomously switched hybrid system (RASHS) and composite smooth control (CSC) methods (Saranathan and Grant, 2018;Taheri et al., 2020a;Arya et al., 2020).
To begin, let σ denote the activation function from (46c). Using the third equation from (48a), we note that (45) is exactly equivalent to the following constraint: σ(ζ)c(z) = 0.
At the core of the RASHS and CSC methods is an approximation of the binary function σ by a continuous sigmoid functionσ : R n ζ → [0, 1]: where the latter equation uses the theory of hyperbolic tangent smoothing (Taheri and Junkins, 2018). The homotopy parameter κ ∈ [0, ∞) regulates the accuracy of the approximation, with increasing accuracy as κ grows, such that lim κ→∞σ (q) = σ(ζ). Figure 7 illustrates how both sigmoid functions evolve as κ increases. The core idea of RASHS and CSC is to begin with a   (53). A sweep is shown from homotopy parameter κ = 1 (lighter colors) to κ = 100 (darker colors). The blue dashed curve shows the alternative sigmoid for κ = 1 (i.e., CSC for (a) and RASHS for (b)). As κ increases, the sigmoid quickly converges to an accurate approximation of the binary activation function in (46c). While the nature of both sigmoids is similar, for a given κ the CSC sigmoid is more localized around the y-axis, and hence is a closer approximation of a step signal.
small κ where the optimal control problem is continuous and "easy" to solve, and to judiciously increase κ to such a large value that the solution becomes indistinguishable from its MIP counterpart. It is worth noting the specific instances of (45) considered by RASHS and CSC. The former method was developed to compute time-or state-triggered multiphase trajectories, where vehicle dynamics change across phases (e.g., stage separation during rocket ascent) (Saranathan and Grant, 2018). Such a system is also known as a differential automaton (Tavernini, 1987). In this case, we can have m constraints of the form (45), where the kth constraint is: and q k : R → R n k ζ indicates the time interval where the kth dynamics apply. Assuming that the time intervals do not overlap, we can sum the smoothed versions of (54) to obtain a single continuous system dynamics constraint: Note that the new dynamics (55) are a convex combination of the individual dynamics over the m time intervals. As κ is increased, the approximation becomes more accurate, and the correct f k functions begin to dominate their respective intervals. Moreover, using (55) instead of (54) has the algorithmic advantage of replacing a multi-point BVP with a TPBVP.
The CSC method, on the other hand, considers systems with fixed dynamics but multiple control inputs or constraints (Taheri et al., 2020a). In both cases, the overall control input can be expressed as a function of m "building block" inputs u k , such that: where q k : R nx × R nu × R → R n k ζ are mutually exclusive indicators of when the kth building block input applies. Like for RASHS, the following equation provides a smooth approximation of the control, from which CSC derives its name: Clearly, RASHS, CSC, and STCs are all approaching the same problem of efficiently handling (45) from subtly different angles. It is worth noting that for the moment, both RASHS and CSC can only handle the AND combination (49a) of trigger and constraint functions. Most recently,  showed that a similar homotopy framework can handle OR combinations of trigger functions. This opens up an interesting research avenue to develop a unifying homotopy method that handles all the logical combinations in (49).

Model Predictive Control
The preceding sections focused on solving one instance of Problem 4. We now place ourselves in the context of a control system whose architecture is illustrated in Figure 8. Two important algorithm categories that are part of a control system are so-called feedforward and feedback (Lurie and Enright, 2000), and optimization-based methods can potentially be applied to both. In the feedback path, the current state estimate of the system is used to continually update the control signal, which means that Problem 4 must be re-solved many times. This is the core idea of model predictive control (MPC).
In its most basic form, an MPC formulation of Problem 36 can be expressed as follows: (58d) Figure 9 illustrates how Problem 58 can be used to control a dynamical system. Note that Problem 58 is a parametric optimization problem because it depends on the current state estimatex ∈ R nx . The first optimal control input u * 1 for Problem 58 becomes u(t) in Figure 8. The weight matrices Q 0 and R 0 in the running cost and the terminal weight matrix Q f 0 are chosen to get a desired response. Together with the terminal constraint (58d), these choices must ensure stability and recursive feasibility in closed-loop operation (i.e., the problem must be feasible the next time that it is solved).
The main advantage of MPC is that it is arguably the most natural methodology for handling system constraints in a feedback controller (Mayne et al., 2000). However, because MPC operates in a feedback loop, stability and performance are both critical and strongly dependent on uncertainty robustness and execution frequency (Lurie and Enright, 2000;Skogestad and Postlethwaite, 2005). Troves of information have been compiled on the subject, which remains an active research area. Numerous surveys on MPC cover general and future methods (Mayne, 2014), robustness (García et al., 1989;Bemporad and Morari, 2007;Mayne, 2015), computational requirements (Alessio and Bemporad, 2009), and industrial applications Mao et al., 2018a;Di Cairano and Kolmanovsky, 2018;Qin and Badgwell, 2003). For space vehicle applications in particular, where onboard computation is limited, we single out so-called explicit MPC (Bemporad et al., 2002;Rawlings et al., 2017;Borrelli et al., 2017). The concept is to pre-compute a lookup table for the solution of Problem 58. This turns out to be possible to do exactly when the MPC problem is a QP, and approximately in more general cases up to MICP . When onboard storage and problem dimensionality permit, explicit MPC yields a much faster and computationally cheaper algorithm in which onboard optimization is replaced by a static lookup table.

Applications
This section describes the application of optimization methods from the previous section to state-of-the-art space vehicle control problems. The following subsections cover the following key areas of spaceflight. Section 3.1 discusses rocket powered descent for planetary landing. Section 3.2 covers spacecraft rendezvous and Section 3.3 covers a closely related problem of small body landing. Section 3.4 talks about attitude reorientation. Endo-atmospheric ascent and entry are surveyed in Section 3.5. Last but not least, orbit transfer is discussed in Section 3.6.

Powered Descent Guidance for Rocket Landing
Powered descent guidance (PDG) is the terminal phase of EDL spanning the last few kilometers of altitude. The goal is for a lander to achieve a soft and precise touchdown on a planet's surface by using its rocket engine(s). PDG technology is fundamental for reducing cost and enabling access to hazardous yet scientifically rich sites (Starek et al., 2015;Steinfeldt et al., 2010;Braun and Manning, 2006;Jones, 2018;Robertson, 2017;Europa Study Team, 2012;NASA Science, 2019;Robinson, 2018). The modern consensus is that iteration-based algorithms within the CGC paradigm, rather than System MPC Disturbance Reference (e.g. 0) State estimatê x at time t 1 Control input u * 1 over Figure 9: Block diagram illustration of an MPC controller. At each time step t k , MPC computes the optimal control input u * 1 by using a mathematical model of the system and solving Problem 58, which is a receding horizon optimal control problem. Note the three states drawn in the diagram: the actual statex, the estimated statex, and the internally propagated MPC state x. Each state may be slightly different due to estimation error, model uncertainty, and disturbances.
closed-form solutions, are required for future landers (Lu, 2017;. The survey of applications in this section demonstrates that optimization offers a compelling iteration-based solution method due to the availability of real-time algorithms that can enforce relevant PDG constraints. To place state-of-the-art PDG into context, let us briefly mention some key heritage methods. Initial closed-form algorithms are known as explicit guidance, which is characterized by directly considering the targeting condition each time the guidance command is generated (Lu, 2020). Early algorithms solved a version of the following OCP: Here, r(t) ∈ R 3 denotes position, a(t) ∈ R 3 is the acceleration control input, g ∈ R 3 is the gravitational acceleration vector and t f is the flight duration. Position and velocity boundary values are fixed. The optimal solution to Problem 59 is known as the E-Guidance (EG) law (Cherry, 1964;D'Souza, 1997): where t go t f − t is the time-to-go and: are the zero-effort-miss and zero-effort-velocity terms (Furfaro et al., 2011;Song et al., 2020). Nominally, (60) results in an affine acceleration profile. If instead one allows the acceleration profile to be quadratic, an additional DoF appears, which is fixed by setting the terminal acceleration a(t f ) = a f . This results in the Apollo guidance (APG) law, which flew on the historic Lunar missions (Klumpp, 1974): The concept of an acceleration profile behind EG and APG has since been extended and generalized, resulting in a polynomial guidance family of algorithms.  augment the cost (59a) with a surface collision-avoidance term. (Guo et al., 2013) formulate a QP to solve for an intermediate waypoint that augments collision-avoidance capabilites as well as enforces actuator saturation for thrust-and power-limited engines. (Lu, 2019) develops a general theory for polynomial guidance laws that contains EG and APG as special cases. For one of the best modern explanations of polynomial guidance methods, the reader should consult (Lu, 2020). Unfortunately, closed-form polynomial guidance is unable to handle many operational constraints (Lu, 2018) and is not fuel optimal since the cost (59a) rather penalized energy.
To overcome these limitations, research has long sought to characterize and eventually solve the more general fuel-optimal PDG problem. The first milestone towards fuel-optimal PDG was a closed-form single-DoF vertical descent solution (Meditch, 1964), illustrated in Figure 12a. Evidence suggests that Apollo landings came close to this optimum (Klumpp, 1974;Mindell, 2008). The maximum principle (Pontryagin et al., 1986) played a key role back then, and continues to do so in the present day.
Seeking to generalize the single-DoF result, Lawden formulated the necessary conditions of optimality for 3-DoF PDG (Lawden, 1963;Marec, 1979). However, solving the necessary conditions requires shooting methods, which are typically too computationally expensive and sensitive to the initial guess to allow efficient onboard implementation (Betts, 1998). More recently, (Topcu et al., 2005(Topcu et al., , 2007 extended the results from (Lawden, 1963) to the case of angular velocity control, and compared the solution quality of fuel-optimal 3-DoF PDG to the necessary conditions of optimality. However, the aim of the work was not realtime onboard implementation, so nonlinear programming (SQP) was used via the GESOP solver.
After decades of research into problem characterization, a watershed moment for problem solution came in the mid 2000s with the papers Ploen, 2005, 2007). The authors solved the following 3-DoF PDG problem, illustrated in Figure 10, via the process of lossless convexification described in Section 2.3: (63a) Figure 10: Illustration of the basic powered descent guidance solved by Problem 63 via lossless convexification. The goal is to safely bring the rocket lander to standstill on the landing pad while satisfying the thrust magnitude constraints and maintaining a minimum glideslope.
Unlike the classical Problem 59, Problem 63 readily handles several important operational constraints, including thrust bounds (63d) and glide slope (63e). Through numerical simulations for a prototype Mars lander, (Açıkmeşe and Ploen, 2007) confirmed that the optimal thrust has a max-min-max profile as shown in Figure 11. This profile was proven to be optimal for 3-DoF PDG in (Lawden, 1963;Topcu et al., 2007).
Over the course of the next decade, the method was expanded to handle fairly general nonconvex input sets (Açıkmeşe and Blackmore, 2011), minimumerror landing and thrust pointing constraints (Blackmore et al., 2010;Carson III et al., 2011a;Açıkmeşe et al., 2013c), classes of affine and quadratic state constraints (Harris and Açıkmeşe, 2013b, classes of nonlinear (mixed-integer) dynamics , certain binary input constraints Harris, 2021), and conservative conic obstacles (Bai et al., 2019).
The maturity of a method can be gauged by the availability of a precise statement of its limits, similar to the role played by the Bode integral in frequencydomain control (Skogestad and Postlethwaite, 2005;Lurie and Enright, 2000). Such a characterization appeared for lossless convexification in the form of constrained reachable or controllable sets (Eren et al., 2015b;Dueri et al., 2016;Dueri, 2018) or "access" conditions (Song et al., 2020). These sets, obtained numerically and with arbitrarily high precision, define the boundary conditions for which versions of Problem 63 are feasible.
The practicality of lossless convexification-based PDG methods was demon-  (Meditch, 1964); (b) corresponds to 3-DoF translation-only landing from (Lawden, 1963;Açıkmeşe and Ploen, 2007); (c) corresponds to planar landing with rotation from . The thrust profile for general 6-DoF PDG with translation and rotation is an open problem. In particular, there is no theory to guarantee that it should be bang-bang, thus (d) shows a profile with no clear structure.  (Scharf et al., 2015), and complete off-line lookup table generation (Açıkmeşe et al., 2008). The resulting algorithm, G-FOLD , solves a tailored version of Problem 63 using a custom C-language SOCP solver called Bsocp . G-FOLD is able to compute rocket landing trajectories in 100 ms on a 1.4 GHz Intel Pentium M processor. Further evidence of real-time performance was presented by , where Bsocp ran on a radiation-hardened BAE RAD750 PowerPC. Despite the significant flight envelope expansion afforded by lossless convexification Wolf et al., 2012;Carson III et al., 2011b), an inherent limitation of 3-DoF PDG is that the computed trajectory cannot incorporate attitude constraints other than those on the thrust vector, which serves as an attitude proxy. An extensive simulation campaign is required to validate the 3-DoF trajectory to be executable by a fundamentally 6-DoF lander system Kamath et al., 2020). Thus, recent PDG research has sought 6-DoF formulations that are able to incorporate attitude dynamics and constraints.
The SCP family of methods, discussed in Section 2.4, has emerged as an effective approach to transition from a fully convex 3-DoF problem to a 6-DoF problem with some nonconvexity. Some of the popular SCP algorithms include SCvx (Mao et al., 2018b), penalized trust region , and GuSTO (Bonalli et al., 2019. Some other algorithms based around similar ideas have also emerged, such as ALTRO which is based on iterative LQR (Howell et al., 2019).
A vast number of flavours of SCP exist, however, since it is a nonlinear optimization technique that works best when tailored to exploit problem structure. In certain cases, lossless convexification is embedded to remove some nonlinearity. (Liu, 2019) convexifies an angle-of-attack (AoA) constraint relating to an aerodynamic control capability, (Simplício et al., 2019) solve a version of Problem 63 in a first step and passes the solution to a second step involving SCP, while (Li et al., 2020;Wang et al., 2019a;Szmuk et al., 2016) use the classical convexification result for the thrust magnitude constraint (38).
Since the mid 2010s, SCP technology enabled the expression of quadratic aerodynamic drag and thrust slew-rate constraints (Szmuk et al., 2016), attitude dynamics (Szmuk et al., 2017), variable time-of-flight , and an ellipsoidal drag model that allows aerodynamic lift generation along with variable ignition time . Several papers on SCP "best practices" have also appeared, including thrust input modeling (Szmuk et al., 2017), the effect of discretization on performance (Malyuta et al., 2019), and using dual quaternions to alleviate nonconvexity in the constraints by offloading it into the dynamics . Practical details on real-time implementation are also available , where the SCP solution is compared to the globally optimal trajectory for a planar landing problem . Most recently, a comprehensive tutorial paper with open-source code was published, and describes the algorithmic and practical aspects of SCP methods and of lossless convexification . Figure 12 summarizes the dominant directions of PDG development since 2005. Starting from the classical vertical-descent result by (Meditch, 1964), Figure 12a, the early breakthrough for practical onboard PDG solution was achieved in 2007 by (Açıkmeşe and Ploen, 2007), Figure 12b. Since then, 3-DoF PDG methods have been extended and flight tested, Figure 12c. In particular, more complicated effects such as aerodynamic drag force were added by these extensions, which are listed in the preceding paragraph. Perhaps the biggest modern shift in PDG technology development has been to consider attitude dynamics, which is motivated by the inability to impose non-trivial attitude constraints in a 3-DoF formulation . This has led to a family of so-called 6-DoF PDG algorithms, Figure 12d, that often rely on SCP methods. To compare how close SCP comes to the global optimum, recent work found optimal solutions for "planar" PDG (Reynolds and Mesbahi, 2020b), Figure 12e. This work restricts the landing trajectory to a 2D plane, but does include attitude dynamics. Therefore it represents both a generalization of Figure 12b and a restriction of Figure 12d, and provides new insight into the 6-DoF PDG optimal solution structure. Today, PDG research evolves along the following broad directions: guaranteeing real-time performance, convergence, Figure 12: Progression of PDG problem complexity. The red, green, and blue arrows denote thrust, torque, and aerodynamic force respectively. The red region denotes the feasible thrust set. and solution quality, handling binary constraints, and incorporating uncertainty as shown in Figure 12f.
One exciting development for SCP in recent years has been the advent of state-triggered constraints, introduced in Section 2.5.1. This allows real-time capable embedding of if-then logic into the guidance problem. To demonstrate the capability,  imposed a velocity-triggered AoA constraint,  imposed a distance-triggered line-of-sight constraint,  imposed a collision-avoidance constraint, and  imposed a slant-range-triggered line-of-sight constraint. In particular, the latter two works develop a theory of compound STCs that apply Boolean logic to combine multiple trigger and constraint functions, as shown in (49). The impact of STCs on the ability to compute solutions in real-time is discussed in Reynolds et al., 2019a).
Simultaneously with the development of SCP for PDG, the pseudospectral discretization community has produced a rich body of work investigating the solution quality benefits of that method. Building on foundational early work (Rao, 2010;Fahroo and Ross, 2002;Garg et al., 2010;Kelly, 2017), it was demonstrated for a variant of Problem 63 that pseudospectral methods yield greater solution accuracy with fewer temporal nodes (Sagliano, 2018b). However, as discussed in Section 2.2.4, pseudospectral methods traditionally yield slower solution times because they generate non-sparse matrices for the discretized equations of motion (Malyuta et al., 2019). By using an hp-adaptive scheme inspired by the finite element method (Darby et al., 2010), it was shown that this can be somewhat circumvented (Sagliano, 2018a(Sagliano, , 2019. Furthermore, it was shown that pseudospectral discretization within an SCP framework yields solutions up to 20 times faster than using sequential quadratic programming (Wang and Cui, 2018). As deterministic PDG algorithms mature, research is becoming increasingly interested in making the trajectory planning problem robust to various sources of uncertainty. One approach is to design a feedback controller to correct for deviations from the nominal trajectory, such that the overall control input is given by: wherex(t) andū(t) are the nominal state and control respectively, and K(t) ∈ R nu×nx is a feedback gain matrix. In (Ganet-Schoeller and Brunel, 2019; Scharf et al., 2017), the feedback controller is designed separately from the nominal trajectory. However, incorporating feedback law synthesis into the nominal trajectory generation problem can achieve more optimal solutions (Garcia-Sanz, 2019). This "simultaneous" feedback-feedforward design was done via multidisciplinary optimization in , desensitized optimal control in (Shen et al., 2010;Seywald and Seywald, 2019), chance-constrained optimization in (Ono et al., 2015), and covariance steering in Tsiotras, 2018, 2019). Other work in this domain includes open-loop robust trajectory design via Chebyshev interval inclusion (Cheng et al., 2019c), and a posteriori statistical analysis through linear covariance propagation  and Monte Carlo simulation . PDG methods based on lossless convexification and SCP are in most cases implicit guidance methods. In this setup, the targeting condition (e.g., soft touchdown on the landing pad) is met by tracking a reference trajectory that yields the correct terminal state. Functionally, PDG methods are most often situated in the FF block of Figure 8, and they generate a complete trajectory upfront that is tracked by a feedback controller. From a systems engineering perspective, this has a clear advantage of allowing heritage control methods to perform the intricate and critical control of the actual vehicle. However, it was mentioned at the start of this section and in Section 2.6 that continually resolving for the PDG trajectory can offer additional robustness. In contrast to traditional polynomial guidance, some modern approaches aim to leverage this robustness and also satisfy system constraints via model predictive control. (Cui et al., 2012) show how to leverage MPC for landing with an uncertain state and variable gravitational field, while  show how to ensure recursive feasibility and a bounded guidance error by executing a nominal and relaxed optimization problem in parallel. In both methods, the full trajectory optimization problem is solved from the current state to the final landing location, thus the MPC horizon "shrinks" throughout the PDG maneuver. A more traditional approach is taken by (Lee and Mesbahi, 2017), where the prediction horizon extends for a finite duration beyond the current state. The authors also show that difficult constraints on sensor line-of-sight and spacecraft Figure 13: Illustration of a basic rendezvous scenario. Roughly speaking, the goal is for a chaser spacecraft to use the thrust T from its reaction control system to dock with a target while avoiding debris and respecting constraints such as staying within the approach cone.
attitude are convex using a dual quaternion representation. Numerical performance of MPC for PDG on an embedded ARM platform was documented in (Pascucci et al., 2015).

Rendezvous and Proximity Operations
Let us now switch contexts from the final stages of planetary landing to the realm of orbital spaceflight. A key task for a spacecraft in orbit is to perform rendezvous and proximity operations (RPO). The goal is to bring an actively controlled chaser vehicle and a passively controlled target vehicle to a prescribed relative configuration, in order to achieve mission objectives such as inspection or docking. A detailed overview of RPO history and technology development can be found in (Fehse, 2003;Goodman, 2006;Woffinden and Geller, 2007;Luo et al., 2014). This section focuses on the challenges and developments in RPO using convex optimization-based solution methods.
Throughout this section we consider the following RPO trajectory optimization problem, illustrated in Figure 13: where r(t),r(t) ∈ R 3 denote the positions of the chaser and target spacecraft in the inertial frame. The basic objective in (65a) is to minimize fuel consumption (Park et al., 2013). Other choices include sparsification of the control sequence (Hartley et al., 2013), trading off flight duration with fuel consumption (Hu et al., 2018), encouraging smoothness of the control sequence (Li and Zhu, 2018a), and reducing the sensitivity to sensing and control uncertainties (Jin et al., 2020). We note that Problem 65 only characterizes the last phase of RPO. The reader is referred to (Hartley et al., 2012;Sun et al., 2019) for examples of multi-phase RPO trajectory optimization.
Since the quantity of interest in RPO is the relative motion between the chaser and the target, it is commonplace to express the dynamics (65b) in a different reference frame. Examples include the local-vertical local-horizontal (LVLH) frame centered at the target, or a line-of-sight polar reference frame (Li and Zhu, 2017). Based on this choice, different models of relative dynamics have been studied, and are surveyed in (Sullivan et al., 2017). For near-circular orbits, linear time-invariant Hill-Clohessy-Wiltshire (HCW) equations are the most popular model (Clohessy and Wiltshire, 1960). For elliptical orbits, the linear time-varying Yamanaka-Ankerson (YA) state transition matrix is the usual choice (Yamanaka and Ankersen, 2002). Perhaps a cleaner approach is to avoid relative dynamics by working in the inertial frame, as done in (65b). (Lu and Liu, 2013;Liu and Lu, 2014) showed that fast and reliable trajectory optimization is still possible in this case, by applying the same lossless convexification as in Problem 63 to the constraints (65c) and (65d) and successively linearizing the dynamics (65b). (Benedikter et al., 2019b) further proposed a filtering technique for updating the linearization reference point to improve the algorithm robustness. The advantage of this approach is its compatibility with general Keplerian orbits and perturbations like J 2 harmonic and aerodynamic drag.
One key challenge in RPO is to avoid collision with external debris or part of the target vehicle itself, which is described by constraint (65e). One approach to enforcing (65e) is to pre-compute a so-called virtual net of trajectories that allows to avoid obstacles in real-time via a simple graph search (Frey et al., 2017;Weiss et al., 2015b). The pre-computation procedure, however, may be prohibitively computationally demanding. In comparison, solving Problem 65 directly can avoid virtual net construction altogether if an efficient solution method is available. To this end, the keep-out zone B(t) is usually chosen to be a polytope, an ellipsoid, or the union of a mix of both if multiple keep-out zones are considered (Hu et al., 2018). As shall be seen below, polytope approximation methods yield better optimality, while ellipsoidal methods yield better computational efficiency. The distinction goes back to Sections 2.4 and 2.5, because polytope methods often rely on MIP programming while ellipsoidal methods tend to use SCP.
For the case where B(t) is a polytope, (Schouwenaars et al., 2001;Richards et al., 2002) first proposed to write (65e) as a set of mixed-integer constraints defined by the polytope facets. The resulting trajectory optimization can be solved using MIP methods discussed in Section 2.5. How, 2003a,b, 2006) apply this approach in the context of MPC with a variable horizon trajectory.
For the case where B(t) is an ellipsoid, (65e) is typically enforced by checking for collision using a conservative time-varying halfspace inclusion constraint: where H(t) is a halfspace. Three methods belonging to this family have been used. The first is a rotating hyperplane method, proposed by (Park et al., 2011;Di Cairano et al., 2012). Here, (65e) is replaced by a pre-determined sequence of halfspaces that are tangent to the ellipsoid and rotate around it at a fixed rate. This approach was first applied to a 2D mission, and later extended to 3D (Weiss et al., 2012(Weiss et al., , 2015a. A variation was introduced in (Park et al., 2016) and further studied in (Zagaris et al., 2018), where the rotating sequence is replaced by just two halfspaces tangent to the obstacle and passing through the chaser and target positions. This method requires to pre-specifying which of the two halfspaces the chaser belongs to at each time instant. Fixing the halfspace sequence enables the first two approaches to retain convexity. However, a third and most natural approach is to impose (66) directly by linearizing the ellipsoidal obstacle. This approach is taken in (Liu and Lu, 2014), and has also been applied to multiple moving obstacles (Jewison et al., 2015;. Because convexity is not maintained, SCP solution methods are used as discussed in Section 2.4. (Zagaris et al., 2018) provide a detailed comparison of the three methods.
Another challenge in RPO is the thrust constraint (65d). This constraint allows the thrust magnitude to take any value in the continuous interval [0, ρ]. In reality, however, control is often realized by a reaction control system (RCS) that produces short-duration pulses of constant thrust. Therefore, in many applications it makes more sense to consider an impulse constraint of the form: where ∆v(t) ∈ R approximates the instantaneous change in the chaser's velocity following a firing from the RCS jets. Realistic RCS thrusters have a minimum impulse-bit (MIB) performance metric that governs the magnitude of the smallest possible velocity change ∆v min > 0. Because (67) is a nonconvex disjoint constraint of the form (2), it has been historically challenging to handle. Indeed, (Larsson et al., 2006) suggest that MIP is necessary in general, but in certain cases the LP relaxation ∆v(t) ∈ [0, ∆v max ] of (67) suffices. This happens, for example, when the velocity measurement noise exceeds the MIB value. More recently, it was shown that the impulsive rendezvous problem can be solved via polynomial optimization (Arzelier et al., 2011(Arzelier et al., , 2013. Using results on non-negative polynomials,  showed that impulsive rendezvous with linear path constraints can be solved as an SDP. This formulation was further embedded in a glideslope guidance framework for RPO (Ariba et al., 2018) and in an MPC approach (Gilz et al., 2019). Distinct from polynomial optimization,  proved that in some special cases the constraint (67) can be losslessly convexified using techniques similar to those in Section 2.3. For problems where lossless convexification is not possible, (Malyuta and  showed that SCP with a numerical continuation scheme is an effective solution method. Yet another approach was presented in (Wan et al., 2019), where an alternating minimization algorithm was proposed for the case ∆v min = ∆v max , in other words when the control is bang-bang. The impulsive rendezvous model (67) considers an instantaneous firing duration. The model's accuracy can be improved by explicitly considering the finite firing duration, leading to a representation of the actual pulse-width modulated (PWM) thrust signal. PWM rendezvous was first studied in (Vazquez et al., 2011(Vazquez et al., , 2014, where an optimization similar to Problem 65 was first solved, then the optimal continuous-valued thrust signal was discretized using a PWM filter and iteratively improved using linearized dynamics. This approach was later embedded in MPC (Vazquez et al., 2015(Vazquez et al., , 2017. A subtly different approach is presented in (Li et al., 2016;Li and Zhu, 2018b), called pulse-width pulse-frequency modulation (PWPF). Instead of iteratively refining the thrust signal, PWPF passes the continuous-valued thrust signal to a Schmitt trigger that converts it into a PWM signal. It was shown that this can save fuel and that stability is maintained. However, a potential implementation disadvantage is that the duration of each period in the resulting PWM signal varies continuously, which conflicts with typical hardware where this period is an integer multiple of a fixed duration. An SCP approach was recently used to account for this via state-triggered constraints from Section 2.5.1 (Malyuta et al., 2020).
Although RPO literature tends to focus on the relative chaser-target position using a 3-DoF model, relative attitude control also plays an important role, especially if the target is tumbling Dong et al., 2020). Thanks to advances in the speed and reliability of optimization solvers as mentioned in Section 2.1, there has been an increasing interest to optimize 6-DoF RPO trajectories with explicit consideration of position-attitude coupling through constraints such as plume impingement and sensor pointing (Ventura et al., 2017;Zhou et al., 2019). The resulting 6-DoF RPO trajectory optimization, however, is much more challenging to solve due to the presence of attitude kinematics and dynamics. Nevertheless, a special case with field of view and glideslope constraints was presented in (Lee and Mesbahi, 2014), where 6-DoF RPO was solved as a convex quadratically constrained QP by using a dual quaternion representation of the dynamics, effectively establishing a convexification.
For more general RPO problems, nonlinear programming software has been used frequently. For example, (Ventura et al., 2017) used SNOPT (Gill et al., 2005) after parameterizing the desired trajectory using polynomials. A B-spline parameterization was used in (Sanchez et al., 2020), and the resulting nonlinear optimization was solved by the IPOPT software (Wächter and Biegler, 2005). MATLAB-based packages were also used in (Malladi et al., 2019;Volpe and Circi, 2019). Recently, SCP techniques discussed in Section 2.4 were applied to 6-DoF RPO trajectory optimization. (Zhou et al., 2019) considered both collision avoidance and sensor pointing constraints. (Malyuta et al., 2020) further considered integer constraints on the PWM pulse width in order to respect the RCS MIB value, and constraints on plume impingement, by using statetriggered constraints. The algorithm was improved in  by making the solution method faster and more robust. The approach uses homotopy ideas from Section 2.5.2 to blend the PTR sequential convex programming method with numerical continuation into a single iterative solution process.
The operation of two spacecraft in close proximity naturally makes RPO a safety-critical phase of any mission. Thus, trajectory optimization that is robust to modeling errors, disturbances, and measurement noise has been an active research topic. MPC has been a popular approach in this context, as it allows efficiently re-solving Problem 65 with online updated parameters using hardware with limited resources (Hartley and Maciejowski, 2014;Goodyear et al., 2015;Park et al., 2013). (Hartley, 2015) provides a tutorial and a detailed discussion. Among the many different approaches that have been developed to explicitly address robustness, we may count feedback corrections (Baldwin et al., 2013), the extended command governor (Petersen et al., 2014), worst-case analysis Xu et al., 2018), stochastic trajectory optimization (Jewison and Miller, 2018), chance constrained MPC (Gavilan et al., 2012;Zhu et al., 2018), sampling-based MPC (Mammarella et al., 2020), tube-based MPC (Mammarella et al., 2018;Dong et al., 2020), and reactive collision avoidance (Scharf et al., 2006). In addition to various uncertainties, anomalous system behavior such as guidance system shutdowns, thruster failures, and loss of sensing, also poses unique challenges in RPO. In order to ensure safety in the presence of these anomalies, (Luo et al., 2007c(Luo et al., ,b, 2008 used a safety performance index to discourage collision with the target, and (Breger and How, 2008) considered both passive and active collision avoidance constraints in online trajectory optimization. (Zhang et al., 2015a) considered passive safety constraints together with field of view and impulse constraints. Aside from optimization-based methods, artificial potential functions (Dong et al., 2017;Liu and Li, 2019) and sampling-based methods (Starek et al., 2017) have also been applied to achieve safety in RPO.

Small Body Landing
A maneuver similar to RPO is that of small body landing, where the target spacecraft is replaced by a small celestial object such as an asteroid or a comet. Trajectory optimization for small body landing has gathered increasing levels of attention, spurred by recent high-profile asteroid exploration missions including Hayabusa (Kawaguchi et al., 2008), Hayabusa2 (Crane, 2019), and OSIRIS-REx (Berry et al., 2013;Lauretta et al., 2017). Unlike planetary rocket landing from Section 3.1, small body landing dynamics are highly nonlinear due to the irregular shape, density, and rotation of the small body (Werner and Scheeres, 1997;Scheeres et al., 1998). Landing must furthermore ensure a small touchdown velocity, possible plume impingement requirements, and collision avoidance. These aspects pose unique challenges for trajectory optimization. This section reviews recent developments in convex optimization-based small body landing algorithms. Alternative trajectory optimization methods have also been studied for this problem and which we do not cover, such as indirect methods (Yang and Baoyin, 2015;Chen et al., 2019).
The prototypical small body landing OCP is illustrated in Figure 14 and can be summarized as follows: Note the similarity between Problems 63, 65, and 68. Compared to Problem 63, small body landing is expressed in the rotating frame of the target. Thus, the main difference is in the dynamics (68b) that contain a general nonlinear gravity term g(r(t)) and inertial forces from the non-negligible angular velocity ω of the small body. The glideslope constraint (68e) is also shared with the approach cone in RPO (65f).
Early work by (Carson III and Açıkmeşe, 2006;Carson III et al., 2008) ignored the mass dynamics (68c), while (68b) was linearized to solve for acceleration rather than a thrust profile. The resulting tube MPC algorithm includes a pre-determined feedback controller optimized using SDP and tracking a feedforward trajectory from an SOCP in a robust and recursively feasible manner. Some time later, (Pinson and Lu, 2015) solved for a fixed-duration trajectory by applying lossless convexification to (68d) and successive linearization to (68b), resulting in an SCP solution method consisting of a sequence of SOCP subproblems. (Pinson and Lu, 2018) further combine this solution procedure with Brent's line search method to solve for the minimum-fuel flight duration, which is similar to the use of golden-section search in the PDG context (Blackmore et al., 2010).  combined convexification with classic Runge-Kutta discretization to improve the solution accuracy. (Yang et al., 2017) showed how to solve the minimum-time landing problem as a sequence of convex optimization problems. As a byproduct, they showed that for timeoptimal and short-duration minimum-landing-error versions of Problem 68, the thrust stays at its maximum value, in which case the lower bound in (68d) can be removed and (68c) simplified. Constraint (68e) is the most basic type of collision avoidance constraint. The heuristic reasoning behind (68e) is that if the lander stays approximately above a minimum glideslope, then it will avoid nearby geologic hazards. An alternative two-phase trajectory optimization was introduced in Liao-McPherson et al., 2016) by splitting the landing maneuver into a circumnavigation and a landing phase. During circumnavigation, the spacecraft is far away from the landing site and (68e) is replaced by collision avoidance constraint with the small body. In the same manner as Section 3.2, the small body is wrapped in an ellipsoid and a rotating hyperplane constraint is used Liao-McPherson et al., 2016;Sanchez et al., 2018). (Reynolds and Mesbahi, 2017) introduced an optimal separating hyperplane constraint that also generates auxiliary setpoints for MPC tracking that converge to the landing site. Once in close proximity to the landing site, the spacecraft enters the landing phase where constraint (68e) is enforced to facilitate pinpoint landing.
Most small body landing work is 3-DoF in the sense that it considers point mass translational dynamics. However, recently (Zhang et al., 2020) studied a two-phase variable landing duration 6-DoF problem. The motivation was to impose a field of view constraint for a landing camera. The resulting nonconvex optimization trajectory problem was solved using SCP as covered in Section 2.4.
Parameters of the small body, such as ω and g, are often subject to inevitable uncertainty, requiring judicious trajectory design. As a result, many aforementioned works use MPC to cope with the uncertainty in small body landing (Reynolds and Mesbahi, 2017;Sanchez et al., 2018). Application examples include tube MPC (Carson III and Carson III et al., 2008) and input observers to compensate for gravity modeling errors Liao-McPherson et al., 2016). (Hu et al., 2016) also proposed to jointly minimize fuel and trajectory dispersion described by closed-loop linear covariance. For a detailed discussion on achieving robustness in small body landing, we refer interested readers to the recent survey (Simplício et al., 2018).

Constrained Reorientation
Scientific observation satellites commonly need to execute large angle reorientation maneuvers while ensuring that their sensitive instruments, such as cryogenically cooled infrared telescopes, are not exposed to direct sunlight or heat.
Famous examples include the Cassini spacecraft, the Hubble Space Telescope, and the upcoming James Webb Space Telescope (Singh et al., 1997;Long, 2004;Downes and Rose, 2001). This section discusses the challenges of constrained reorientation as a trajectory optimization problem, and focuses on how convex optimization methods have been leveraged to address these challenges.
A basic constrained reorientation OCP is illustrated in Figure 15 and can be formulated as follows: The set of constraints (69d) encodes conical keep-out zones for n stars, similarly to the illustration in Figure 15 for one star. The parameters M i ∈ R 4×4 are symmetric matrices that are not positive semidefinite, as introduced in Section 2.3. The main challenge of solving Problem 69 stems from the fact that (69d) and the attitude dynamics (69b)-(69c) are nonconvex. (Kim and Mesbahi, 2004) were the first to prove that (69d) can be losslessly replaced by convex quadratic constraints, provided q(t) 2 = 1. Based on this observation, (Kim and Mesbahi, 2004) proposed to greedily optimize one discretization point at a time instead of the entire trajectory jointly. The method was further extended to the case of integral and dynamic pointing constraints in (Kim et al., 2010).
Although the method of (Kim and Mesbahi, 2004) is computationally efficient, it is inherently conservative and may fail to find a feasible solution to Problem 69 by greedily optimizing one discretization point at a time. As a result, several attempts have been made to improve its performance. For example, (Tam and Lightsey, 2016) propose to replace constraint (69d) with penalty terms in the objective function in order to ensure that a feasible trajectory can be found. Binary logical variables were also introduced in (69d) to account for redundant sensors. (Hutao et al., 2011) showed how the convexification of constraints (69d) should be adjusted when optimizing an entire trajectory, rather than a single time step as originally done in (Kim and Mesbahi, 2004). Put into an MPC framework, the resulting trajectory optimization yields less conservative performance. Alternatively, (Eren et al., 2015a) proposed to first optimize a quaternion sequence without kinematic and dynamic constraints, and then to compute the corresponding angular velocity and torque using the resulting quaternions. A hyperplane approximation of the unit sphere is used during quaternion optimization to ensure dynamic feasibility, and is imposed via MIP. Recently, (McDonald et al., 2020) proposed an SCP method with a line search step that helps convergence, which provides a potential real-time solution to Problem 69. Figure 16: Illustration of a basic ascent scenario. The goal is to find an optimal angle-of-attack θ trajectory to transfer the launch vehicle's payload from the planet's surface to orbit, while minimizing fuel and satisfying structural integrity constraints.
Aside from the quaternion representation in Problem 69, which is the most popular choice, a direction cosine matrix representation of attitude was also used by (Walsh and Forbes, 2018) to solve an equivalent problem. The resulting trajectory optimization can be approximated as an SDP using successive linearization and relaxing (69d).
Due to its challenging nature, Problem 69 has inspired many optimization solutions other than those based on convex optimization. Pseudospectral methods and NLP optimization software have all been used to solve Problem 69 directly (Xiaojun et al., 2010;Lee and Mesbahi, 2013). An indirect shooting method was used in Phogat et al., 2018a), and a differential evolution method was used in (Wu et al., 2017;Wu and Han, 2019). Compared with convex optimization based methods, these methods typically require more computational resources to achieve real-time implementation.

Endo-atmospheric Flight
Launching from or returning to a planet with an atmosphere are integral parts of many space missions. These problems concern launch vehicles, missiles, and entry vehicles such as capsules, reusable launchers, and hypersonic gliders. Significant portions of launch and entry occur at high velocities and in the presence of an atmosphere, making aerodynamics play a large role. Aerodynamics and thermal heating are indeed the core differentiating factors between endoatmospheric flight and PDG from Section 3.1. For the latter problem, small velocities and thinness of the atmosphere make aerodynamic effects negligible in many cases (Eren et al., 2015b). This section summarizes recent contributions to endo-atmospheric trajectory planning using convex optimization-based methods. In particular, Section 3.5.1 discusses ascent and Section 3.5.2 discusses entry.

Ascent Flight
The optimal ascent problem seeks to transfer a launch vehicle's payload from a planet's surface to orbit while minimizing a quantity such as fuel. Naturally, optimal control theory from Section 2.1 has found frequent applications in ascent guidance, and we refer the reader to (Hanson et al., 1994) for a survey. Heritage algorithms date back to the iterative guidance mode (IGM) of Saturn rockets (Chandler and Smith, 1967;Horn et al., 1969;Adkins, 1970;Haeussermann, 1970) and the powered explicit guidance (PEG) of the Space Shuttle (McHenry et al., 1979). A simple yet relevant optimal control problem describing an orbital launch scenario is known as the Goddard rocket problem (Betts, 2010;Bryson Jr. and Ho, 1975). A version with variable gravity and no atmospheric drag is stated as follows: Problem 70 models a three-dimensional point mass moving in a gravity field under the influence of thrust. As such, it also applies to orbit transfer problems which we discuss later in Section 3.6. The vector r(t) ∈ R 3 is the position vector, T (t) ∈ R 3 is the thrust vector, and m(t) ∈ R is the vehicle mass. The vector function ψ : R 6 → R k imposes k ≤ 6 terminal conditions. In ascent and orbit transfer applications, ψ usually acts to constrain k orbital elements while leaving the other 6 − k orbital elements free.
A key issue when solving Problem 70 using an indirect method is to resolve the transversality conditions of the resulting TPBVP (Pontryagin et al., 1986;Berkovitz, 1974): where x(t) (r(t);ṙ(t)), p(·) ∈ R 6 are the costates relating to the position and velocity, and ν p ∈ R k is a Lagrange multiplier vector. Unfortunately, ν p has no physical or exploitable numerical interpretation, and the magnitudes of its elements can vary wildly (Pan et al., 2013). This causes a lot of difficulty for the solution process in terms of numerics, robustness, and initial guess selection. Traditionally, the problem is overcome by converting (71) into a set of 6 − k socalled reduced transversality conditions, which are equivalent (Lu et al., 2003): The linearly independent vectors y i ∈ R 6 are known as the reduced transversality vectors, and are a function of x(t f ). If they are known analytically, then (72b) can replace (71), which eliminates ν p from the problem and simplifies the solution process considerably. However, solving for y i symbolically is a difficult task, and the resulting expressions can be complicated (Brown and Johnson, 1967). An alternative approach was introduced in (Pan et al., 2013) where the authors provide an easy to use "menu" of the 6 − k constraints in (72b) that are needed. This is achieved by considering Problem 70 specifically and exploiting the structure offered by the classical orbital elements. It is only assumed that the terminal constraint function ψ fixes exactly k of the 6 orbital elements, and leaves the other orbital elements free.
The Goddard rocket problem in Problem 70 assumes no atmosphere. When there is an atmosphere, a popular classical method is the gravity turn maneuver, which maintains a low angle-of-attack so as to minimize lateral aerodynamic loads. However, the general ascent problem with an atmosphere is complicated due to strong coupling of aerodynamic and thrust forces . Thus, ascent is typically performed via open-loop implicit guidance, in the sense that feedback control is used to track a pre-computed ascent trajectory stored onboard as a lookup table. However, this approach cannot robustly handle offnominal conditions, aborts, and contingencies, which motivates research into closed-loop ascent techniques (Brown and Johnson, 1967;Lu, 2017).
A notable strategy in this context is to include aerodynamics in an onboard ascent solution via a homotopy method starting from an optimal vacuum ascent initial guess (Calise et al., 1998;Gath and Calise, 2001;Calise and Brandt, 2004). Another approach was developed in (Lu et al., 2003(Lu et al., , 2005Lu and Pan, 2010;Pan and Lu, 2010) using indirect trajectory optimization. Here, a finitedifference scheme is proposed to solve for the necessary conditions of optimality for ascent with an atmosphere. In particular, fixed-point formulations were considered (Lu et al., 2005), primer vector theory was invoked to determine trajectory optimality , and a generalization to arbitrary numbers of burn and coast arcs was developed . Finally, indirect methods relying on control smoothing via trigonometrization have been developed to address problems with bang-bang input and singular arcs (Mall et al., 2020). The Epsilon-Trig method (Mall and Grant, 2017), which is an example of such an approach, was applied to the Goddard maximum altitude ascent problem to obtain its bang-singular-bang optimal solution. See Section 2.1.1 for a brief description of these approaches.
Modern improvements in convex optimization have made direct optimization methods attractive for ascent guidance. To this end, consider the following illustrative ascent problem for a two-stage launch vehicle, as shown in Figure 16: Problem 73 is planar and formulated in an Earth-centered inertial (ECI) frame. Control is performed using the angle-of-attack θ, which determines the direction of an otherwise pre-determined thrust profile Li et al., 2020;Liu and Lu, 2014). The major aerodynamic forces are those of lift L and drag D, each of which may be complex expressions of state and control. Note that in (73b) we used the shorthand T [t], L[t], and D[t] from Section 1 to abstract away the possible state and control arguments. The atmospheric density is denoted by ρ, which varies during ascent as a nonlinear function of the position r. An example is given later in (75c). Important constraints on the dynamic pressure (73e) and bending moment (73f) are used to ensure the vehicle's structural integrity . The target orbit is prescribed by the vector function ψ in (73g), which is the same as in (70d) and specifies some or all of the target orbital elements. (Benedikter et al., 2020) chose boundary conditions based on the radius and inclination of a circular target orbit. A final nuance is that, if the rocket is assumed to be a two-stage vehicle, a stage separation event must be scheduled at a pre-determined time t 1 via (73h). At the separation instant, the mass variable experiences a discontinuous decrease that amounts to the dry weight of the first stage (Benedikter et al., 2019a(Benedikter et al., , 2020. A related constraint for stage separation requires θ(t 1 ) = 0 in order to reduce lateral load (Zhengxiang et al., 2018). Furthermore, the splashdown location of a burnt-out separated stage can also be constrained (Benedikter et al., 2020).
Due to the presence of strong nonlinearities, convex optimization-based solution algorithms for Problem 73 typically use SCP from Section 2.4. However, several manipulations have been helpful to make the problem less nonlinear. Conversion of the system dynamics (73b) to control-affine form, at times by choosing an independent variable other than time, followed by the use of lossless convexification within an SCP framework has been a common approach.  obtained a control-affine form by assuming the AoA to be small and defining u 1 = θ, u 2 = θ 2 as the new control variables. This choice makes drag a linear function of the control, while the constraint u 2 1 = u 2 is relaxed to u 2 1 ≤ u 2 via lossless convexification. Similarly, (Benedikter et al., 2019a(Benedikter et al., , 2020 chose thrust direction as input and losslessly convexified the unit norm constraint on the thrust direction to a convex inequality. (Cheng et al., 2017) considered a 3D problem with AoA and bank angle as control inputs, and applied lossless convexification to a constraint of the form u 2 1 + u 2 2 + u 2 3 = 1. Furthermore, their choice of altitude as the independent variable simplified the convexification of constraints involving density, since it renders the density a state-independent quantity. In particular, during collocation over a known grid within an altitude interval, the density value is known at each node. This fails to be the case when collocation is performed over time. The choice of altitude as independent variable was also explored in (Liu et al., 2016b).
The two-agent launch problem is an interesting and relevant modern-day extension of Problem 73 (Ross and D'Souza, 2004). In this case, the launch vehicle first stage is not just an idle dropped mass, but is a controlled vehicle that must be brought back to Earth. This is the case for the SpaceX Falcon 9 rocket, whose first stage is recovered by propulsive landing after a series of post- Figure 17: Illustration of a basic atmospheric entry scenario. The goal is to find an optimal bank angle σ trajectory to dissipate energy while meeting structural integrity and targeting requirements.
separation maneuvers (Blackmore, 2016). It was shown in (Ross and D'Souza, 2005) how hybrid optimal control can be used to solve the problem via mixedinteger programming. More generally, hybrid optimal control has also found applications in low-thrust orbit transfer using solar sails (Stevens and Ross, 2005;Stevens et al., 2004).

Atmospheric Entry
Atmospheric entry, also known as reentry, is fundamentally a process of controlled energy dissipation while meeting targeting and structural integrity constraints (Lu, 2014). Computer-controlled entry guidance dates back to the Gemini and Apollo projects, and (D'Souza and Sarigul-Klijn, 2014) provide a comprehensive survey of existing methods. Good documentation is available for Mars Science Laboratory's entry guidance, which is based on Apollo heritage (Way et al., 2007;Mendeck and Craig, 2011;Steltzner et al., 2014).
A large body of work is available on predictor-corrector methods for entry guidance Johnson et al., 2020Johnson et al., , 2018Johnson et al., , 2017Lu, 2014) and for aerocapture . These methods are based on root-finding algorithms, or variations thereof, and some versions are grounded in solving the necessary conditions of optimality (Lu, 2018). We refer the reader to (Lu, 2008) for further details. In addition to reentry trajectory generation, mission analysis tools for generating landing footprints have also been developed Eren et al., 2015b).
Guidance schemes based on univariate root-finding, which are near-optimal for reentry (Lu, 2014) and optimal for aerocapture , have also been developed. Reentry applications exploit the quasi-equilibrium glide condition (QEGC), while aerocapture leverages the bang-bang nature of the control solution obtained via the maximum principle. By using the knowledge of where the input switches, univariate root-finding can approximate the optimal solution in each phase to high accuracy. Such an approach, though based on an indirect method, avoids directly solving the TPBVP. Recall that lossless convexification, discussed in Section 2.3, is another approach where clever reformulation of the optimal control problem and application of the maximum principle yields an efficient solution strategy. This ties back to the last paragraph of Section 2.1.2, which states that the fusion of indirect and direct solution methods often yields more efficient solution algorithms than using any one method in isolation. Because root-finding algorithms do not involve an explicit call to an optimizer, we do not survey them here. Instead, this section focuses on contributions by convex optimization-based methods to the problem of entry trajectory computation.
Another methodology that simplifies the typical strategy in indirect methods is the RASHS approach (Saranathan and Grant, 2018). As discussed in Section 2.5.2, RASHS converts a multi-phase optimal control problem into a single-phase problem by using sigmoid functions of state-based conditions to instigate smooth transitions between phases. As a consequence, the multi-point BVP from Pontryagin's maximum principle is reduced to a TPBVP. The complete entry, descent, and landing (EDL) problem is one example that can be solved effectively via this technique.
Consider a basic entry guidance problem illustrated in Figure 17, which is formulated as follows: Problem 74 is planar and formulated in an ECI frame like Problem 73. Aerodynamic forces are governed by the lift, drag, and atmospheric density, which are expressed as follows: The lift and drag coefficients are given by c L and c D while ρ 0 , h 0 , and r 0 denote the reference density, reference altitude, and planet radius. R θ corresponds to a counter-clockwise rotation by θ radians. The control input u(t) = cos(σ(t)) is the cosine of the bank angle, and serves to modulate the projection of the lift vector onto the plane of descent, known also as the pitch plane. Entry is an extremely stressful maneuver for much of the spacecraft's hardware, therefore structural integrity constraints are placed on dynamic pressure (74d) and aerodynamic load (74e). The objective is to minimize the peak heating rate, given by the Detra-Kemp-Riddell stagnation point heating correlation (Detra et al., 1957;Garrett and Pitts, 1970), which is appropriate for an insulative reusable thermal protection system (TPS) such as on the Space Shuttle and SpaceX Starship:: Problem 74 is a simple example that gives a taste for the reentry problem. We now survey variants of this problem that have been explored in the literature. First of all, many other objectives have been proposed in place of (74a). These include a minimum heat load (Wang and Grant, 2017;Han et al., 2019), minimum peak normal load (Wang, 2019a,b), minimum time-of-flight Han et al., 2019), minimum terminal velocity (Wang and Grant, 2017), maximum terminal velocity (Wang and Grant, 2019), minimum phugoid oscillation , and minimum cross-range error (Fahroo et al., 2003a,b). In the problem of aerocapture, where a spacecraft uses the planet's atmosphere for insertion into a parking orbit, minimum velocity error (Zhang et al., 2015b) and minimum impulse, time-of-flight, or heat load  were studied. Minimizing the total heat load, which is equivalent to the average heating rate, is particularly relevant for ablative TPS that work by carrying heat away from the surface through mass loss. This has been the method of choice for Apollo, SpaceX Crew Dragon, and almost all interplanetary entry systems, because it can sustain very high transient peak heating rates (Hicks, 2009).
Problem 74 is expressed in the pitch plane and without regard for planetary rotation. To account for rotation and aspects like cross-range tracking, other formulations have been explored. This includes a pitch plane formulation with rotation (Chawla et al., 2010), a 3D formulation without rotation (Zhao and Song, 2017), and a 3D formulation with rotation (Wang and Grant, 2017;Wang et al., 2019a;Wang, 2019a;Wang and Grant, 2018a;Wang et al., 2019b;Liu et al., 2016aHan et al., 2019).
The two main path constraints present in Problem 74 are on the dynamic pressure (74d) and aerodynamic load (74e). The heating rate is also indirectly constrained since (74a) must achieve a lower value than the maximum heating rateQ max , otherwise the computed trajectory melts the spacecraft. Since these three constraints are critical for structural integrity, they permeate much of reentry optimization literature Grant, 2017, 2019;Wang et al., 2019a;Wang and Grant, 2018a;Wang, 2019a;Zhao and Song, 2017;Wang et al., 2019b;Liu et al., , 2016aHan et al., 2019;Sagliano and Mooij, 2018). Some researchers have also included no-fly zone (NFZ) constraints, as illustrated in Figure 17 (Zhao and Song, 2017;Liu et al., 2016a). A bank angle reversal constraint has also been considered in (Zhao and Song, 2017;Han et al., 2019;Liu et al., , 2016a. This is a nonconvex constraint of the form: A common approach to handle (77) is to define u 1 (t) cos(σ(t)) and u 2 (t) sin(σ(t)), and to impose: where the nonconvex equality constraint is subsequently losslessly convexified to u 1 (t) 2 + u 2 (t) 2 ≤ 1 Liu et al., 2016a. The bank angle with a prescribed AoA profile is a popular control input choice for reentry, dating back to Apollo (Rea and Putnam, 2007). Some works have considered bank angle rate as the input, which improves control smoothness Grant, 2017, 2019;Wang, 2019a;Wang et al., 2019a;Wang and Grant, 2018a;Wang et al., 2019b). However, banking is not the only possible control mechanism for reentry, and several other choices have been explored. (Chawla et al., 2010) use the AoA as input and omit bank and heading. (Fahroo et al., 2003a) use AoA, bank angle, and altitude, assuming the aforementioned QEGC with a small flight-path angle between the velocity vector and the local horizontal. (Zhao and Song, 2017) use bank angle and a normalized lift coefficient as inputs.
High frequency oscillation in the control signal, known as jitter, is a common issue in entry trajectory optimization. Several works explicitly address this issue (Szmuk et al., 2017;Wang and Grant, 2017;Liu et al., 2016a. Jitter is believed to be caused by the nonlinear coupling of state and control constraints (Wang and Grant, 2017), and it appears to be reduced by a controlaffine reformulation of the dynamics (Liu et al., 2016a). Other strategies to remove jitter have been to apply the reparametrization (78) or to filter the control signal. The latter approach includes the aforementioned use of bank angle rate as the control, or using a first-order low-pass filter .
Aside from fixing jitter, efforts have been devoted to simplifying the SCPbased solution methods, and to improving their convergence properties. Reformulating the dynamics using energy as the independent variable, in a similar way to how altitude was used for optimal ascent, is one tactic that achieves the former (Lu, 2014;Liu et al., 2016a). Such a parametrization eliminates the differential equation for airspeed, and instead yields an algebraic approximation for airspeed in terms of energy. (Fahroo et al., 2003a,b) applied a related elimination process by considering energy as state variable. Apart from this, it is worth noting the heuristics proposed for improving the convergence of the SCP-based approaches. Wang and Grant, 2019;Wang, 2019b) used backtracking line search at each SCP iteration to reduce constraint violation. It was found that with the line search, the number of iterations required for convergence reduced by half . (Zhang et al., 2015b) constrained the SCP iterates to form a Cauchy sequence.  proposed a dynamic trust region update scheme that is tailored for hypersonic reentry. In particular, the trust region update accounts for the linearization error due to each state instead of the typical approach of considering the average linearization error.
Aside from using SCP to optimize the entire entry trajectory, another popular approach for entry guidance is via MPC from Section 2.6. Some important recent MPC-based developments are the dynamic control allocation scheme  Figure 18: Illustration of a basic orbit transfer scenario. The goal is to use thrust T to transfer the spacecraft state from an initial orbit to a target orbit under the influence of gravity, while minimizing a quantity such as fuel or time. (Luo et al., 2007a) and the application of model predictive static programming (MPSP). The approach by (Luo et al., 2007a) centers around posing an SQP problem as a linear complimentarity problem, while the principle behind MPSP is to combine MPC and approximate dynamic programming though a parametric optimization formulation (Halbe et al., 2014(Halbe et al., , 2010Chawla et al., 2010). (van Soest et al., 2006;Recasens et al., 2005) corroborate the effectiveness of MPCbased approaches by comparing the performance of constrained MPC with that of PID control applied to feedback-linearized reentry flight. Last but not least, we conclude by discussing pseudosectral discretization from Section 2.2.4 as a popular methodology in a variety of reentry problem formulations. The method is appealing for its ability to yield accurate solutions with a relatively sparse temporal collocation grid, and recent results on the estimation of costates with spectral accuracy provide a strong theoretical grounding (Françolin et al., 2014;Gong et al., 2010). (Rea, 2003;Fahroo et al., 2003a,b) applied direct Legendre collocation to generate an entry vehicle footprint by solving a nonconvex NLP. (Sagliano and Mooij, 2018) used Legendre-Lobatto collocation and lossless convexification to generate an optimal profile for the heritage drag-energy guidance scheme. In addition to these approaches, which rely solely on the direct method, a combination of direct and indirect methods was discussed in (Josselyn and Ross, 2003) for verifying optimality of reentry trajectories using the DIDO solver (see Table 1). In (Tian and Zong, 2011), a feedback guidance law through an indirect Legendre pseudospectral method was developed to track a reference generated using a direct pseudospectral method. Finally, akin to explicit MPC, (Sagliano et al., 2017) developed a pre-computed interpolation-based multivariate pseudospectral technique that is coupled with a subspace selection algorithm to generate nearly optimal trajectories in real-time for entry scenarios in the presence of wide dispersions at the entry interface.

Orbit Transfer and Injection
A usual task in a space mission is to attain a certain orbit, or to change orbits. The goal of the so-called orbit transfer and injection (OTI) problem is to transport a low-thrust space vehicle from an initial to a target orbit while minimizing a quantity such as time or fuel. Unsurprisingly, the problem is as old as spaceflight itself, with the earliest bibliographic entry dating to the late 1950s (Faulders, 1958). Traditionally, the problem has been solved using optimal control theory from Section 2.1, and for this we can cite the books (Longuski et al., 2014;Lawden, 1963;Bryson Jr. and Ho, 1975;Kirk, 1970;Conway, 2014). Numerous solution methods have been studied, including methods based on primer vector theory (Russell, 2007;Petropoulos and Russell, 2008;Restrepo and Russell, 2017), direct methods based on solving an NLP (Betts, 2000;Arrieta-Camacho and Biegler, 2005;Ross et al., 2007;Starek and Kolmanovsky, 2012;Rao, 2015, 2016), and indirect methods (Alfano and Thorne, 1994;Fernandes, 1995;Kechichian, 1995;Haberkorn et al., 2004;Gong et al., 2008;Gil-Fernandez and Gomez-Tierno, 2010;Zimmer et al., 2010;Pan et al., 2012;Pontani and Conway, 2013;Cerf, 2016;Taheri et al., 2016Taheri et al., , 2017Rasotto et al., 2015;Lizia et al., 2014). Some recent advances for indirect methods include homotopy methods Pan and Pan, 2020;Cerf et al., 2011), optimal switching surfaces (Taheri and Junkins, 2019), the RASHS and CSC approaches from Section 2.5.2 (Saranathan and Grant, 2018;Taheri et al., 2020a,b), and simultaneous optimization (also known as co-optimization) of the trajectory and the spacecraft design parameters (Arya et al., 2020).
In a similar way to the previous sections, improvements in convex optimization technology have prompted an increased interest in applying the direct family of methods to OTI. For example, (Betts and Erb, 2003) solved a minimum-fuel Earth to Moon transfer using a solar electric propulsion system, which is a complex problem with a transfer duration of over 200 days. The problem is highly nonconvex, and the optimization algorithm is based on SQP. This section discusses some of the recent developments for solving OTI using convex optimization-based methods, and their extensions to optimal exo-atmospheric launch vehicle ascent.
A basic optimal OTI problem is illustrated in Figure 18 and can be formulated as follows: Just like in Problem 70, the vector function ψ in (79f) describes the final orbit insertion constraints, usually in the form of orbital elements. Note that (79b)-(79d) are identical to Problem 65. Naturally, we may hope that previously developed lossless convexification and SCP techniques from (Lu and Liu, 2013) apply for OTI. The main novelty is the nonlinear insertion constraint (79f). (Liu and Lu, 2014) showed that (79f) can be linearized with a second-order correction term, and Problem 79 can be solved via SCP as a sequence of SOCPs. The method is efficient and reliable, even for extremely sensitive cases like McCue's orbit transfer problem (McCue, 1967). Using similar convexification techniques, a 3D minimum-fuel OTI problem was considered in (Wang and Grant, 2018b). Similarly, a 2D minimum-time OTI problem was studied in (Wang and Grant, 2018c), where the dynamics were parametrized by transfer angle (i.e., orbit true anomaly) instead of time as the independent variable. Both works consider circular orbits, where (79f) can be linearized using spherical or polar coordinates. (Tang et al., 2018) solved a minimum-fuel orbit transfer problem by combining SCP with lossless convexification and pseudospectral discretization. (Song and Gong, 2019) studied a minimum-time interplanetary solar sail mission, where the thrust is replaced by solar radiation force, and optimized the trajectory via SCP as a sequence of SOCP problems.
The above paragraph mentions works that deal mainly with orbit transfer. A companion problem is that of orbit injection, where the vehicle is taken from a non-orbiting state to a target orbit. This occurs, for instance, in the last stage of rocket ascent. (Liu and Lu, 2014) showed that Problem 79 can also model the optimal exo-atmospheric ascent flight of a medium-lift launch vehicle. In this case, the initial condition (79e) typically denotes burnout of the launch vehicle's previous stage. In (Liu and Lu, 2014), constraint (79f) denotes the radius and velocity at the perigee of the target circular orbit. ) considered a similar optimal ascent problem where the thrust magnitude is constant, and constraint (79f) describes the orbital elements of a general elliptical orbit. Using pseudospectral discretization and SCP, this optimal ascent problem is solved as a sequence of SOCPs. (Li et al., 2020) further considered optimal ascent flight in the case of a power system fault. In this case, depending on the severity of the fault, (79f) describes progressively relaxed insertion constraints. constraints. Once the spacecraft is in orbit, an adjacent task is to avoid debris crossing its path. (Armellin, 2021) develops a real-time collision avoidance algorithm based on lossless convexification and SCP, and provides a detailed statistical analysis corroborating the method's effectiveness.
Mission planning often sits one layer above the OTI problem. For example, a mission plan may consist of a series of planetary flyby and gravity assist maneuvers. A mission, then, can be viewed as a sequence of OTI solutions that minimizes an overall objective such as fuel usage or travel time. A modern approach to mission planning is through hybrid optimal control, and some methods were already mentioned at the end of Section 3.5.1 (Ross and D'Souza, 2005;Stevens and Ross, 2005;Stevens et al., 2004). Evolutionary optimization using genetic algorithms offers an alternative solution for mission planning (Conway, 2014). This approach was used to plan several complex missions: a Galileo-type mission from Earth to Jupiter, a Cassini-type mission from Earth to Saturn, and an OSIRIS-REx type mission from Earth to the asteroid Bennu (Englander et al., 2012). The Saturn mission is almost identical to that used by the actual NASA/ESA Cassini mission, but is obtained fully automatically at a fraction of the time and cost. The algorithm, known as the evolutionary mission trajectory generator (EMTG), has been made available by NASA Goddard as an open-source software package (Englander, 2020).

Outlook
This paper surveyed promising convex optimization-based techniques for next generation space vehicle control systems. We touched on planetary rocket landing, small body landing, spacecraft rendezvous, attitude reorientation, orbit transfer, and endo-atmospheric flight including ascent and reentry. The discussion topics were chosen with a particular sensitivity towards computational efficiency and guaranteed functionality, which are questions of utmost importance for spaceflight control. We conclude by listing in Section 4.1 some of the most popular optimization software now available to the controls engineer, and outlining in Section 4.2 future research directions to which the reader may wish to contribute.

Optimization Software
Success in any computational engineering discipline owes in large part to the availability of good software. Table 1 lists modern optimization software packages that facilitate the implementation of algorithms discussed in Section 3. This list is by no means complete, and should be understood to merely indicate some of the popular optimization software packages that are quite mature and already available today.

Future Directions
We conclude this survey paper by listing some interesting and important future directions for optimization-based space vehicle control.

Guaranteed Performance
When proposing a new control algorithm for a real system, it is sobering to remember that the vehicle's survival, along with that of its occupants, literally hangs in the balance (Stein, 2003). The modern controls engineer has immense responsibility both to mission success and to upholding the foundation of trust created by the high reliability of traditional control methods. If we cannot guarantee an equal or greater level of reliability, then new optimization-based control methods will quite certainly be relegated to a ground support role .
Consequently, a direction of great importance for optimization-based space vehicle control is to rigorously certify that optimization-based algorithms converge to solutions that yield safe and robust operation in the real world. Active research is being done in the area, but general results are limited and many promising optimization-based methods lack proper guarantees. Today, Name Reference D ir ec t In di re ct R ea l-ti m e O pe n-so ur ce Ps eu dosp ec tr al C la ss B ac k-en d La ng ua ge
¶ Internal SQP-type methods as well as interface to external solver qpOASES. * * Written in C++ with a MATLAB interface available. ■ Written in Fortran 90 with a MATLAB interface available. SLLSQP, SNOPT, and others. ∆ Fortran 77, Ada 95, C, and MATLAB.
Conic operator splitting method (based on ADMM). Table 1: Summary of popular optimization software packages. The columns Direct and Indirect specify which solution method the software uses, as discussed in Section 2.1. The column Real-time denotes if the software is destined for real-time onboard use. Open-source identifies free-to-download packages with viewable source code. Pseudospectral identifies software that is compatible with pseudospectral discretization. Class describes the most general class of problems that the software can solve. Back-end lists which low-level optimizers are used, and Language lists the implementation and front-end interface languages of the package. Certain classifications that do not apply to the "generic parsers" software category are indicated by an empty cell background.
researchers are looking at real-time performance Malyuta et al., 2020), optimality , and convergence rates (Mao et al., 2018b;Bonalli et al., 2019). Perhaps the most important yet difficult guarantee is that the algorithm terminates in finite time, which is imperative for control. In the convex setting, algorithms with guaranteed convergence are available and have been flight-tested Scharf et al., 2017), so one direction to explore is how to convexify more general types of nonlinearity Liu and Lu, 2014;Lee and Mesbahi, 2014). For more difficult nonlinearities that are not convexifiable, an emerging subject of funnel libraries is being investigated (Majumdar and Tedrake, 2017;Reynolds, 2020;Açıkmeşe et al., 2008). The idea, akin to explicit MPC, is to pre-compute a lookup table of trajectories and invariant controllers in order to replace onboard optimization with a search algorithm followed by, in some cases, numerical integration. This can result in a substantially simpler onboard implementation at the expense of a higher storage memory footprint.

Machine Learning
Impressive advances in machine learning, and particularly in reinforcement learning (RL), could not side-step space vehicle control without due consideration (Tsiotras and Mesbahi, 2017). The main advantage of RL is that it is able to optimize over a stochastic data stream rather than assuming a particular description of a dynamic model (Buşoniu et al., 2018;Arulkumaran et al., 2017). As an optimization tool for nonlinear stochastic systems, it is not surprising that the RL method is attractive for aerospace control.
Although RL for space vehicle control is less than a decade old, a certain amount of literature is now available that addresses almost all of the applications presented in Section 3. The reader is referred to (Izzo et al., 2019) for a dedicated survey. In powered descent guidance, (Cheng et al., 2019a) use deep RL (DRL) for lunar landing,  improve ZEM/ZEV guidance via DRL, and (Gaudet and Linares, 2019) use recursive RL for Mars landing. For spacecraft rendezvous, (Scorsoglio et al., 2019) use actor-critic RL (ACRL) in near-rectilinear orbits, (Gaudet et al., 2018) consider cluttered environments, and (Doerr et al., 2020;Linares and Raquepas, 2018) use inverse RL to learn the target's behavior. In reentry guidance, (Shi and Wang, 2020) aim for real-time computation by training a deep neural network (DNN) to learn the functional relationship between state-action pairs obtained from a high-fidelity optimizer. Alternatively, (Cheng et al., 2020a) use a DNN to provide a numerical predictor-corrector guidance algorithm with a range prediction based on the current vehicle state. This method improves runtime performance by replacing traditional propagation-based trajectory prediction with a neural network. A different line of work is presented in (Jin et al., 2016), where the attitude of a reentry vehicle with model uncertainty and external disturbances is controlled by a robust adaptive fuzzy PID-type sliding mode controller designed using a radial basis function neural network. For small body landing, (Gaudet et al., 2020b use RL meta-learning for greater adaptability, and (Cheng et al., 2020c) train several DNNs to approximate a nonlinear gravity field as well as the optimal solution obtained using an indirect method. Another interesting approach was proposed in (Cheng et al., 2020b), where DNNs are used to supply good costate initial guesses, while an accurate trajectory is obtained by a downstream shooting method and a homotopy process. In orbit insertion and transfer applications, (Cheng et al., 2019b) develop a multiscale DNN architecture to approximate the optimal solution for a solar sail mission, (Holt et al., 2021) use ACRL for low-thrust trajectory optimization under changing dynamics, (LaFarge et al., 2020) use RL for libration point transfer in lunar applications, and  use proximal policy optimization.
A promising modern direction for spacecraft trajectory RL is to learn a small number of "behind the scenes" parameters (called solution hyperparameters) that govern the optimal solution, instead of directly learning the highdimensional optimal state-input map. Most importantly, the relationship between these parameters and the control policy is much more predictable, and hence can be learned more easily and with less training data. This survey paper makes it clear that most if not all spaceflight trajectory generation problems can be formulated as a variant of the optimal control Problem 4. Hence, the solution hyperparameters are often the maximum principle costates, or combinations thereof, that completely define the optimal control policy. Among these, we find aforementioned concepts of a primer vector (Açıkmeşe and Ploen, 2007;Lu and Pan, 2010;Lawden, 1963), and switching functions for bang-bang control (Taheri et al., 2020a). This RL approach was shown to be effective for 3-DoF PDG in (You et al., 2020a,b), where the authors learned 10 hyperparameters instead of the map from a 7D state to a 3D input. Most importantly, only ≈ 10 3 training trajectories were required. In comparison, the state-input map learning approach of (Sánchez-Sánchez and Izzo, 2018) also achieved good results, but required ≈ 10 7 training samples. A slightly different approach was taken for 3-DoF small body landing in (Cheng et al., 2020b), where homotopy and coordinate transforms were used to learn a 5D costate vector instead of the map from a 7D state to a 3D input. The DNN's output was then used to provide accurate initial guesses and to improve the convergence of a downstream shooting method. To summarize, the fact that learning hyperparameters works better than learning the optimal state-input mapping is just an observation of the fact that application domain knowledge can go a long way towards improving learning performance (Tabuada and Fraile, 2020). In the case of spacecraft trajectory optimization, this knowledge often comes from applying Pontryagin's maximum principle.
As discussed in Section 4.2.1, performance guarantees for an RL-based controller will have to be provided before serious onboard consideration. This may be harder to achieve for RL, since controllers are typically based on neural networks whose out-of-sample performance is still very difficult to characterize. Nevertheless, RL and other machine learning approaches are appealing for adaptive control systems. Future research will likely see the aerospace control community search for the right opportunities where RL can be embedded to improve traditional control systems.
to Taylor P. Reynolds for his invaluable review of an early draft and for his knowledge of space vehicle control at large. Arzelier, D., Louembet, C., Rondepierre, A., Kara-Zaitri, M., 2013 Carson III, J., Açıkmeşe, B., 2006. A model-predictive control technique with