M9-T4 — Powered-Descent Guidance by Convexification

Learning objectives

After this lesson, you should be able to:

• State the fuel-optimal powered descent problem with thrust floor, ceiling and glide slope.
• Introduce the slack variable and show that the relaxed problem is a second-order cone program.
• Show from the minimum principle that the relaxation is lossless.
• Solve the problem for Hopper and read the bang-bang structure of the thrust.

The intuition

UAV M5-T8: Model Predictive Control solves constrained problems with a convex cost; the lander problem adds a norm lower bound that the earlier lessons treat as a defect. The explicit guidance of M8-T3: Explicit Guidance to a Terminal State found the least-energy descent with no limits, and the minimum principle of M9-T1: The Optimal Control Problem and the Minimum Principle the bang-bang structure of a bounded problem. This lesson puts the two together for a vehicle that has to hold a thrust between a floor and a ceiling.

The problem

Hopper has a rocket engine that, once lit, cannot throttle below a floor: here 720 N, 30% of the maximum 2400 N. It must go from 250 m up and 150 m uprange, moving at (10, −25) m/s, to the pad at rest. It must stay inside a glide-slope cone above the pad, a wedge that keeps it from sweeping low over the ground short of the pad, and it must use as little propellant as it can. Propellant flow is proportional to thrust, so the cost is the integral of the thrust magnitude, $\int\|\mathbf{T}\|dt$.

Where it is not convex

Dynamics (constant gravity), the glide slope, the terminal conditions and the thrust ceiling are all linear or conic in the thrust vector. The floor, $\|\mathbf{T}\| \ge T_{\min}$, excludes the disk of small thrusts: the feasible set is a ring, as in the figure of the previous lesson.

The relaxation

Introduce a slack $\Gamma$ for the magnitude and write $\|\mathbf{u}\| \le \Gamma$ with $\Gamma_{\min} \le \Gamma \le \Gamma_{\max}$. That is convex: a cone and two linear bounds. It is also a relaxation: it admits controls with $\|\mathbf{u}\| < \Gamma$, in particular thrusts below the floor that the original problem does not allow. The point of lossless convexification is that the optimum does not use that freedom. Because the cost is $\int\Gamma\,dt$, there is no reason to take $\Gamma$ larger than needed, and the minimum principle shows that at every instant $\|\mathbf{u}\| = \Gamma$: the relaxed optimum is feasible for the original problem and therefore optimal for it.

What the optimum looks like

For Hopper from 250 m, the solutions for different descent times show that fuel is not monotonic in time: 126.7 m/s at 8 s, 120.1 at 9 s, 126.9 at 10 s, 144.5 at 12 s, 182.6 at 16 s, and the problem is infeasible at 7 s. Short descents need more thrust to stop; long ones lose to gravity. At the best time of 9 s, the thrust sits at the floor of 720 N for about 4.5 s (12 steps; the vehicle falls almost freely with a little thrust) and then at the ceiling of 2400 N for the last 4 s to stop. This is the bang-bang structure of the minimum principle, now found by a solver with the glide slope included.

Replace the thrust norm by a slack bounded between the floor and the ceiling: the problem is a convex cone program. The optimum has norm equal to slack at every instant, so the relaxation costs nothing, and the answer is bang-bang between the floor and the ceiling.

Definitions

Thrust acceleration $\mathbf{u}$

$\mathbf{u} = \mathbf{T}/m$: thrust per unit mass; the acceleration of the vehicle is $\mathbf{u} + \mathbf{g}$.

Slack $\Gamma$

A variable bounding the thrust acceleration magnitude from above, with bounds $\Gamma_{\min} = T_{\min}/m$ and $\Gamma_{\max} = T_{\max}/m$.

Glide-slope constraint

$z \ge \tan\gamma\,\|(x, y)\|$: the position stays inside a cone above the target, with half-angle measured from the ground.

Lossless convexification

A relaxation of a nonconvex constraint whose optimal solutions satisfy the original constraint.

Mass profile

The variation of mass along the descent that converts the thrust bounds to acceleration bounds; here a fixed nominal profile.

Derivation: why the relaxation is lossless

Apply the minimum principle to the relaxed problem. With position $\mathbf{r}$, velocity $\mathbf{v}$ and the costates $\boldsymbol\lambda_r, \boldsymbol\lambda_v$, the Hamiltonian is

