跳到主要内容

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

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

Starship-like vehicle dynamics

x˙=Vxz˙=VzVx˙=−Tsin⁡(θ+δ)+DxmVz˙=Tcos⁡(θ+δ)+Dzm−g0θ˙=ωω˙=MT+MDJm˙=−TIspg0\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}

DxD_x and DzD_z are the total aerodynamic drag in the xx and zz directions, respectively, described as

Dx=−CLDVx2+Vz2⋅VxDz=−CLDVx2+Vz2⋅Vz\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.

J=112m(6rs2+ls2)MT=−lcg⋅T⋅sin⁡δMD=−(lcp−lcg)⋅(Dx⋅cos⁡θ+Dz⋅sin⁡θ)\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:

x(0)=x0, x(tf)=xtf, m(tf)≥mdry\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:

Tmin≤T≤Tmaxδmin≤δ≤δmax\begin{split} T_{\text{min}} \leq T \leq T_{\text{max}} \\ \delta{\text{min}} \leq \delta \leq \delta_{\text{max}} \\ \end{split}

气动参数

参数mdrym_{\text{dry}}, m0m_0L,rL, rg0g_0CLD0C_{LD0}IspI_{\text{sp}}
值8.5e4, 1e5 (kg)50, 3 (m)(m)9.822.16330
参数LcgL{\text{cg}}LcpL{\text{cp}}TmaxT_{\text{max}}throttlemax\text{throttle}_{\text{max}}throttlemin\text{throttle}_{\text{min}}
值2020.52.21e6(N)1.00.4
信息

对于飞行器而言, 气动系数一般是随着速度、高度等环境因素变化的, 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​

starship_full_traj starship_full

[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.