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}$$

气动参数

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

info

对于飞行器而言, 气动系数一般是随着速度、高度等环境因素变化的, OPTIMake将有专门处理 LookupTable函数支持这类输入, 并实现与求解器的高效集成.

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.