1. $$H = \Gamma + \boldsymbol\lambda_r^{\mathsf{T}}\mathbf{v} + \boldsymbol\lambda_v^{\mathsf{T}}(\mathbf{u} + \mathbf{g})$$

   The running cost is the slack (the fuel); the dynamics are the double integrator with thrust acceleration u. The glide-slope constraint, when inactive, does not enter.

2. $$\dot{\boldsymbol\lambda}_r = \mathbf{0},\qquad \dot{\boldsymbol\lambda}_v = -\boldsymbol\lambda_r$$

   The costate equations: $\boldsymbol\lambda_r$ is constant, so $\boldsymbol\lambda_v$ is a vector that is linear in time: $\boldsymbol\lambda_v(t) = \mathbf{a} - \mathbf{b}t$.

3. $$\min_{(\mathbf{u}, \Gamma)}\ \Gamma + \boldsymbol\lambda_v^{\mathsf{T}}\mathbf{u}\quad\text{s.t.}\ \|\mathbf{u}\| \le \Gamma,\ \Gamma_{\min}\le\Gamma\le\Gamma_{\max}$$

   At each instant the optimal (u, Γ) minimizes the Hamiltonian over the allowed set: the part of H that depends on the controls.

4. $$\min_{\|\mathbf{u}\|\le\Gamma}\boldsymbol\lambda_v^{\mathsf{T}}\mathbf{u} = -\|\boldsymbol\lambda_v\|\,\Gamma\quad\text{at } \mathbf{u} = -\Gamma\,\boldsymbol\lambda_v/\|\boldsymbol\lambda_v\|$$

   For a given Γ the best direction is opposite to the costate and the best magnitude is the largest allowed: the norm of u equals Γ, whenever the costate is nonzero.

5. $$\min_{\Gamma}\ \Gamma\,(1 - \|\boldsymbol\lambda_v\|) \;\Rightarrow\; \Gamma = \Gamma_{\min}\ \text{if } \|\boldsymbol\lambda_v\| < 1,\quad \Gamma = \Gamma_{\max}\ \text{if } \|\boldsymbol\lambda_v\| > 1$$

   Substitute the best u and minimize over Γ: a linear function of Γ with coefficient $1 - \|\boldsymbol\lambda_v\|$; the minimum is at a bound.

6. $$\|\boldsymbol\lambda_v(t)\| = 1 \text{ only at isolated instants}$$

   The norm of a vector that is a straight line in time equals 1 on an interval only if the line is constant with norm 1, which a nontrivial (normal) problem excludes: for a generic problem it crosses 1 at most at two points, which are the switch times.

The Hamiltonian is linear in the slack, and the costate is a straight line in time, so the optimal thrust magnitude is at the floor or the ceiling and switches where the costate has length one. The optimum satisfies the thrust floor, so the relaxation loses nothing.

A worked example

Hopper: start (−150, 250) m at (10, −25) m/s, 35° glide slope, thrust 720 to 2400 N, 24 steps, Isp 150 s, propellant flow proportional to thrust. The thrust bounds are converted to acceleration bounds with a fixed nominal mass profile (mean thrust 1400 N).

1. Feasibility and fuel against time. The problem is infeasible at 7 s. The fuel proxy $\int\Gamma\,dt$ is 126.67 m/s at 8 s, 120.11 at 9 s, 126.88 at 10 s, 144.51 at 12 s and 182.63 at 16 s: a minimum at about 9 s. Each solve took 10 barrier stages.
2. Losslessness. At 9 s the largest of $\Gamma_k - \|\mathbf{u}_k\|$ over the steps is 5×10⁻¹⁰ m/s²: the cone constraint is active at every step, so the relaxed solution is a true solution of the nonconvex problem.
3. Structure. The thrust is 1785 N at the first step, then 720 N (the floor) for 12 steps, then 2288 N and 2400 N (the ceiling) for the last 11 steps. The thrust angle from the vertical goes from +39.5° at the start through 0.8° in the middle to −25.4° at the end.
4. Glide slope. The path stays above the cone and touches it only at the pad: the smallest margin at the nodes is zero.
5. The mass profile is an approximation. With the true mass history (thrust $m\Gamma$ with the actual mass after burning), the propellant used is 9.43 kg against 8.56 kg assumed, the thrust exceeds the 2400 N ceiling by at most 39.8 N (1.7%) and never falls below the floor. The residual is what a log-mass formulation, or re-solving in flight, removes.

