# Can someone give me a 'trajectory optimisation for dummies' answer?

So, I am currently doing a project on trajectory optimisation. I get the principle of it, minimise fuel cost for a journey, or try to maximise the spacecraft weight at target destination, and that spits out an optimal trajectory to do so. I'm just not sure.. how to do it really? I'm reading a lot on things like 'primer vectors' and shape based methods and even Hamiltonians (which scare the **** out of this poor engineer) and I'm just not sure how to even start on the construction of an optimal trajectory?

Is there anyway to explain this for dummies, or a book to explain it for dummies?

• While there are currently four up votes, there are also three votes to close and yet no helpful comments how to improve the question, so I'll try to help. The two votes to close are for the reason that the question is "too broad", and that is probably true. It doesn't mean that it is a bad question, but that it is not a good match for Stack Exchange. I'd recommend you try to narrow this down to a more specific question, something that can easily be answered in a few sentences to a few paragraphs at most. Right now, the entire hypothetical book "trajectory optimisation for dummies" – uhoh May 4 '18 at 4:36
• is really needed to answer this question and that's way too much for one question. You could ask one question about prime vectors for example; quote a passage that uses the term and give a link to the source, and ask a question specific to it. You could also look at the dozens of questions and answers here about "trajectory optimization" (note the alternate spelling) space.stackexchange.com/search?q=trajectory+optimization and then ask about a specific topic that isn't fully explained sufficiently for you. – uhoh May 4 '18 at 4:36
• cheers! Appreciate the help, very new to this site! will try be more narrow in future! – Harvey Rael May 4 '18 at 14:08
• I have found this site to be the most useful if I can break it into individual questions. Put the simplest version of the question into the question title, and leave clarification for the question body. – Dan Sorensen May 19 '18 at 16:33
• @uhoh It seems like a book recommendation would, in fact, be useful. Do you have such a book that you know of? – Doctor Mohawk Jul 31 '18 at 14:38

I'm in search of a similar answer, so I thought I'd start one here and let others offer input. Here's what I've found so far.

I'll start to compile some general links here.

# Trajectory design steps

## Bootstrap

To plan a mission, you need to identify some general mission constraints. If it's a new problem, you might need to bootstrap a solution in order to understand what constraints might exist.

To do that, you might pick a specific mission start date when you think a trajectory might be feasible, as well as some other specific constraints (some amount of available fuel, engine Isp, etc). According to Ocampo (Chapter 4 of Spacecraft Trajectory Optimization, edited by Conway), who created the NASA Copernicus trajectory optimization program, there are four basic steps:

1. Find a feasible impulsive solution (instantaneous delta-V).
2. Optimize the impulsive solution (e.g. by minimizing delta-V).
3. Convert impulsive maneuvers to finite burns (non-instantaneous burns) using suboptimal control and optimize again (this time by maximizing final mass or by minimizing initial wet mass required).
4. Reconstruct thrust vectors using optimal control and optimize a final time.

I think steps 1-2 are all you really need for generating C3 or DV plots. After you identify your general constraints (below), you can pick a specific set of mission dates and do all of the steps.

If reference trajectories already exist, you can skip steps 1-2, generate your C3/DV plots, and then do all four steps once you've determined your departure and arrival dates and such.

It's worth noting that Paul has provided some really good information on how the optimization itself is done in another answer for this question. I won't try to improve upon his answer right now.

## Identify general constraints

Where are you going? What are the minimum and maximum mission durations? When is the earliest you could leave? The latest? Are you interested in low-thrust trajectories or something more Hohmann-like?

### Interplanetary

Use a porkchop plot for extremely approximate trajectories that meet your selected bounds.

For Earth–Mars, there's a web-based tool for generating these. This tool seems to use Lambert targeting, so it does direct transfers, and probably isn't going to help much with low-energy trajectories involving gravity assists.

