不变流形与庞加莱截面

不变流形(稳定/不稳定)是周期轨道附近的渐近轨道族,是低能量转移设计的 基本构件。本页介绍流形计算、庞加莱截面工具,以及基于流形拼接的低能 转移流水线。

不变流形计算

种子生成原理:周期轨道单值矩阵 M(沿轨道传播一周的 STM)的实特征值给出 流形方向,稳定流形取 abs(λ) < 1 的实特征向量,不稳定流形取 abs(λ) > 1 的 实特征向量; 单位圆上的特征值(含 λ=1 的周期方向)不构成双曲方向,予以剔除。 沿轨道均匀取 n_points 个相位点,用从首点到各相位的 STM 把特征向量转运到 该相位,位置部分归一化后施加 ±ε 的无量纲扰动得到种子;稳定流形反向积分、 不稳定流形正向积分得到流形管。

InvariantManifold 的用法:

from e2m2e.algorithm.manifold import InvariantManifold, ManifoldKind

# orbit 为周期轨道:须关联 system 且 period 已知(可只存首点)
epsilon = 50.0 / 384405.0   # 无量纲扰动幅度,典型取 50 km / DU
manifold = InvariantManifold(orbit, ManifoldKind.STABLE, "-", epsilon)

# 相位扫掠种子,形状 (n_points, 6)
seeds = manifold.seeds(12)

# 批量传播流形弧(稳定流形反向积分,t_span 取绝对值)
tube = manifold.propagate(4.0)
print(f"流形弧数: {len(tube.trajectories)}")

branch"+""-",对应扰动的两个方向,分别走向轨道两侧。 返回的 ManifoldTube 携带轨道引用、 流形类型、分支与 ε,trajectories 为流形弧列表(无量纲 CR3BP 态)。

propagate 传入 section 参数时,每条弧在首次穿越截面处截断, 并把求精后的穿越态追加为弧的末点:

from e2m2e.algorithm.manifold import PoincareSection

section = PoincareSection.periapsis("earth", orbit.system)
tube = manifold.propagate(4.0, section=section)

庞加莱截面

PoincareSection 由标量函数 s(state) 的 零等值面定义,提供两类构造:

  • PoincareSection.plane(axis, value):平面截面 s = state[axis] − value, axis 为状态分量索引(0=x, 1=y, 2=z, 3=vx, 4=vy, 5=vz)

  • PoincareSection.periapsis(center, system):近拱点截面 s = r·v, r 为相对 center 天体(主/次天体名称,不区分大小写)的位置

穿越检测采用事后方案:传播时密采样(流形传播默认步长 0.005 无量纲时间), 逐采样点求截面函数值,符号变化区间内对分段线性插值态用 Brent 法求精。 平面截面穿越残差可达 1e-10 以下。

crossings() 检测流形管中所有弧的穿越,返回 SectionCrossings

crossings = section.crossings(tube)
print(crossings.states.shape)          # (k, 6),插值求精后的穿越态
print(crossings.times.shape)           # (k,),穿越时刻
print(crossings.trajectory_index)      # (k,),每个穿越点所属弧的索引

crossings() 是事后检测(先传播、再在采样点上找穿越)。若要在积分 过程中检测穿越(例如首次到达截面即停),用 event() 生成 scipy 语义的 事件函数传给 Dynamics.propagate(events=...)

event = section.event(direction=-1, terminal=True)
result = dynamics.propagate(y0, (0.0, 10.0), events=[event])

详见 动力学 的「事件检测」一节。

流形拼接与低能转移

两条流形管在同一截面上的穿越点两两配对,可拼接成转移初猜。 patch_manifolds() 按加权拼接代价 w_r·|Δr| + w_v·|Δv| 升序输出候选(PatchCandidate):

from e2m2e.algorithm.transfer import patch_manifolds

# 出发轨道不稳定流形 + 目标轨道稳定流形,同一截面
candidates = patch_manifolds(tube_a, tube_b, section, weights=(1.0, 1.0))
best = candidates[0]
print(f"|Δr|={best.delta_r:.4e}, |Δv|={best.delta_v:.4e}")