import math
import gnckit as g

G0, ISP = 9.81, 150.0
alpha = 1.0 / (ISP * G0)                       # propellant flow per newton of thrust, kg / (N s)
M0, TMAX, TMIN = 120.0, 2400.0, 720.0
r0, v0 = (-150.0, 250.0), (10.0, -25.0)         # downrange x and height z, metres and metres per second
GAMMA = math.radians(35.0)                      # glide-slope cone half angle from the vertical complement: z >= tan(35 deg) |x|
N = 24
mdot_nom = alpha * 1400.0                      # nominal mean thrust used to fix the mass profile

def solve(T, verbose=False):
    h = T / N
    nv = 3 * N                                  # ux, uz, Gamma at each step
    mk = [M0 - mdot_nom * (k + 0.5) * h for k in range(N)]
    def pos_row(ax, k):                        # position of axis ax at node k as a linear function of the controls
        row = [0.0] * nv
        for j in range(k):
            row[ax * N + j] = h * ((k - j) * h - 0.5 * h)
        return row
    def vel_row(ax, k):
        row = [0.0] * nv
        for j in range(k):
            row[ax * N + j] = h
        return row
    gvec = (0.0, -G0)
    Aeq, beq = [], []
    for ax in (0, 1):
        Aeq.append(vel_row(ax, N)); beq.append(0.0 - v0[ax] - gvec[ax] * T)
        Aeq.append(pos_row(ax, N)); beq.append(0.0 - r0[ax] - v0[ax] * T - 0.5 * gvec[ax] * T * T)
    G, hh = [], []
    for k in range(1, N):                       # glide slope at the nodes: z - tan x >= 0 and z + tan x >= 0
        tk = k * h
        px, pz = pos_row(0, k), pos_row(1, k)
        cx = r0[0] + v0[0] * tk
        cz = r0[1] + v0[1] * tk - 0.5 * G0 * tk * tk
        for sgn in (1.0, -1.0):
            row = [-(pz[i] - sgn * math.tan(GAMMA) * px[i]) for i in range(nv)]
            G.append(row); hh.append(cz - sgn * math.tan(GAMMA) * cx)
    for k in range(N):                          # Gamma_min(k) <= Gamma_k <= Gamma_max(k)
        row = [0.0] * nv; row[2 * N + k] = 1.0
        G.append(row); hh.append(TMAX / mk[k])
        row = [0.0] * nv; row[2 * N + k] = -1.0
        G.append(row); hh.append(-TMIN / mk[k])
    socs = []
    for k in range(N):
        F = [[0.0] * nv, [0.0] * nv]; F[0][k] = 1.0; F[1][N + k] = 1.0
        q = [0.0] * nv; q[2 * N + k] = 1.0
        socs.append((F, [0.0, 0.0], q, 0.0))
    c = [0.0] * (2 * N) + [h] * N
    x, info = g.solve_convex(c, A=Aeq, b=beq, G=G, h=hh, socs=socs)
    return x, info, mk, h

def replay(x, h, mk, true_mass=True):
    ux, uz, gam = x[:N], x[N:2 * N], x[2 * N:]
    r = [r0[0], r0[1]]; v = [v0[0], v0[1]]; m = M0
    for k in range(N):
        mm = m if true_mass else mk[k]
        F = (mm * ux[k], mm * uz[k])
        Fn = math.hypot(*F)
        a = (F[0] / m, F[1] / m - G0)
        r = [r[0] + v[0] * h + 0.5 * a[0] * h * h, r[1] + v[1] * h + 0.5 * a[1] * h * h]
        v = [v[0] + a[0] * h, v[1] + a[1] * h]
        m -= alpha * Fn * h
    return r, v, m

