地球至金星最小推进剂转移

背景

低推力电推进能够以较少的推进剂提供很大的累计速度增量，但较小的 加速度必须持续作用数月。因此，推力方向和节流量需要在整条轨迹上 统一规划。本例要求一艘初始质量为 \(1{,}500\ \mathrm{kg}\) 的航天器 在固定的 1,000 天内从地球转移至金星，并使推进剂消耗最小。

日心轨迹采用改进春分点轨道根数（modified equinoctial elements， MEE）表示。与经典轨道根数不同，MEE 对圆形、赤道顺行轨道保持正则， 适合构造直接配点模型。

问题定义

日心笛卡尔边界数据采用天文单位和儒略年：

边界向量	数值	单位
\(\boldsymbol r_0\)	\([0.9708322,\ 0.2375844,\ -1.671055\times10^{-6}]\)	AU
\(\boldsymbol v_0\)	\([-1.598191,\ 6.081958,\ 9.443368\times10^{-5}]\)	AU/year
\(\boldsymbol r_f\)	\([-0.3277178,\ 0.6389172,\ 2.765929\times10^{-2}]\)	AU
\(\boldsymbol v_f\)	\([-6.598211,\ -3.412933,\ 0.3340902]\)	AU/year

\(\boldsymbol r_0\) 的第一个分量为正号，这是有意的，并且与 2005-10-07 当天地球的日心经度一致。这个基准有时会写成负号，但该 符号与给定速度和日期不一致。

状态向量和优化控制向量为

\[ \boldsymbol x=[p,f,g,h,k,L,m]^{\mathsf T}, \qquad \boldsymbol q=[d_r,d_t,d_n,\rho]^{\mathsf T}. \]

其中，\(p\) 是半通径，\(f,g,h,k\) 描述偏心率和轨道倾角， \(L\) 是真经度，\(m\) 是质量。向量 \(\boldsymbol d=[d_r,d_t,d_n]^{\mathsf T}\) 位于带保护量的单位球内， 非负变量 \(\rho\) 是节流量。实际归一化 RTN 指令为 \(\boldsymbol u=\rho\boldsymbol d=[u_r,u_t,u_n]^{\mathsf T}\)。

定义

\[ \begin{aligned} w &= 1+f\cos L+g\sin L, & s^2 &= 1+h^2+k^2,\\ \chi &= h\sin L-k\cos L, & \gamma &= \sqrt{\frac{p}{\mu}},\\ [a_r,a_t,a_n]^{\mathsf T} &=\frac{T_{\max}}{m}[u_r,u_t,u_n]^{\mathsf T}. \end{aligned} \]

MEE 和质量动力学为

\[ \begin{aligned} \dot p &= \frac{2p\gamma}{w}a_t,\\ \dot f &= \gamma\left[ a_r\sin L+\frac{((w+1)\cos L+f)a_t-\chi g a_n}{w}\right],\\ \dot g &= \gamma\left[ -a_r\cos L+\frac{((w+1)\sin L+g)a_t+\chi f a_n}{w}\right],\\ \dot h &= \frac{\gamma s^2\cos L}{2w}a_n,\\ \dot k &= \frac{\gamma s^2\sin L}{2w}a_n,\\ \dot L &= \sqrt{\mu p}\left(\frac{w}{p}\right)^2 +\frac{\gamma\chi}{w}a_n,\\ \dot m &= -c\rho, \qquad c=\frac{T_{\max}}{I_{\mathrm{sp}}g_0}. \end{aligned} \]

发动机参数为 \(T_{\max}=0.33\ \mathrm N\) 和 \(I_{\mathrm{sp}}=3{,}800\ \mathrm s\)。方向与节流约束为

\[ d_r^2+d_t^2+d_n^2\le0.9999^2,\qquad 0\le\rho\le1. \]

因此

\[ u_r^2+u_t^2+u_n^2 =\rho^2(d_r^2+d_t^2+d_n^2) \le0.9999^2\rho^2\le\rho^2. \]

很小的相对保护量可使数值约束误差仍落在物理锥内，同时不需要正的 节流下限：\(\rho=0\) 仍精确对应零推力和零质量流率。

初始质量固定，并在质量方程中使用同一个 \(\rho\)，因此最小化

\[ J=\int_0^{t_f}c\rho\,\mathrm dt=m(0)-m(t_f) \]

就等价于最小化推进剂消耗。终端 MEE 固定为金星状态，终端质量自由。 真经度额外展开三个完整周次，使配点轨迹沿预期的连续分支变化。

变量与单位

符号	含义	单位
\(t\)	从出发时刻起算的时间	儒略年
\(p\)	半通径	AU
\(f,g,h,k\)	描述轨道形状和方向的春分点根数	-
\(L\)	展开的真经度	rad
\(m\)	航天器质量	kg
\(d_r,d_t,d_n\)	带保护量的 RTN 方向向量	-
\(u_r,u_t,u_n\)	归一化 RTN 推力指令	-
\(\rho\)	发动机节流量	-