design_low_energy_transfer() 把上述步骤 串成流水线:出发轨道不稳定流形与目标轨道稳定流形(± 分支四种组合全局 取最优)传播到次天体近拱点截面,取最优拼接候选,出发弧直接用流形弧, 拼接点之后由 ThreeBodyLambert 打靶闭合到目标轨道。脉冲由三段构成:出发脉冲(上出发流形)、拼接脉冲 (截面处)、到达脉冲(入目标流形)。

from e2m2e.algorithm.transfer import OrbitTerminal, design_low_energy_transfer
from e2m2e.data.templates import ConvergenceState

sol = design_low_energy_transfer(OrbitTerminal(departure_orbit), target_orbit)

if sol.status == ConvergenceState.CONVERGED:
    print(f"弧段数: {len(sol.arcs)}")           # 2
    print(f"出发脉冲: {sol.arcs[0].delta_v:.6f} km/s")
    print(f"拼接脉冲: {sol.arcs[1].delta_v:.6f} km/s")
    print(f"到达脉冲: {sol.arrival_delta_v:.6f} km/s")
    print(f"总脉冲: {sol.total_delta_v:.6f} km/s")
    print(f"转移时间: {sol.transfer_time:.1f} s")

返回两段弧的 TransferSolution,物理单位。 当前仅支持 CR3BP 模型;星历转换(CR3BP 闭合解 → 星历模型)尚未接入, epoch 参数为其预留入口。端到端基准见 tests/algorithm/transfer/test_low_energy.py: L1 Lyapunov 族内中间轨道到大幅值轨道,拼接脉冲在几十 m/s 量级。

不变流形计算模块

提供 CR3BP 周期轨道不变流形(稳定/不稳定)的种子生成与批量传播功能。

算法要点

单值矩阵 M(沿轨道传播一周的 STM)的实特征值给出流形方向: 稳定流形取模小于 1 的实特征向量,不稳定流形取模大于 1 的实特征向量。 沿周期轨道取 n_points 个相位点,用从轨道首点到各相位的 STM 把特征向量 转运到该相位,位置部分归一化后施加 ±ε 的无量纲扰动得到种子; 稳定流形反向积分、不稳定流形正向积分得到流形管。

数值内核在 Rust(特征分解、STM 转运、种子扰动、批量传播调度); Python 侧只做参数校验、领域对象组装与可选的事后截面截断。

class e2m2e.algorithm.manifold.manifolds.ManifoldKind(*values)[源代码]

基类:Enum

不变流形类型枚举

STABLE = 'stable'
UNSTABLE = 'unstable'
class e2m2e.algorithm.manifold.manifolds.ManifoldTube(orbit, kind, branch, epsilon, trajectories=<factory>)[源代码]

基类:object

不变流形管(一族流形弧)。

变量:
参数:
orbit: Orbit
kind: ManifoldKind
branch: str
epsilon: float
trajectories: list[Orbit]
class e2m2e.algorithm.manifold.manifolds.InvariantManifold(orbit, kind, branch, epsilon)[源代码]

基类:object

周期轨道的不变流形

变量:
  • orbit -- 周期轨道(可只存首点,但 period 须已设置)

  • kind -- 流形类型(稳定/不稳定)

  • branch -- 扰动分支 "+" / "-"

  • epsilon -- 无量纲扰动幅度(位置方向长度,典型取 50 km / DU)

  • dynamics -- CR3BP_Dynamics 对象

参数:
SAMPLE_DT = 0.005
__init__(orbit, kind, branch, epsilon)[源代码]

初始化不变流形

参数:
  • orbit (Orbit) -- 周期轨道(period 须已知;只存首点时 seeds 会先传播一周采样相位)

  • kind (ManifoldKind) -- ManifoldKind.STABLE 或 ManifoldKind.UNSTABLE

  • branch (Literal['+', '-']) -- 扰动方向分支,"+""-"

  • epsilon (float) -- 无量纲扰动幅度(正值)

抛出:

ValueError -- 轨道周期未知、未关联 system、branch 或 epsilon 非法

返回类型:

None

seeds(n_points)[源代码]

生成相位扫掠种子

沿周期轨道均匀取 n_points 个相位点,把首点处的流形特征向量 用各段 STM 转运到该相位,位置部分归一化后施加 ±ε 扰动。

参数:

n_points (int) -- 相位点个数

返回:

种子状态数组,形状 (n_points, 6)

返回类型:

ndarray

propagate(t_span, section=None, n_workers=1)[源代码]