print(f'planar powered descent for Hopper: start ({r0[0]:.0f}, {r0[1]:.0f}) m at ({v0[0]:.0f}, {v0[1]:.0f}) m/s, glide-slope half-angle {math.degrees(GAMMA):.0f} deg from the ground, thrust {TMIN:.0f} to {TMAX:.0f} N, {N} steps')
print('  descent time   status      fuel proxy (integral of Gamma, m/s)   barrier stages')
best = None
for T in (7.0, 8.0, 9.0, 10.0, 12.0, 16.0):
    x, info, mk, h = solve(T)
    if info['status'] == 'optimal':
        print(f'  {T:8.1f} s   optimal   {info["cost"]:14.2f}                          {info["iterations"]}')
        if best is None or info['cost'] < best[0]: best = (info['cost'], T, x, mk, h)
    else:
        print(f'  {T:8.1f} s   infeasible')
cost, T, x, mk, h = best
ux, uz, gam = x[:N], x[N:2 * N], x[2 * N:]
print(f'best of these: T = {T:.0f} s')
print(f'lossless: largest of Gamma_k - ||u_k|| = {max(gam[k] - math.hypot(ux[k], uz[k]) for k in range(N)):.2e} m/s2 (zero means the relaxed constraint is active at every step)')
Th = [mk[k] * gam[k] for k in range(N)]
ang = [math.degrees(math.atan2(ux[k], uz[k])) for k in range(N)]
print(f'thrust history (N): ' + ', '.join(f'{v:.0f}' for v in Th[::4]) + f' ... at the limits: {sum(1 for v in Th if v > TMAX - 5)} steps at {TMAX:.0f}, {sum(1 for v in Th if v < TMIN + 5)} steps at {TMIN:.0f}')
print(f'thrust angle from the vertical at the start {ang[0]:+.1f} deg, in the middle {ang[N // 2]:+.1f}, at the end {ang[-1]:+.1f}')
# the thrust limits were imposed with the nominal mass profile; check them against the true mass history
m = M0; worst_hi = 0.0; worst_lo = 0.0
for k in range(N):
    F = m * gam[k]
    worst_hi = max(worst_hi, F - TMAX); worst_lo = max(worst_lo, TMIN - F)
    m -= alpha * m * math.hypot(ux[k], uz[k]) * h
print(f'true propellant used: {M0 - m:.2f} kg (nominal profile assumed {mdot_nom * T:.2f} kg); with the true mass the thrust exceeds the {TMAX:.0f} N limit by at most {worst_hi:.1f} N and falls below the {TMIN:.0f} N floor by at most {worst_lo:.1f} N')
# glide-slope margin along the trajectory
r = [r0[0], r0[1]]; v = [v0[0], v0[1]]; marg = 1e9
for k in range(N):
    a = (ux[k], uz[k] - G0)
    r = [r[0] + v[0] * h + 0.5 * a[0] * h * h, r[1] + v[1] * h + 0.5 * a[1] * h * h]
    v = [v[0] + a[0] * h, v[1] + a[1] * h]
    marg = min(marg, r[1] - math.tan(GAMMA) * abs(r[0]))
print(f'smallest height above the glide-slope cone at the nodes: {max(marg, 0.0):.2f} m (the path touches the cone at the end)')

planar powered descent for Hopper: start (-150, 250) m at (10, -25) m/s, glide-slope half-angle 35 deg from the ground, thrust 720 to 2400 N, 24 steps
  descent time   status      fuel proxy (integral of Gamma, m/s)   barrier stages
       7.0 s   infeasible
       8.0 s   optimal           126.67                          10
       9.0 s   optimal           120.11                          10
      10.0 s   optimal           126.88                          10
      12.0 s   optimal           144.51                          10
      16.0 s   optimal           182.63                          10
best of these: T = 9 s
lossless: largest of Gamma_k - ||u_k|| = 5.37e-10 m/s2 (zero means the relaxed constraint is active at every step)
thrust history (N): 1785, 720, 720, 720, 2400, 2400 ... at the limits: 10 steps at 2400, 12 steps at 720
thrust angle from the vertical at the start +39.5 deg, in the middle +0.8, at the end -25.4
true propellant used: 9.43 kg (nominal profile assumed 8.56 kg); with the true mass the thrust exceeds the 2400 N limit by at most 39.8 N and falls below the 720 N floor by at most 0.0 N
smallest height above the glide-slope cone at the nodes: 0.00 m (the path touches the cone at the end)

In robotics

Lossless convexification is the basis of the onboard guidance of landing rockets, and of research on planetary landers.

