Fuel-Efficient 6-DoF Mars Powered Descent

Background

During terminal descent, a lander must remove horizontal position and velocity errors, keep its thrust direction and attitude within safe envelopes, and arrive upright without consuming excessive propellant. Translation, attitude, mass depletion, and thrust-vector moments are coupled, making this a useful nonlinear guidance benchmark.

This example plans a 16-second Mars descent from 100 m altitude. It uses a three-component modified Rodrigues parameter (MRP) attitude state instead of a redundant quaternion, and scales large translational and propulsion quantities before transcription.

Rigid-body and propulsion model

Coordinates are ordered as up, east, and north. Let \(\boldsymbol r\) and \(\boldsymbol v\) be inertial position and velocity, \(m\) the remaining mass, \(\boldsymbol p\) the MRP vector, \(\boldsymbol\omega\) the body angular velocity, and \(\boldsymbol T_b\) the body-frame thrust command. The body \(x\) axis points through the engine.

With body-to-inertial rotation \(R(\boldsymbol p)\), exhaust velocity \(v_e\), body inertia \(I_b\), and engine offset \(\boldsymbol\rho_e\), the physical dynamics are

\[ \begin{aligned} \dot m&=-\frac{\lVert\boldsymbol T_b\rVert_2}{v_e},\\ \dot{\boldsymbol r}&=\boldsymbol v,\\ \dot{\boldsymbol v} &=\frac{R(\boldsymbol p)\boldsymbol T_b}{m} -g_M\boldsymbol e_u,\\ \dot{\boldsymbol p} &=\frac14\left[ (1-\boldsymbol p^{\mathsf T}\boldsymbol p)I +2[\boldsymbol p]_\times +2\boldsymbol p\boldsymbol p^{\mathsf T} \right]\boldsymbol\omega,\\ \dot{\boldsymbol\omega} &=I_b^{-1}\left[ \boldsymbol\rho_e\times\boldsymbol T_b -\boldsymbol\omega\times(I_b\boldsymbol\omega) \right]. \end{aligned} \]

The MRP-to-rotation conversion is generated through

\[ q_0=\frac{1-\lVert\boldsymbol p\rVert_2^2} {1+\lVert\boldsymbol p\rVert_2^2}, \qquad \boldsymbol q=\frac{2\boldsymbol p} {1+\lVert\boldsymbol p\rVert_2^2}, \]

followed by the corresponding rotation matrix. The quaternion components are only an algebraic construction; they are not optimization states.

Objective and endpoint conditions

The fixed horizon is \(T=16\ \mathrm s\). Since mass flow is proportional to thrust magnitude, the scaled minimum-thrust-impulse objective is equivalent to minimum propellant:

\[ \min J=\int_0^T \frac{\lVert\boldsymbol T_b(t)\rVert_2}{F_s}\,\mathrm dt, \qquad m(0)-m(T)=\frac{F_sJ}{v_e}, \qquad F_s=2000\ \mathrm N. \]

The initial conditions are

\[ \begin{aligned} m(0)&=2000\ \mathrm{kg},\\ \boldsymbol r(0)&=(100,30,-15)^{\mathsf T}\ \mathrm m,\\ \boldsymbol v(0)&=(-12,-4,2)^{\mathsf T}\ \mathrm{m/s},\\ \boldsymbol p(0)&=(0,0.04,-0.03)^{\mathsf T},\\ \boldsymbol\omega(0)&=\boldsymbol 0. \end{aligned} \]

At touchdown, mass is free while position, velocity, attitude, and rate are fixed:

\[ \boldsymbol r(T)=\boldsymbol 0,\qquad \boldsymbol v(T)=(-0.5,0,0)^{\mathsf T}\ \mathrm{m/s},\qquad \boldsymbol p(T)=\boldsymbol 0,\qquad \boldsymbol\omega(T)=\boldsymbol 0. \]

Thus the contact speed is a prescribed soft vertical descent, not a zero-speed hover endpoint.

Path constraints

The physical limits are