脚本先将 \(T_{\max}\)、\(I_{\mathrm{sp}}\) 和 \(g_0\) 的 SI 数值转换 为 AU/year/kg 单位，再构造符号动力学。

建模选择

笛卡尔状态到 MEE 的转换直接使用角动量向量和偏心率向量构造。这样 避免了连续使用反余弦、手动修正象限，以及在接近圆轨道时除以很小的 偏心率。求解前，脚本会把两个边界状态从 MEE 转回笛卡尔坐标进行 往返检查。

把 \(\rho\) 与 \(\boldsymbol d\) 分开也是有意的。若直接把质量流率 写成 \(\sqrt{u_r^2+u_t^2+u_n^2}\)，函数在零推力处不可微。乘积 \(\boldsymbol u=\rho\boldsymbol d\) 和二次方向球都是光滑的，保留了 推力锥内部，并改善锥顶的数值性质：当 \(\rho=0\) 时，无论方向变量 为何值，物理指令都精确为零。

求解先在 \(12\times3\) 的 Radau 网格上生成可靠热启动。两种模式随后 都适配到 120 个单点 Radau 区间；快速模式在此结束，默认模式再求解 一次 320 个单点区间。单点 Radau 控制在每个区间内为常数，状态为 线性函数。因此节流量和方向边界在节点之间仍然成立，而 \(\dot m=-c\rho\) 使插值质量单调不增。最终低阶表示避免了高阶 Radau 控制多项式在未受约束的区间右端发生外推。

每次最终求解后，脚本都会在 10,001 个等间隔时刻计算插值结果。只要 节流量离开 \([0,1]\)、实际 RTN 指令离开锥、质量增加，或者已有的 \(p\)、\(m\)、\(w\) 路径边界被违反，就会显式抛出异常。绘图直接使用 这些通过验证的插值结果，不会裁剪不可行控制。

运行示例

在仓库根目录执行经过验证的默认求解：

python -m examples.orbit_transfer

使用较小的最终网格：

python -m examples.orbit_transfer --quick

保存默认结果且不打开绘图窗口：

python -m examples.orbit_transfer --save orbit-transfer.png --no-show

核心实现

system = System(0, fastmath=True)
phase = system.new_phase(
    ["p", "f", "g", "h", "k", "longitude", "mass"],
    ["direction_radial", "direction_transverse", "direction_normal", "rho"],
)
d_r, d_t, d_n, rho = phase.u
u_r, u_t, u_n = rho * d_r, rho * d_t, rho * d_n

# 完整脚本定义了上文给出的六个 MEE 导数。
phase.set_dynamics([dot_p, dot_f, dot_g, dot_h, dot_k, dot_L, -c * rho])
phase.set_integral([c * rho])

direction_squared = d_r**2 + d_t**2 + d_n**2
phase.set_phase_constraint(
    [p, mass, w, rho, d_r, d_t, d_n, direction_squared],
    [0.20, 500.0, 0.10, 0.0, -0.9999, -0.9999, -0.9999, 0.0],
    [2.00, 1500.0, np.inf, 1.0, 0.9999, 0.9999, 0.9999, 0.9999**2],
)
phase.set_boundary_condition(
    [*INITIAL_MEE, 1500.0],
    [*TARGET_MEE, None],
    0.0,
    1000.0 / 365.25,
)
phase.set_discretization(12, 3)
system.set_phase([phase])
system.set_objective(phase.I[0])

# 热启动求解完成后：
phase.set_discretization(120, 1)  # 快速模式的最终网格
# 默认模式还会适配到 phase.set_discretization(320, 1)。

完整脚本还包括坐标转换、物理常量、向内螺旋初始猜测、分阶段求解、 10,001 点验证和绘图。

验证结果

经过验证的结果为：

模式	最终网格	推进剂消耗	终端质量
快速	\(120\times1\)	\(209.872834\ \mathrm{kg}\)	\(1290.127166\ \mathrm{kg}\)
默认	\(320\times1\)	\(209.720375\ \mathrm{kg}\)	\(1290.279625\ \mathrm{kg}\)

默认解的全部 10,001 个样本均满足 \(0\le\rho\le1\)。最小物理锥裕量 \(\rho-\lVert\boldsymbol u\rVert_2\) 为 \(3.162\times10^{-12}\)，最小二次锥裕量为 \(5.608\times10^{-23}\)，输出精度下的最大方向模为 \(0.9999\)。 每 0.1 天的最大质量变化为 \(-6.821\times10^{-13}\ \mathrm{kg}\)，没有任何采样区间出现质量增加。

密集检查得到半通径范围 \([0.723300547,1.000381785]\ \mathrm{AU}\)，质量范围 \([1290.279625,1500]\ \mathrm{kg}\)，最小半径分母为 \(w=0.949954686\)。输出精度下的最大终端 MEE 误差为零。在本次验证 机器上，快速和默认命令分别耗时约 137 秒和 157 秒；实际耗时会随编译 缓存和硬件变化。

源代码

完整可运行示例见： examples/orbit_transfer.py。