• Flight use. Convex powered-descent guidance, derived from this formulation, was developed for planetary landers and has been flown on test vehicles; the key facts are the one-solution guarantee and the bounded solve time of the second-order cone program.
• Variable mass. The full problem tracks mass through its logarithm, whose derivative is minus the thrust over the exhaust speed, a linear relation, so that the thrust bounds become convex after a second-order expansion. The simplification used here, a fixed mass profile, leaves a 1.7% error in the thrust limit which a re-solve each second corrects.
• Time of flight. The convex problem is solved for a fixed time, and the best time is found by a one-dimensional search over the descent time, as the table of fuel against time shows: a smooth, single-minimum function.
• More constraints. A cone on the thrust direction (pointing), a limit on speed, an exclusion of the region below the glide slope are all convex and can be added without changing the class. A constraint requiring the vehicle to avoid an obstacle is not.
• Combining with feedback. The onboard cycle solves the problem from the current estimate, flies the first part, and solves again: the receding-horizon version, which is the subject of the next lesson.

The lab M9-L2: Solve the Landing as a Convex Program poses a divert-and-land as a quadratic program and re-solves it after a gust. The explicit law of M8-T3: Explicit Guidance to a Terminal State is the unconstrained version; compare its thrust demand with the constrained one, and the capstone M12-L2: Capstone: Land Hopper uses a guidance law of this family.

Common misconceptions

Key takeaways

1. Fuel-optimal landing has linear dynamics, a thrust-magnitude ceiling and a glide-slope cone, all convex, and a thrust floor, which is not.
2. The slack $\Gamma$ with $\|\mathbf{u}\| \le \Gamma$ and $\Gamma_{\min} \le \Gamma \le \Gamma_{\max}$ makes it a second-order cone program.
3. The minimum principle shows $\|\mathbf{u}\| = \Gamma$ at the optimum, because the Hamiltonian is linear in Γ and the costate is a straight line in time: the relaxation is lossless.
4. Hopper's fuel proxy has a minimum near 9 s (120.1 m/s) and the problem is infeasible at 7 s.
5. The optimal thrust is at the floor for 12 steps and at the ceiling for the last 11: bang-bang.
6. The fixed mass profile leaves 1.7% on the thrust ceiling; a log-mass formulation or re-solving removes it.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• B. Açıkmeşe and S. R. Ploen, "Convex programming approach to powered descent guidance for Mars landing", Journal of Guidance, Control, and Dynamics 30(5), 2007, the lossless convexification. Paid (AIAA).
• B. Açıkmeşe, J. M. Carson and L. Blackmore, "Lossless convexification of nonconvex control bound and pointing constraints of the soft landing optimal control problem", IEEE Transactions on Control Systems Technology 21(6), 2013. Paid (IEEE).
• L. Blackmore, "Autonomous precision landing of space rockets", The Bridge 46(4), 2016, National Academy of Engineering, an accessible account. Free to read.
• CVXPY documentation, tutorial , for posing second-order cone programs like this one with a modelling language. Free to read; licence Apache 2.0.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M9-T3: Convex Optimization for Guidance

  Recognise linear, quadratic and second-order cone programs, explain why convexity gives a bounded solve time and a certificate, and pose a guidance problem as one.

• Lab
  M9-L2: Solve the Landing as a Convex Program

  Pose a planar divert-and-land as a quadratic program with thrust and glide-slope limits, solve it with the course solver, and re-solve it in a receding horizon when a gust displaces the vehicle.

Other courses on this site treat neighbouring topics from their own angle:

• UAV M5-T8: Model Predictive Control
• UAV M5-T9: Trajectory generation
• UAV M5-T12: Kinodynamic planning
• Navigation M5-T6: LQR and MPC for Path Tracking
• Navigation M5-T3: Curvature, Velocity and Acceleration Limits
• Navigation M5-T2: Smoothing: Shortcuts, Splines and Gradient Methods
• Robot Learning M7-T5: Model Predictive Control with Learned Models
• Robot Learning M4-T1: Markov Decision Processes

This course leaves the following to the pages that already teach them:

• Minimum-snap polynomial trajectories and differential flatness of a quadrotor
• Path smoothing and velocity profiles for ground robots
• Sampling-based planning with a learned model (random shooting, CEM)
• Reinforcement learning