\[ \begin{aligned} 1400&\le m(t)\le2000\ \mathrm{kg},\\ 4000&\le\lVert\boldsymbol T_b(t)\rVert_2\le30000\ \mathrm N,\\ T_{b,x}&\ge \cos(15^\circ)\lVert\boldsymbol T_b\rVert_2,\\ \sqrt{r_e^2+r_n^2}&\le\tan(30^\circ)r_u,\qquad r_u\ge0,\\ R_{ux}(\boldsymbol p)&\ge\cos(25^\circ),\\ \lVert\boldsymbol\omega\rVert_2&\le0.70\ \mathrm{rad/s},\\ \lVert\boldsymbol p\rVert_2&\le\tan(120^\circ/4). \end{aligned} \]

The third inequality is the \(15^\circ\) engine-gimbal cone, the fourth is a \(30^\circ\) glide-slope cone about the landing site, and \(R_{ux}\) is the alignment between inertial up and the body \(x\) axis. The transcription applies a 10 N guard to the minimum thrust at nodes; this compensates for interpolation between nodes, while dense validation still uses the physical 4000 N limit.

Variables, parameters, and units

Symbol	Meaning	Value or unit
\(m\)	Lander mass	kg
\(\boldsymbol r,\boldsymbol v\)	Up-east-north position and velocity	m, m/s
\(\boldsymbol p\)	Modified Rodrigues attitude vector	dimensionless
\(\boldsymbol\omega\)	Body angular velocity	rad/s
\(\boldsymbol T_b\)	Body-frame thrust vector	N
\(g_M\)	Mars gravitational acceleration	\(3.711\ \mathrm{m/s^2}\)
\(v_e\)	Effective exhaust velocity	\(2250\ \mathrm{m/s}\)
\(I_b\)	Body inertia diagonal	\((4000,2500,2500)\ \mathrm{kg\,m^2}\)
\(\boldsymbol\rho_e\)	Engine offset from center of mass	\((-1.5,0,0)\ \mathrm m\)

Internally, length, velocity, mass, and force are divided by \(100\ \mathrm m\), \(10\ \mathrm{m/s}\), \(2000\ \mathrm{kg}\), and \(2000\ \mathrm N\), respectively. All input data, constraints, diagnostics, and plot labels remain in SI units.

Modeling choices and limits

The MRP state is minimal and avoids a quaternion unit-norm constraint. It is a local chart, not a global attitude coordinate: the example does not switch to the MRP shadow set. Bounding the principal attitude angle by \(120^\circ\) keeps the entire trajectory far from the \(360^\circ\) singularity. The model is therefore intended for a controlled terminal descent, not a tumbling entry.

The model assumes constant gravity, constant inertia and center of mass, a single ideal throttling engine, instantaneous thrust-vector commands, and no aerodynamic force, plume interaction, terrain variation, navigation error, or uncertainty. The result is a deterministic open-loop reference trajectory. The default transcription uses a \(40\times2\) Lobatto mesh.

Run the example

From the repository root:

python -m examples.rocket_powered_descent

Save the figure without opening a window:

python -m examples.rocket_powered_descent --save rocket-powered-descent.png --no-show

Verified result

The script reconstructs physical states and thrust at 4,001 times. It checks touchdown position, velocity, attitude, and angular rate; dry mass; thrust; gimbal, tilt, glide-slope, altitude, angular-rate, and local-chart bounds; and agreement between integrated thrust impulse and mass depletion.

Quantity	Default-solve value
Flight time	\(16.000\ \mathrm s\)
Propellant used	\(62.219\ \mathrm{kg}\)
Vertical contact velocity	\(-0.500000\ \mathrm{m/s}\)
Dense thrust range	\(4005.25\) to \(30000.00\ \mathrm N\)
Peak gimbal angle	\(9.0869^\circ\)
Peak vehicle tilt	\(11.4496^\circ\)
Minimum glide-cone margin	\(0.000000\ \mathrm m\)

The zero reported glide-cone margin means that part of the locally optimal trajectory rides the conic boundary; it does not mean the constraint was omitted. At default accuracy, terminal position and velocity errors must each remain below \(5\times10^{-5}\) in their SI units, and the dense path check must remain below \(5\times10^{-4}\).

Source code

See the complete runnable example: examples/rocket_powered_descent.py.