航天器轨迹优化可视化异常：仅显示起止点无轨迹求助

问题背景 #

正在开展航天器运动优化研究，航天器初始位于相对O点的近地无扰圆轨道（轨道半径6871km），地球引力常数为3.9860044×10^14 m³/s²。任务要求航天器在86400秒内，从初始位置r0=(10 km, 100 km, -5 km)^T、初始速度v0=(1 m/s, -10 m/s, 3 m/s)^T飞抵原点（rf=0，vf=0）。航天器初始质量1000kg，推进系统比冲220秒（2157.463 m/s），推力范围0.362~100N，目标是最小化推进剂质量完成轨道转移与交会对接。

使用GEKKO编写优化代码后，求解得到结果，但绘图时仅能看到起始点和终点，无法显示飞行轨迹。

问题原因 #

1. 尺度不匹配：地球模型的尺度（半径6871km）远大于航天器的相对轨迹（最大初始位置仅100km），轨迹线被地球模型完全覆盖，肉眼无法分辨。
2. 单位不一致：代码中位置数据以米为单位，但坐标轴标签标注为千米，显示逻辑混乱。
3. 时间维度未映射：代码使用归一化时间[0,1]，未对应实际任务时间86400秒，虽不影响绘图，但降低了物理意义的明确性。

解决方案 #

• 调整绘图坐标系范围，聚焦航天器的相对轨迹；
• 统一单位，将位置数据转换为千米后绘图；
• 可选：将归一化时间还原为实际任务时间，强化物理意义；
• 缩小地球模型比例（或直接隐藏），避免覆盖轨迹。

修改后的完整代码 #

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create GEKKO model
m = GEKKO()

# Time discretization - 还原为实际任务时间86400秒
nt = 101
t_final = 86400
m.time = np.linspace(0, t_final, nt)

# Constants
mu = 3.9860044e14
R = 6871e3  # 地球半径，单位m
Isp = 2157.463
g0 = 9.81

# Initial conditions - 单位m
r0 = np.array([10e3, 100e3, -5e3])
v0 = np.array([1, -10, 3])
rf = np.array([0, 0, 0])
vf = np.array([0, 0, 0])
m0 = 1000

# Variables
r = [m.Var(value=r0[i], name=f'r{i}') for i in range(3)]
v = [m.Var(value=v0[i], name=f'v{i}') for i in range(3)]
mass = m.Var(value=m0, name='mass')
thrust = m.MV(value=50, lb=0.362, ub=100, name='thrust')
thrust.STATUS = 1
theta = [m.Var(value=0, lb=-np.pi, ub=np.pi, name=f'theta{i}') for i in range(3)]

# Equations of motion
for i in range(3):
    m.Equation(r[i].dt() == v[i])
    m.Equation(v[i].dt() == (thrust * m.cos(theta[i])) / mass - (mu / R**3) * r[i])

# Mass flow rate
m_flow = thrust / (Isp * g0)
m.Equation(mass.dt() == -m_flow)

# Final conditions
for i in range(3):
    m.fix(r[i], pos=nt-1, val=rf[i])
    m.fix(v[i], pos=nt-1, val=vf[i])

# Objective: Minimize the final mass (等价于最大化推进剂消耗)
m.Obj(mass)

# Solve the optimization problem
m.options.IMODE = 6
m.solve(disp=True)

# Display results
print(f'Final mass: {mass.value[-1]:.2f} kg')
for i in range(3):
    print(f'Final position r{i}: {r[i].value[-1]:.2f} m')
    print(f'Final velocity v{i}: {v[i].value[-1]:.2f} m/s')

# Extract trajectory data - 转换为km便于绘图
r_data_km = np.zeros((nt, 3))
for i in range(3):
    r_data_km[:, i] = np.array(r[i].value) / 1000

# Plot the trajectory in 3D
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# Plot the trajectory
ax.plot(r_data_km[:, 0], r_data_km[:, 1], r_data_km[:, 2], linewidth=2, label='Trajectory')

# 可选：绘制缩小比例的地球模型，避免覆盖轨迹
earth_scale = 0.1  # 缩小为原尺寸的10%
u, v_sphere = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
x_earth = (R/1000)*earth_scale * np.cos(u) * np.sin(v_sphere)
y_earth = (R/1000)*earth_scale * np.sin(u) * np.sin(v_sphere)
z_earth = (R/1000)*earth_scale * np.cos(v_sphere)
ax.plot_surface(x_earth, y_earth, z_earth, color='w', alpha=0.3)

# Plot the initial and final positions (转换为km)
ax.scatter(r0[0]/1000, r0[1]/1000, r0[2]/1000, color='r', s=50, label='Initial Position')
ax.scatter(rf[0]/1000, rf[1]/1000, rf[2]/1000, color='g', s=50, label='Final Position')

# Labels and legend
ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Z (km)')
ax.set_title('Spacecraft Relative Trajectory')
ax.legend()

# 自动调整坐标轴范围，聚焦轨迹
ax.set_xlim(np.min(r_data_km[:,0])-5, np.max(r_data_km[:,0])+5)
ax.set_ylim(np.min(r_data_km[:,1])-5, np.max(r_data_km[:,1])+5)
ax.set_zlim(np.min(r_data_km[:,2])-5, np.max(r_data_km[:,2])+5)

plt.show()

关键修改说明 #

• 时间维度修正：将m.time从归一化区间改为实际任务时间区间，确保动力学方程的时间尺度正确；
• 单位统一：将位置数据从米转换为千米，与坐标轴标签保持一致；
• 坐标系聚焦：通过set_xlim/set_ylim/set_zlim手动设置坐标轴范围，确保轨迹完全显示；
• 地球模型缩放：缩小地球模型比例，避免其覆盖相对轨迹（若不需要可直接删除地球模型绘制代码）。

内容的提问来源于stack exchange，提问作者Angola XXI