More generally, NASA has a trajectory search engine (which provides a lot of information about this topic). As with the above app, and apparently in general, "The database holds solutions for direct transfers and simple gravity assists on outbound legs. More complex gravity assist maneuvers and other strategies may be used to find alternative trajectories not listed in the database."

It's not uncommon for porkchop plots to include big empty spaces where no trajectory information is shown. For example, consider this one from NASA's Trajectory Browser:

It may be that the authors didn't compute trajectories in the white area because they knew the DV would be too high, or it may be that they just aren't displayed because they're orders of magnitude higher than those shown and skew the colorbar. The take-away here is that the porkchop plot should only show relevant trajectories.

### Lunar

I haven't seen any porkchop plots for lunar missions (so far). However, the idea seems to be the same: you identify your constraints and choose a class of missions (e.g. direct transfers) to look at. You then consider either the specific energy (C3) or the DV.

For example, consider this chapter from JPL's Low-Energy Lunar Transfer Design book. They start by looking at three different classes of missions: direct transfer, short low-energy, long low-energy. They vary the longitude of the ascending node (which is a function of departure date) and show how the delta-V for lunar insertion changes for each of these.

So, in general, they've figured out constraints for their missions and then determined which variables exert the most control over DV (or C3).

## Develop informed constraints

Once you've produced your porkchop (or delta-V) plots, you can identify the constraints which simultaneously meet your mission requirements and provide the most flexibility (e.g. long launch window, most tolerance for failures, etc.).

And I'll stop here for now and come back and add more later (as I learn).