批量传播流形弧

积分方向由 kind 决定:STABLE 反向积分、UNSTABLE 正向积分, t_span 取绝对值作为积分时长。

参数:
  • t_span (float) -- 积分时长(无量纲时间,符号被忽略)

  • section (PoincareSection | None) -- 可选庞加莱截面;给定时每条弧在首次穿越截面处截断, 并把求精后的穿越态追加为弧的末点

  • n_workers (int) -- 并行 worker 数(>1 启用 Rayon)

返回:

流形管,含全部流形弧

返回类型:

ManifoldTube

庞加莱截面工具模块

提供平面截面与近拱点截面的定义,以及两类截面穿越检测:

  • 积分中检测PoincareSection.event() 生成事件函数 g(t, state) -> float (携带 direction/terminal 属性),传给 Dynamics.propagate(events=...),由积分器在步内定位穿越,无需密采样。

  • 事后检测crossings/detect_crossings,传播时密采样 t_eval → 逐采样点求截面函数 s(t)(平面:state[axis]-value;近拱点: r·v,r 为相对中心天体的位置)→ 符号变化区间内对线性插值态用二分 法求精,穿越态残差可达 1e-10 以下。

class e2m2e.algorithm.manifold.sections.SectionCrossings(section, states, times, trajectory_index)[源代码]

基类:object

截面穿越结果容器。

变量:
参数:
section: PoincareSection
states: ndarray
times: ndarray
trajectory_index: ndarray
e2m2e.algorithm.manifold.sections.detect_crossings(times, states, section_fn)[源代码]

事后截面穿越检测

逐采样点求截面函数,符号变化区间内对分段线性插值态用二分法求精。 正向与反向(时间递减)积分均适用。

参数:
  • times (ndarray) -- 采样时刻,形状 (n,),单调(递增或递减)

  • states (ndarray) -- 采样状态,形状 (n, 6)

  • section_fn (Callable[[ndarray], float]) -- 截面函数,输入 6 维状态返回标量,零点即截面

返回:

穿越列表,每项为 (穿越时刻, 穿越态, 区间左端采样索引)

返回类型:

list[tuple[float, ndarray, int]]

class e2m2e.algorithm.manifold.sections.PoincareSection(section_fn, description='')[源代码]

基类:object

庞加莱截面

截面由标量函数 s(state) 的零等值面定义。提供两类常用构造: 平面截面(某一坐标分量等于给定值)与近拱点截面(相对中心天体 位置 r 与速度 v 的点积为零,即 r·v = 0)。

参数:
  • section_fn (Callable[[np.ndarray], float])

  • description (str)

event(direction=0, terminal=False)[源代码]

生成事件函数,供积分中检测使用

返回的 callable g(t, state) -> float 可直接传给 Dynamics.propagate(events=...);函数对象携带 direction/terminal 属性。截面函数只依赖前 6 维物理状态, STM 增广传播(42 维状态)时自动截取前 6 维。

参数:
  • direction (int) -- 穿越方向过滤。> 0 只记上行穿越(s 由负到正), < 0 只记下行,0 双向

  • terminal (bool) -- True 时首次触发即终止积分

返回:

事件函数,零点即截面

返回类型:

Callable[[float, ndarray], float]

classmethod plane(axis, value)[源代码]

构造平面截面 s = state[axis] - value

参数:
  • axis (int) -- 状态分量索引(0=x, 1=y, 2=z, 3=vx, 4=vy, 5=vz)

  • value (float) -- 平面取值

返回:

PoincareSection 实例

返回类型:

PoincareSection

classmethod periapsis(center, system)[源代码]

构造近拱点截面 s = r·v

r 为相对 center 天体的位置。CR3BP 会合系中主天体固定于 x=-μ, 次天体固定于 x=1-μ。center 与 system 的主/次天体名称 (不区分大小写)匹配,"earth"/"moon" 分别回退为主/次天体。

参数:
  • center (str) -- 中心天体名称(如 "earth" / "moon")

  • system (CR3BP_System) -- CR3BP_System 对象(提供 mu 与天体名称)

返回:

PoincareSection 实例

返回类型:

PoincareSection

crossings(tube)[源代码]

检测流形管中所有流形弧的截面穿越

参数:

tube (ManifoldTube) -- 流形管(或任何带 trajectories 列表的容器, 每条轨迹须含 timesstates

返回:

全部穿越点(插值求精后)

返回类型:

SectionCrossings

流形拼接与低能转移流水线模块。

低能转移初猜生成(郑越、赵敏 2023 流程的产品化):

  1. 出发轨道不稳定流形与目标轨道稳定流形各自传播到同一庞加莱截面;

  2. patch_manifolds() 把两管穿越点交给默认 Rust 数值核完成两两配对、 位置/速度差和加权排序;

  3. design_low_energy_transfer()ThreeBodyLambert 打靶闭合弧段。

Python 侧保留流形对象与高层编排;backend="python" 只作显式等价性参照, Rust 不可用时不自动回退。

class e2m2e.algorithm.transfer.low_energy.PatchCandidate(i_a, i_b, state_a, state_b, delta_r, delta_v, cost)[源代码]

基类:object

流形管拼接候选(同一截面上两个穿越点的配对)。

变量:
  • i_a (int) -- 管 A 穿越点索引

  • i_b (int) -- 管 B 穿越点索引

  • state_a (numpy.ndarray) -- 管 A 穿越态(无量纲六维向量)

  • state_b (numpy.ndarray) -- 管 B 穿越态(无量纲)

  • delta_r (float) -- 位置差(无量纲)

  • delta_v (float) -- 速度差(无量纲)

  • cost (float) -- 加权拼接代价,w_r·位置差 + w_v·速度差

参数:
i_a: int
i_b: int
state_a: ndarray
state_b: ndarray
delta_r: float
delta_v: float
cost: float
e2m2e.algorithm.transfer.low_energy.patch_manifolds(m_a, m_b, section, weights=(1.0, 1.0), *, backend='rust', parallel=None, n_workers=None, progress_callback=None)[源代码]

两流形管在同一截面的穿越点两两配对,按拼接代价升序输出。

默认 backend='rust':配对、范数、代价和排序均由 Rust 执行;Rust 扩展缺失时直接报错。backend='python' 仅保留作显式等价性对照,绝不 作为运行时回退。截面对象和流形管遍历仍由 Python 编排层负责。

参数:
  • m_a (ManifoldTube) -- 管 A(如出发轨道不稳定流形)

  • m_b (ManifoldTube) -- 管 B(如目标轨道稳定流形)

  • section (PoincareSection) -- 庞加莱截面

  • weights (tuple[float, float]) -- (w_r, w_v) 位置/速度差权重(无量纲量纲下的相对权重)

  • backend (Literal['rust', 'python']) -- 'rust'``(默认)或显式参照 ``'python'

  • parallel (bool | None) -- Rust 路径是否使用 Rayon;None 由环境变量决定

  • n_workers (int | None) -- Rust Rayon 线程数;None 使用全局线程池

  • progress_callback (Callable[[int], Any] | None) -- 每完成一个截面态配对调用,Rust 侧可合并增量

返回:

cost 升序排列的候选列表;任一管无穿越时返回空列表

返回类型:

list[PatchCandidate]

e2m2e.algorithm.transfer.low_energy.design_low_energy_transfer(departure, target, epoch=None, model='cr3bp')[源代码]

低能转移流水线:流形拼接初猜 + CR3BP 打靶闭合。

流程:出发轨道不稳定流形 + 目标轨道稳定流形(± 分支都试)→ 次天体近拱点截面拼接取最优候选 → 出发弧(流形弧,精确)+ ThreeBodyLambert 把拼接点之后闭合到目标轨道。

脉冲构成:出发脉冲(上出发流形)+ 拼接脉冲(截面处)+ 到达脉冲(入目标流形)。

参数:
  • departure (OrbitTerminal) -- 出发轨道终端

  • target (Orbit) -- 目标周期轨道(须关联 system 且 period 已知)

  • epoch -- 参考历元(预留;星历转换未接入,见下)

  • model (Literal['cr3bp']) -- 动力学模型,当前仅支持 "cr3bp"

返回:

TransferSolution,两段弧(物理单位)

抛出:

ValueError -- model 不支持,或两流形管在截面上无穿越点

返回类型:

TransferSolution

TODO: 星历转换(CR3BP 闭合解 → 星历模型)未接入;接入时复用设计链路

e2m2e.algorithm.design 的 Rust 多重打靶修正,epoch 参数即为其入口。