Reusable Rocket Landing Problem

Problem Description

The objective function of starship-like landing problem is defined as the fuel-optimal, and its performance index can be expressed as

$$\begin{split} \rm{J} &= -m(t_f) \end{split}$$

Starship-like vehicle dynamics

$$\begin{split} \dot{x}&=V_x \\ \dot{z}&=V_z \\ \dot{V_x}&=\tfrac{-T \sin(\theta + \delta) + D_x}{m} \\ \dot{V_z}&=\tfrac{T \cos(\theta + \delta) + D_z}{m} - g_0 \\ \dot{\theta}&=\omega \\ \dot{\omega}&=\tfrac{M_T + M_D}{J} \\ \dot{m}&=-\tfrac{T}{I_{\text{sp}} g_0} \end{split}$$

$D_x$ and $D_z$ are the total aerodynamic drag in the $x$ and $z$ directions, respectively, described as

$$\begin{split} D_x &= -C_{LD} \sqrt{V_x^2 + V_z^2} \cdot V_x \\ D_z &= -C_{LD} \sqrt{V_x^2 + V_z^2} \cdot V_z \\ \end{split}$$

The terminal state constraints that satisfy the fixed-point vertical soft landing as follows.

$$\begin{split} J &= \tfrac{1}{12} m (6 r_s^2 + l_s^2) \\ M_T &= -l_{cg} \cdot T \cdot \sin \delta\\ M_D &= -(l_{cp} - l_{cg}) \cdot (D_x \cdot \cos \theta + D_z \cdot \sin \theta) \end{split}$$

The initial state constraints and the terminal state constraints that satisfy the fixed-point vertical soft landing as follows:

$$\begin{split} x(0) &= x_0, \ x(t_f) = x_{t_f}, \ m(t_f) \geq m_{\text{dry}} \\ \end{split}$$

For the control constraints, there are thrust amplitude constraints and engine nozzle swing angle amplitude constraints:

$$\begin{split} T_{\text{min}} \leq T \leq T_{\text{max}} \\ \delta{\text{min}} \leq \delta \leq \delta_{\text{max}} \\ \end{split}$$

Aerodynamic Parameters

Parameter	$m_{\text{dry}}$, $m_0$	$L, r$	$g_0$	$C_{LD0}$	$I_{\text{sp}}$
Value	8.5e4, 1e5 (kg)	50, 3 $(m)$	9.82	2.16	330
Parameter	$L{\text{cg}}$	$L{\text{cp}}$	$T_{\text{max}}$	$\text{throttle}_{\text{max}}$	$\text{throttle}_{\text{min}}$
Value	20	20.5	2.21e6(N)	1.0	0.4

info

For aerospace vehicles, aerodynamic coefficients generally vary with environmental factors such as speed and altitude. OPTIMake will provide dedicated LookupTable function support for this type of input, enabling efficient integration with the solver.

Modeling

NN = 100
prob = multi_stage_problem('starship_full', NN + 1)

m0 = 1e5 # initial mass of rocket (kg)
mdry = 8.5e4 # dry weight of rocket (kg)

x0, y0, theta0, vx0, vy0, dtheta0 = prob.parameters(['x0', 'y0', 'theta0', 'vx0', 'vy0', 'dtheta0'],  stage_dependent=False)

x = prob.variable('x')
y = prob.variable('y', hard_lowerbound=0.0)
theta = prob.variable('theta')
vx = prob.variable('vx')
vy = prob.variable('vy')
dtheta = prob.variable('dtheta')
mass = prob.variable('mass', hard_lowerbound=mdry/m0) # normalized mass

tf = prob.variable('tf', hard_lowerbound=0.0, hard_upperbound=50.0)

angle_bound = 20.0 / 180.0 * pi
trust_angle = prob.variable('trust_angle', hard_lowerbound=-angle_bound, hard_upperbound=angle_bound)
trust = prob.variable('trust', hard_lowerbound=0.4, hard_upperbound=1.0)

""" ode """
hs = 7110
rou0 = 1.225
beta = -1 / hs
rou = rou0 * exp(beta * y)
CLD = 2.16 * rou # drag coefficient
Dx = -CLD * vx * sqrt(vx**2 + vy**2) # add 1e-3 to avoid singularity when vx = vy = 0 (backward euler, trapzoid)
Dy = -CLD * vy * sqrt(vx**2 + vy**2)

L = 50 # length of rocket (m)
r = 3 # radius of rocket (m)
g = 9.806 # gravity (m/s^2)
Isp = 330 # specific impulse (s)
I = 1/12*(mass * m0)*(3 * r**2 + L**2)  # inertia
max_trust = 2.21e6 # max trust (N)
T = max_trust * trust

Lcg = 20 # center of mass position (m)
Lcp = 22.5 # center of pressure position (m)
MT = -Lcg * T * sin(trust_angle)
MD = -(Lcp - Lcg) * (-Dx * cos(theta) + Dy * sin(theta))
Mz = MT + MD

Fx = T * sin(trust_angle + theta) + Dx
Fy = T * cos(trust_angle + theta) + Dy
dxdt = [vx, vy, dtheta, Fx/(mass * m0), Fy/(mass * m0) - g, Mz/I, -T/(Isp*g*m0), 0]
ts = tf / NN
ode = differential_equation(
    state=[x, y, theta, vx, vy, dtheta, mass, tf],
    state_dot=dxdt, 
    stepsize=ts, 
    discretization_method='forward_euler')
prob.equality(ode)

obj = general_objective(-mass)
# minimize fuel consumption
prob.end_objective(obj)

seq = general_equality([x - x0, y - y0, theta - theta0, vx - vx0, vy - vy0, dtheta - dtheta0, mass - 1])
prob.start_equality(seq)

eeq = general_equality([x - 0, y - 0, theta - 0, vx - 0, vy - 0, dtheta - 0])
prob.end_equality(eeq)

option = codegen_option()
codegen = code_generator()
codegen.codegen(prob, option)

Solution

---

[1]. Chen, H.; Ma, Z.; Wang, J.; Su, L. Online Trajectory Optimization Method for Large Attitude Flip Vertical Landing of the Starship-like Vehicle. Mathematics 2023, 11, 288.