In its broadest form, Trajectory Optimization is simply finding the best path (and control values to produce this path) that satisfy the dynamics. The dynamics, in this case, are either the orbital dynamics (e.g. kepler's 2 body problem with proper ephemeris), rocket ascent dynamics or ballistic/lifting reentry dynamics. The cost function can be created in many different ways, depending on what you want in the end. For example, you may want to optimize the dynamics to produce a particular state at the end of a fixed time interval. Or you may want to minimize the fuel used in the trip (which maximizes the payload mass). Or you may want to minimize the time it takes to reach a particular end state using the least amount of (throttleable/vectored thrust).

Once you establish what your dynamic equations are and what your "cost function" are, you can translate that into a set of equations which you can solve. This process is called "transcription" in trajectory optimization nomenclature. Broadly speaking, there are two ways to transcribe the problem. The "indirect" method suggests that you should formulate the optimality conditions first analytically, then discretize those conditions over some number of timesteps N. This will result in a non-linear system of equations with $M\cdot N$ variables, where M is the number of variables to describe a state at a given time (e.g. the number of differential equations); plus $P\cdot N$ control variables, where P is the number of controls that are applied in each timestep. For example, in an 2-body problem in cartestian coordinates, you will have four state variables at each timestep (x-position, x-velocity, y-position, and y-velocity). The control may be a throttleable thrust which always points in the velocity direction. So, there would be only 1 control at each timestep. Often times, this can be tedious to do since a lot of the work has to be done by hand.

The alternative is referred to as the "direct" method. In this case, you would discretize the equations first and then formulate the conditions for optimality. So instead of formulating the conditions on an infinite space (indirect method), you're formulating them only for the finite dimensional space (a finite number of state + control variables). The advantage of this method is that the problem now becomes a non-linear constrained optimization problem (which is easier to program directly and solve using software packages).

Under the umbrella of direct methods, there are couple of different ways to discretize the problem. One way is called collocation, where we force the dynamics to agree at certain points in time (while assuming some sort of interpolation scheme between these points, such as linear or quadratic interpolation). The other approach is called shooting, which is kind of like a trial and error approach where the discrepancy (a.k.a. the defect constraint) at each timestep is discretized and also minimized. I haven't used shooting methods for trajectory optimization, but I have used collocation method so I will address it directly.

The key trick to the collocation method is knowing how to compute the non-linear equality constraints using the dynamics of the problem. The key idea is that on each time step, one must constrain the predicted value from the dynamics to be equal to the true value at the end of the timestep. Suppose we are considering the ith timestep, which has state variables $x_i$ and $x_{i+1}$ respectively at timesteps $t_i$ and $t_{i+1}$. Let's also suppose that the dynamic equation can be written in the form $\dot{x}=f$. If we integrate both sides of this equation with respect to time over the interval from $t_i$ to $t_{i+1}$, we get $$x_{i+1}-{x_i}=\int_{t_i}^{t_{i+1}}f dt.$$

One can discretize the integral in many ways. For simplicity, let's choose something like the trapezoid rule. Then we can obtain $$x_{i+1}={x_i}+\frac{f_{i}+f_{i+1}}{2} \Delta t.$$

Let's define the value predicted by the dynamics as $x_{i+1}^p={x_i}+\frac{f_{i}+f_{i+1}}{2} \Delta t.$ The actual value at the end of this timestep is $x_{i+1}$. So, to impose a dynamic constraint, we simply set the predicted value equal to the true value at the end of this time interval. Mathematically, we write: $x_{i+1}^p-x_{i+1} =0$. We write it this way because most nonlinear programming packages require you to write equality constraints with the RHS equal to zero.

That's the majority of the hard work. You simply need to create a function which computes a vector containing all of these constraints at each timestep simultaneously. Similarly, you will have to write your objective function in terms of the discrete variables $x_i$, which may include a discretized integral over time and/or a function of the end state.

Once you have the objective function and the equality constraint function (both of which are functions of the state variables and controls at all timesteps), you will also have to provide lower and upper bounds on each of the state and control variables at all timesteps. These bounds are usually physically motivated, like mass must be positive or the orbit radius cannot exceed 2 astronomical units, etc... Additionally, you often have to provide an initial guess that satisfies the lower and upper bounds. Sometimes, a trajectory optimization problem may be sensitive to your initial guess so it is a good idea to choose it as close to what you think the true solution should be as possible (e.g. if you have any physical intuition into the problem, now is the time to use it).

Once you have all this information, you can just feed it to a nonlinear programming solver (e.g. Matlab's fmincon, ipopt, etc...). There are usually different methods that you can use to solve nonlinear constrained optimization problems, such as sequential quadratic programming (SQP), active set methods (ASM), or interior point methods (IP). Of the methods I've used in trajectory optimization for orbital mechanics, I have found interior point methods to be the fastest and most robust. Regardless of which non-linear programming method you choose, make sure that when you check for convergence, check both feasibility (how much constraint violation are we incurring) and first-order optimality (how close are we to being at the bottom of a bowl/valley). Both values need to be small in order to be considered a "converged" solution.

Before you start programming a trajectory optimization code for orbital mechanics, I highly recommend trying a simpler problem with a known solution. The so-called bang-bang control problem is a good example of a 1D optimal control problem involving 1D motion. IN this problem, one seeks the optimal acceleration required to reach a certain distance in the minimum time such that your final velocity is also zero. This is a classic problem with lots of physical intuition (and even an analytical solution). Intuitively, the best control would be to accelerate as fast as you can toward your goal destination up until a certain point and then decelerate as fast as possible until you just barely reach your destination point with a velocity of zero.

You can see how this bang-bang problem can be transcribed using the direct collocation method (complete with matlab code) on Sam Pfrommer's blog:

For additional information on trajectory optimization in general, with more complicated examples such as the cart-pole control problem, I suggest reading/watching Matthew Kelley's tutorial on this topic:

[ This is going to be a draft answer as I learn more and get more time to fill it out and edit with symbolic math -- right now its going to be super rough and this answer will evolve a lot over time ]

Basic Simplified Problem

The generic problem is to take a spacecraft from a set of initial conditions $$\vec{r}(0), \vec{v}(0), m(0)$$ to a set of terminal conditions $$\vec{r}(t_f), \vec{v}(t_f)$$ where $$\vec{u}(t)$$ is a continuous control to be solved for which is the direction of thrust of the rocket, subject to a cost function $$J$$ to be minimized. The cost function is commonly minimum time, or maximum terminal mass for a free end time problem, or maximum orbital energy for a fixed-time problem. There may be other discrete controls like switching times or other continuous controls if throttling is involved, along with path constraints and other complexity which I will ignore for now.

If one is given the solution $$\vec{u}(t)$$ for all t, then the problem will reduce to a simple initial value problem (IVP) with a coupled system of ordinary different equations (ODEs), which could be integrated forward with an ODE solver like Runge–Kutta–Fehlberg (RKF45) method or the Dormand–Prince (RKDP) solver in Matlab's ode45.

Since we are not given $$\vec{u}(t)$$ and must instead solve for that control such that the terminal boundary conditions are met it becomes a two point boundary value problem (TPBVP).

Indirect Method

The "indirect" method uses the Calculus of Variations (CoV) and a branch of mathematics called Optimal Control Theory which utilizes Pontryagin's maximum (or minimum) principle. When applied to the problem of optimizing rocket trajectories, this becomes Lawden's Primer Vector Theory, where much of the mathematical framework has been done for the engineer.

The indirect method typically uses shooting (or multiple shooting) methods to integrate a guess at the initial conditions forward in time to the terminal state. The residual of all the boundary values are constructed and then a root finding tool like Levenberg-Marquardt or a quasi-Newton trust region method are used (vanilla Newton-Raphson will have convergence difficulties with this nonlinear problem).

Omitting the mathematical development, the problem for a rocket in a vacuum, subject to an inverse square central force potential $$g(r) = -\frac{1}{r^3}\vec{r}$$ where the final mass is to be maximized the cost function is simple Mayer problem:

$$J=φ(\vec{r}_f, \vec{V}_f, m_f, t_f) = −m_f$$

[EDIT: need a commment about normalization] The Hamiltonian I will state as:

$$H=\vec{P}^T_r \vec{V} + \vec{P}^T_V[−\frac{1}{r^3} \vec{r} + T \vec{1}_b ] + P_m (−\sqrt{\frac{R_0}{G_0}} \frac{T \cdot m}{I_{sp}} )$$

[...still work-in-progress here...]

The indirect method can be solved using Matlab's fsolve or lsqnonlin functions for root finding, and the ode45 solver for integrating IVPs and shooting. Or the collocation solvers bvp4c/bvp5c can be used to integrate the state and costate and to apply interior constraints.

The limitations of the indirect method are:

1. Guessing the costate becomes very hard and more of an artform as the problem complexity increases.
2. The solution often becomes sensitive to small perturbations in the initial costate and does not converge (particularly with stiff problems like an exponential atmosphere).
3. Path constraints (such as dynamic pressure or angle of attack constraints) become difficult to employ.
4. The complexity of the primer vector equation becomes difficult when considering atmospheric flight and path constraints.
5. The transversality equations are difficult to derive and are an artform.

The advantages of the indirect method are:

1. Fast, particularly when using step-size adjusting integrators or analytic computation of the trajectory
2. Relatively simple for unconstrained solutions in a vacuum requiring only an ODE solver and a root finding algorithm

Direct Method

The direct method uses some for of orthogonal collocation and discretization to convert the TPBVP into a nonlinear programming problem. This eliminates the costate and the costate equations and the transversality conditions from the problem. There will generally not be one unique solution to the problem which satisfies the terminal conditions, but the cost function $$J(t)$$ is directly applied to find the solution which satisfies the terminal constraints.

[... huge number of details omitted which I'm still only learning about...]

The Matlab function fmincon or the Matlab bindings to ipopt can be used to solve this kind of problem, SNOPT is a commonly used commercial package used to solve this kind of problem which also has Matlab bindings.

The advantages of the direct method include:

1. No guessing of the initial costate values, much less sensitivity to the initial control guess
2. No costate equations
3. No transversality conditions
4. Much easier to apply interior path constraints directly to the problem