e2m2e.algorithms package

e2m2e算法模块

包含用于轨道设计和优化的各种算法。

e2m2e.algorithms.compute_halo_coefficients(mu, L)[源代码]

计算 Halo 轨道 Richardson 三阶近似所需的全部系数。

根据 Richardson (1980) 的三阶解析构造方法,在共线平动点附近将 CR3BP 运动方程 展开为非线性扰动级数。三阶近似将轨道位移分解为面内(u-v)和面外(w)分量, 用 Fourier 级数表示,包含基频 omega_p 的各阶谐波。

核心系数包括: - a_ij: u 方向(沿主天体连线)的振幅修正系数 - b_ij: v 方向(面内垂直方向)的振幅修正系数 - d_ij: w 方向(面外方向)的振幅修正系数 - k, delta: 与平动点位置相关的符号因子 - kappa1, kappa2: 频率修正系数(用于计算非线性周期)

参数:
  • mu (float) -- 质量比 mu = m2 / (m1 + m2)。

  • L (int) -- 拉格朗日点编号(1=L1, 2=L2)。

返回:

包含所有 Richardson 三阶近似系数的字典,键包括:
  • gamma: 次天体到平动点的距离(L2 时取负值)

  • omega_p: 面内振荡基频

  • c1, c2, c3: Legendre 系数(有效势展开的前三阶)

  • a21~a31, b21~b31, d21~d32: 各方向振幅修正系数

  • k, delta: 符号因子

  • l1~l3, kappa1, kappa2: 频率和周期修正系数

返回类型:

Dict[str, float]

抛出:

ValueError -- 当 L 不为 1 或 2 时。

引用

Richardson, D. L. (1980). Analytic construction of periodic orbits about the collinear points. Celestial Mechanics, 22(3), 303-320.

e2m2e.algorithms.halo_third_order_approximation(mu, Au, Aw, phi, L, tf, N, halo_class=0)[源代码]

计算Halo轨道三阶解析近似

参数:
  • mu (float) -- 质量比

  • Au (float) -- U方向振幅

  • Aw (float) -- W方向振幅

  • phi (float) -- 相位偏移

  • L (int) -- 拉格朗日点 (1=L1, 2=L2)

  • tf (float) -- 终止时间

  • N (int) -- 点数

  • halo_class (int) -- 0=Class I (北), 1=Class II (南)

返回:

状态向量序列 (N, 6),[u, v, w, u_dot, v_dot, w_dot] t: 时间序列 T: 周期

返回类型:

SV_uvw

Reference:

Richardson, D. L. (1980). Analytic construction of periodic orbits about the collinear points. Celestial Mechanics.

e2m2e.algorithms.compute_halo_initial_guess(mu, z_amplitude, L=1, halo_class=0)[源代码]

计算Halo轨道初始猜测参数

使用 Richardson 三阶近似系数生成初始猜测,配合微分修正器使用。 初始状态位于 XZ 平面穿越点(y=0),赤道面穿越处(z=0)。

参数:
  • mu (float) -- 质量比

  • z_amplitude (float) -- Z方向振幅

  • L (int) -- 拉格朗日点 (1=L1, 2=L2)

  • halo_class (int) -- 0=北Halo, 1=南Halo

返回:

  • x0: 初始x坐标

  • y0: 初始y坐标 (0)

  • z0: 初始z坐标 (0)

  • vx0: 初始vx (0)

  • vy0: 初始vy

  • vz0: 初始vz (0)

  • T_half: 半周期

  • Au: U方向振幅

  • Aw: W方向振幅

返回类型:

包含初始猜测参数的字典

class e2m2e.algorithms.DifferentialCorrection(dynamic, target=None, free_vars=None)[源代码]

基类:object

微分修正算法

通过迭代修正初始条件,使轨道满足指定的约束条件(如周期性、对称性等)。

支持的对称性配置: - 2D对称X固定X0: 平面对称周期轨道,固定初始x坐标 - 2D对称X固定T: 平面对称周期轨道,固定轨道周期 - 3D对称X固定X0: 空间对称周期轨道(Halo轨道等) - 3D对称XZ固定X0: 空间XZ对称周期轨道 - 3D对称XZ固定Z0: 空间XZ对称周期轨道,固定Z0

参数:
dynamic

CR3BP_Dynamics对象

target_conditions

目标约束条件字典

free_variables

自由变量列表

tolerance

收敛容差

max_iterations

最大迭代次数

convergence_history

收敛历史

DEFAULT_TOLERANCE = 1e-12
DEFAULT_MAX_ITERATIONS = 50
DEFAULT_DAMPING_FACTOR = 1.0
VALID_SETUP_TYPES = ['2D_symmetric_x_fixed_x0', '2D_symmetric_x_fixed_t', '2D_symmetric_y_fixed_y0', '3D_symmetric_x_fixed_x0', '3D_symmetric_xz_fixed_x0', '3D_symmetric_xz_fixed_z0', 'halo_orbit_fixed_z0', 'halo_orbit_fixed_x0']
__init__(dynamic, target=None, free_vars=None)[源代码]

初始化修正器

参数:
  • dynamic (CR3BP_Dynamics) -- CR3BP_Dynamics对象

  • target (dict[str, Any] | None) -- 目标约束条件字典(可选)

  • free_vars (list[str] | None) -- 自由变量列表(可选)

返回类型:

None

setup_2D_symmetric_x_fixed_x0(x0=0.0)[源代码]

配置平面问题中固定初始x坐标的对称周期轨道搜索

在平面圆形限制性三体问题(PCRTBP)模型中,动力学方程关于会合坐标系的x轴具有对称性。 利用这一性质,周期轨道的搜索可以简化为寻找合适的初始条件: 从x轴上一点垂直出发(y=0, x_dot=0),经过半周期T/2后再次垂直穿越x轴(y=0, x_dot=0)。

本函数针对这种对称性设置微分修正问题,固定初始x坐标x0,将初始y方向速度y_dot0 和半周期T/2作为自由变量进行调整,以满足终点处的垂直穿越条件。

参数:

x0 -- 固定的初始x坐标,轨道从点(x0, 0)垂直出发, 这里将值设置为0.0是因为在这个函数中只需要使用x0进行初始化

返回:

配置好的微分修正器实例

备注

  • 自由变量: [y_dot0, T_half] - 初始y方向速度和半周期时间

  • 目标约束: [y(T/2)=0, x_dot(T/2)=0] - 终点处再次垂直穿越x轴

  • 状态向量索引: [1, 3] 分别对应y坐标和x方向速度

此配置对应于Broucke(1968)等经典文献中寻找对称周期轨道的基本方法, 可用于生成围绕平动点或主天体的各类周期轨道家族。

Reference:

Broucke R A. Periodic orbits in the restricted three body problem with Earth-moon masses[R]. 1968.

setup_2D_symmetric_x_fixed_t(t_half)[源代码]

配置平面问题中固定半周期的对称周期轨道搜索

固定半周期T/2,调整初始条件x0和y_dot0满足约束。

参数:

t_half -- 固定的半周期

返回:

配置好的微分修正器实例

setup_2D_symmetric_y_fixed_y0(y0=0.0)[源代码]

配置平面问题中固定初始y坐标的y轴对称周期轨道搜索

适用于共振轨道(RO)等从y轴出发的周期轨道。 轨道从点(0, y0)出发(x=0, x_dot=0),经过半周期T/2后再次穿越y轴(x=0, x_dot=0)。

参数:

y0 -- 固定的初始y坐标

返回:

配置好的微分修正器实例

备注

  • 自由变量: [x_dot0, T_half] - 初始x方向速度和半周期时间

  • 目标约束: [x(T/2)=0, x_dot(T/2)=0] - 终点处再次穿越y轴

  • 状态向量索引: [0, 3] 分别对应x坐标和x方向速度

setup_3D_symmetric_x_fixed_x0(x0)[源代码]

配置空间问题中固定初始x坐标的对称周期轨道搜索(如Halo轨道)

参数:

x0 (float) -- 固定的初始x坐标

返回:

配置好的微分修正器实例

返回类型:

self

setup_3D_symmetric_xz_fixed_x0(x0)[源代码]

配置空间XZ对称周期轨道搜索,固定X0

参数:

x0 (float) -- 固定的初始x坐标

返回:

配置好的微分修正器实例

返回类型:

self

setup_3D_symmetric_xz_fixed_z0(z0)[源代码]

配置空间XZ对称周期轨道搜索,固定Z0

参数:

z0 (float) -- 固定的初始z坐标

返回:

配置好的微分修正器实例

返回类型:

self

setup_halo_orbit_fixed_z0(z0, libration_point=1)[源代码]

配置 Halo 轨道微分修正,固定初始 Z0(XZ 对称)

Halo 轨道具有 XZ 平面对称性,利用该对称性可以将问题简化为: 从 XZ 平面上一点 (x0, 0, z0) 出发,经过半周期 T/2 后再次到达 XZ 平面。

参数:
  • z0 (float) -- 固定的初始 z 坐标

  • libration_point (int) -- 平动点编号 (1=L1, 2=L2),默认 L1

返回:

配置好的微分修正器实例

返回类型:

self

备注

  • 自由变量: [x0, y_dot0, T_half] - 初始 x 坐标、y 方向速度和半周期时间

  • 目标约束: [y(T/2)=0, x_dot(T/2)=0, z_dot(T/2)=0] - 半周期处再次位于 XZ 平面

  • 状态向量索引: [0, 4, 6] 分别对应 x0、y_dot0 和时间 T_half

  • 注意: z 在半周期时会改变符号,不作为约束

setup_halo_orbit_fixed_x0(x0, libration_point=1)[源代码]

配置 Halo 轨道微分修正,固定初始 X0(XZ 对称)

Halo 轨道具有 XZ 平面对称性,利用该对称性可以将问题简化为: 从 XZ 平面上一点 (x0, 0, z0) 出发,经过半周期 T/2 后再次到达 XZ 平面。

参数:
  • x0 (float) -- 固定的初始 x 坐标

  • libration_point (int) -- 平动点编号 (1=L1, 2=L2),默认 L1

返回:

配置好的微分修正器实例

返回类型:

self

备注

  • 自由变量: [z0, y_dot0, T_half] - 初始 z 坐标、y 方向速度和半周期时间

  • 目标约束: [y(T/2)=0, x_dot(T/2)=0, z_dot(T/2)=0] - 半周期处再次位于 XZ 平面且垂直穿越

  • 状态向量索引: [2, 4, 6] 分别对应 z0、y_dot0 和时间 T_half

iterate_correction(initial_guess, verbose=False)[源代码]

迭代修正主算法(基于STM的牛顿法)

通过状态转移矩阵(STM)构建雅可比矩阵,使用牛顿迭代法修正自由变量, 使终点状态满足目标约束条件,从而找到精确的周期轨道。

参数:
  • initial_guess (Orbit) -- 初始猜测轨道,或初始状态向量

  • verbose (bool) -- 是否打印迭代过程信息

返回:

  • 返回修正后的 Orbit 对象

check_convergence()[源代码]

检查收敛性

返回:

是否收敛

返回类型:

bool

get_convergence_history()[源代码]

获取收敛历史

返回:

收敛历史数据

返回类型:

dict

class e2m2e.algorithms.Continuation(corrector, step=None)[源代码]

基类:object

轨道族延拓

通过延拓算法生成一族周期轨道,支持自然参数延拓和伪弧长延拓。

参数:
corrector

DifferentialCorrection对象

continuation_parameter

延拓参数名称

step_size

延拓步长

family_orbits

轨道族列表

DEFAULT_STEP_SIZE = 0.01
MIN_STEP_SIZE = 1e-06
MAX_STEP_SIZE = 0.1
DEFAULT_PREDICTOR_ORDER = 1
__init__(corrector, step=None)[源代码]

初始化延拓器

Args: - corrector: DifferentialCorrection对象 - param: 延拓参数(如 "energy", "period", "amplitude", "x0", "z0") - step: 初始步长

参数:
返回类型:

None

natural_continuation(seed_orbit, param_range, step_size, verbose=False)[源代码]

自然参数延拓

从种子轨道出发,逐步改变延拓参数,生成一族周期轨道。 支持双向延拓:如果param_range的最小值小于种子轨道参数值,则向小值方向延拓; 如果param_range的最大值大于种子轨道参数值,则向大值方向延拓。

参数:
  • seed_orbit -- Orbit, 种子轨道

  • param_range -- tuple, 延拓参数范围 (param_min, param_max)

  • step_size -- float, 步长(始终为正值,延拓方向由参数范围自动确定)

  • verbose -- 是否打印信息

返回:

包含轨道族的OrbitFamily对象

返回类型:

OrbitFamily

pseudo_arclength_continuation(seed_orbit, n_orbits=50, step_size=0.005, direction='positive', verbose=True, TolPAL=1e-06, TolDiffCorr=1e-06, IterMax=100, dc_scheme='adaptive', libration_point=1, directional_increment=False, target_vector=0, target_direction=1)[源代码]

伪弧长延拓(对应 MATLAB continuation_PAL_CR3BP,plane=13 / XZ 对称)

自由变量 X = [rx, rz, vy, T/2]Xdot = null(dF),PAL 约束 G = [F; (Xnew-X)·Xdot - DeltaS];内层用 Xnew 计算 F``(与 MATLAB 仅用固定 ``X 相比更一致)。每步后用微分修正闭合。

参数:
  • seed_orbit -- 种子轨道(Orbit对象)

  • n_orbits (int) -- 本支要生成的新轨道条数(与 MATLAB 中 N 一致)

  • step_size (float) -- 伪弧长步长 |DeltaS|``(正数;``direction 决定符号)

  • direction (str) -- positivenegative``(双侧延拓请调用两次或使用 ``halo_pseudo_arclength_continuation(..., direction='both')

  • dc_scheme (str) -- 见 Continuation 文档字符串

  • target_vector (int) -- 与 MATLAB TargetVector 对应的 0 基下标:0=rx,1=rz,2=vy,3=T/2

  • verbose (bool)

  • TolPAL (float)

  • TolDiffCorr (float)

  • IterMax (int)

  • libration_point (int)

  • directional_increment (bool)

  • target_direction (int)

返回:

仅含种子 + 本支新轨道(不重复添加种子)

返回类型:

OrbitFamily

generate_halo_seed_orbit(libration_point, amplitude_z, halo_class=0, verbose=False)[源代码]

生成Halo轨道初始猜测并修正

使用Richardson三阶近似生成初始猜测,然后通过微分修正 得到精确的周期轨道。

参数:
  • libration_point (int) -- 拉格朗日点 (1=L1, 2=L2)

  • amplitude_z (float) -- Z方向振幅

  • halo_class (int) -- 0=北Halo, 1=南Halo

  • verbose (bool) -- 是否打印详细信息

返回:

Halo周期轨道

返回类型:

Orbit

generate_halo_family(seed_orbit, n_orbits=50, direction='positive', step_size=0.001)[源代码]

生成Halo轨道族

使用自然参数延拓法生成Halo轨道族。

参数:
  • seed_orbit (Orbit) -- 种子轨道

  • n_orbits (int) -- 目标轨道数量

  • direction (str) -- 延拓方向 ("positive", "negative", "both")

  • step_size (float) -- 步长

返回:

Halo轨道族

返回类型:

List[Orbit]

halo_pseudo_arclength_continuation(seed_orbit, n_orbits=50, direction='both', step_size=0.0045, step_size_negative=None, verbose=True, TolPAL=1e-06, TolDiffCorr=1e-06, IterMax=100, dc_scheme='adaptive', directional_increment=True)[源代码]

Halo 轨道族伪弧长延拓(对齐 CR3BP_MATLAB_Library

对应 continuation_PAL_CR3BP + XZ 对称(X = [rx,rz,vy,T/2]),与 examples/FAMILY_L1Halo_North.m 的 PAL 步一致。

微分修正:MATLAB 脚本在 PAL 后使用 type=1``(固定 ``x0)。当前 Python 在 PAL 初值下 setup_halo_orbit_fixed_x0 易与 STM 牛顿耦合到非物理解,故默认 dc_scheme='adaptive'``(按 Δx/Δz fixed x / fixed z 间切换,与原 ``pseudo_arclength_continuation 一致)。若需对齐 MATLAB 的 fixedX,可设 matlab_halo_type1``(失败时会自动再试 ``fixed_z0)。 - DirectionalIncrementTargetVectorTargetDirection 与脚本一致:

正向支 TargetVector=2``(``rz)、TargetDirection=+1;负向支 TargetVector=1``(``rx)、TargetDirection=-1``(0 基下标见 ``pseudo_arclength_continuation)。

FAMILY_L1Halo_North.m 中正负支 DeltaS 不同(0.0045 与 0.009),可用 step_size_negative 单独指定负向步长模长。

参数:
  • seed_orbit (Orbit) -- 已收敛的 Halo 种子轨道

  • n_orbits (int) -- 每一支的新轨道条数 N``(与 MATLAB 一致;``direction=both 时 正向、负向各生成 n_orbits 条)

  • step_size (float) -- 正向伪弧长步长 ``|DeltaS|``(默认 0.0045)

  • step_size_negative (float | None) -- 负向步长模长,默认与 step_size 相同

  • dc_scheme (str) -- matlab_halo_type1 / matlab_halo_type2 / adaptive

  • direction (str)

  • verbose (bool)

  • TolPAL (float)

  • TolDiffCorr (float)

  • IterMax (int)

  • directional_increment (bool)

返回:

种子 + 各支新轨道(无重复种子)

返回类型:

OrbitFamily

class e2m2e.algorithms.StabilityAnalysis(orbit, dynamics=None)[源代码]

基类:object

轨道稳定性分析

计算轨道的单值矩阵、Floquet乘子、稳定性指数等, 并进行稳定性分类和分岔检测。

参数:
orbit

Orbit对象

dynamic

CR3BP_Dynamics对象

monodromy_matrix

单值矩阵

eigenvalues

特征值

stability_indices

稳定性指数

STABILITY_THRESHOLD = 1e-06
BIFURCATION_TOLERANCE = 1e-08
__init__(orbit, dynamics=None)[源代码]

初始化分析器

参数:
  • orbit (Orbit) -- Orbit对象

  • dynamics (CR3BP_Dynamics | None) -- CR3BP_Dynamics对象(可选,如果orbit关联了system则自动创建)

返回类型:

None

compute_monodromy()[源代码]

计算单值矩阵

通过积分一个完整周期的状态转移矩阵获得单值矩阵。

返回:

6x6 单值矩阵

返回类型:

np.ndarray

compute_floquet_multipliers()[源代码]

计算Floquet乘子(特征值)

返回:

Floquet乘子

返回类型:

np.ndarray

compute_stability_index()[源代码]

计算稳定性指数

Broucke稳定性参数定义为:ν = λ + 1/λ, 其中λ是单值矩阵的特征值。

返回:

稳定性指数字典

返回类型:

dict

classify_orbit()[源代码]

对轨道进行稳定性分类

返回:

稳定性分类结果

返回类型:

dict

analyze_bifurcation()[源代码]

分析分岔类型

通过检查单值矩阵特征值的分布判断是否存在分岔。

返回:

分岔分析结果

返回类型:

dict

full_analysis()[源代码]

执行完整的稳定性分析

返回:

完整分析结果

返回类型:

dict

static detect_bifurcation_in_family(orbits, dynamics, tolerance=1e-08)[源代码]

检测轨道族中的分岔点

遍历轨道族中的每条轨道,计算其单值矩阵特征值, 检测是否有特征值接近 +1(切分岔/saddle-node bifurcation)。

参数:
  • orbits (list[Orbit]) -- Orbit对象列表(轨道族)

  • dynamics (CR3BP_Dynamics) -- CR3BP_Dynamics对象

  • tolerance (float) -- 特征值接近 +1 的容差,默认 1e-8

返回:

分岔点列表,每个元素包含:
  • orbit_index: 轨道在族中的索引

  • orbit: Orbit对象

  • eigenvalues: 特征值数组

  • eigenvalue_diff: |λ - 1| 的最小值

  • bifurcation_type: 分岔类型

返回类型:

List[Dict[str, Any]]

static find_nearest_bifurcation(orbits, dynamics, target_x0=None, tolerance=0.0001)[源代码]

在轨道族中找到最接近目标参数的分岔点

参数:
  • orbits (list[Orbit]) -- Orbit对象列表

  • dynamics (CR3BP_Dynamics) -- CR3BP_Dynamics对象

  • target_x0 (float | None) -- 目标x0坐标(可选)

  • tolerance (float) -- 搜索容差

返回:

分岔点字典,如果未找到则返回None

返回类型:

dict[str, Any] | None

class e2m2e.algorithms.StabilityType(*values)[源代码]

基类:Enum

稳定性类型枚举

STABLE = 'stable'
UNSTABLE = 'unstable'
MARGINALLY_STABLE = 'marginally_stable'
HYPERBOLIC = 'hyperbolic'
ELLIPTIC = 'elliptic'
PARABOLIC = 'parabolic'
class e2m2e.algorithms.BifurcationType(*values)[源代码]

基类:Enum

分岔类型枚举

NONE = 'none'
PERIOD_DOUBLING = 'period_doubling'
SADDLE_NODE = 'saddle_node'
TORUS = 'torus'
PITCHFORK = 'pitchfork'
TRANSCRITICAL = 'transcritical'
SECONDARY_HOPF = 'secondary_hopf'
class e2m2e.algorithms.MultipleShooting(dynamics, n_workers=1, kernel_dir=None)[源代码]

基类:object

多重打靶法(Multiple Shooting)修正器。

将一条轨迹分为 N 个节点、n_seg = N-1 段弧段,对每段独立积分后, 通过匹配相邻段端点状态来构建残差向量,再利用雅可比矩阵(含 STM) 进行最小二乘修正,反复迭代直到残差满足容差。

当 var_time=True 时,时间节点也作为自由变量参与修正(适用于自由时间问题)。

并行策略

  • n_workers=1:串行(默认)。

  • n_workers>1, kernel_dir=None:多线程(ThreadPoolExecutor)。 适合 CR3BP 等纯 Python/NumPy 动力学,但受 GIL 限制,并发收益有限。

  • n_workers>1, kernel_dir=<路径>**多进程**(ProcessPoolExecutor)。 每个子进程重载 SPICE 内核,绕过 GIL,可充分利用多核 CPU。 仅适用于 ``EphemerisDynamics``(需 SPICE 内核)。

__init__(dynamics, n_workers=1, kernel_dir=None)[源代码]
参数:
  • dynamics --

    动力学模型对象,需提供以下接口: - propagate(state, time_span, with_stm=True):

    积分传播,返回含 "states" 和 "stm" 的字典

    • equations_of_motion(t, state): 计算状态导数(右端函数值)

  • n_workers (int) --

    并行工作进程/线程数,默认 1(串行)。 - n_workers>1kernel_dir 已设置:

    使用 ProcessPoolExecutor(多进程,推荐)。

    • n_workers>1kernel_dir=None: 使用 ThreadPoolExecutor(多线程,受 GIL 限制)。

  • kernel_dir (str | None) -- SPICE 内核目录路径 (含 de440.bspnaif0012.tls)。 仅在 n_workers>1 时需要。设置后自动启用多进程模式。 例:kernel_dir=os.environ.get("SPICE_KERNEL_DIR", "../e2m2e/kernels")

返回类型:

None

correct(t_patch, state_patch, var_time=False, max_iter=None, tolerance=None, verbose=False)[源代码]

执行多重打靶修正。

将整条轨迹分为若干弧段,对每段独立积分后检验节点处的状态连续性, 利用状态转移矩阵(STM)组装雅可比矩阵,通过最小二乘求解修正量并迭代。

参数:
  • t_patch (ndarray) -- 初始时间节点数组,长度 N

  • state_patch (ndarray) -- 初始状态量数组,形状 (N, 6),每行 [x, y, z, vx, vy, vz]

  • var_time (bool) -- 是否允许时间节点作为自由变量参与修正

  • max_iter (int | None) -- 最大迭代次数(默认使用 self.max_iter)

  • tolerance (float | None) -- 收敛容差(默认使用 self.tolerance)

  • verbose (bool) -- 是否显示进度条

返回:

包含修正后的时间/状态、收敛标志、迭代次数和残差历史

返回类型:

MultipleShootingResult

参数:
  • n_workers (int)

  • kernel_dir (str | None)

e2m2e.algorithms.sample_patch_points(orbit, n_points)[源代码]

沿周期轨道均匀采样 patch points(打靶节点)。

该方法用于多重打靶法(Multiple Shooting)的前处理,将一条周期轨道 在时间上均匀分割为 n_points 个节点,并通过线性插值获取每个节点处的状态。

参数:
  • orbit -- 轨道对象,需包含以下属性: - period: 轨道周期(归一化时间单位) - times: 时间数组,形状 (M,) - states: 状态数组,形状 (M, 6),每行 [x, y, z, vx, vy, vz]

  • n_points (int) -- 需要采样的节点数量

返回:

包含两个数组的元组:
  • t_patch: 采样时间节点数组,形状 (n_points,),归一化时间单位

  • states: 采样状态数组,形状 (n_points, 6),每行 [x, y, z, vx, vy, vz]

返回类型:

Tuple[np.ndarray, np.ndarray]

抛出:

ValueError -- 当轨道对象没有 period 属性时抛出

备注

  • 采样时间范围为 [0, period),不包含周期终点(endpoint=False)

  • 使用线性插值从原始轨道数据中获取节点状态

  • 适用于 CR3BP 归一化坐标系下的周期轨道采样

e2m2e.algorithms.convert_to_j2000(t_patch_syn, states_syn, syn_j2000, reference_et, tu_days=4.34811305)[源代码]

将 synodic 坐标系下的 patch points 转换到 J2000 惯性坐标系。

该方法用于将 CR3BP 归一化 synodic 坐标系中的轨道节点转换到 J2000 惯性坐标系,以便在星历模型(Ephemeris)中进行高精度轨道修正。

参数:
  • t_patch_syn (ArrayLike) -- synodic 坐标系下的时间节点数组,归一化时间单位(TU)

  • states_syn (ArrayLike) -- synodic 坐标系下的状态数组,形状 (N, 6), 每行 [x, y, z, vx, vy, vz],归一化单位(DU, DU/TU)

  • syn_j2000 -- SynodicJ2000Transformation 对象,提供坐标转换功能

  • reference_et (float) -- 参考历元的 SPICE ephemeris time(ET),单位秒

  • tu_days (float) -- 归一化时间单位(TU)对应的天数,默认值为 4.34811305 天

返回:

包含两个数组的元组:
  • t_patch_j2000: J2000 坐标系下的时间数组,SPICE ET(秒)

  • states_j2000: J2000 坐标系下的状态数组,形状 (N, 6),

    每行 [x, y, z, vx, vy, vz],单位(km, km/s)

返回类型:

Tuple[np.ndarray, np.ndarray]

备注

  • 时间转换公式:t_j2000 = reference_et + t_syn * (tu_days * 86400)

  • 状态转换使用 SynodicJ2000Transformation.batch_synodic_to_j2000() 方法

  • 适用于将 CR3BP 轨道转换到星历模型进行高精度修正的场景

  • 转换后的状态可用于 EphemerisDynamics 进行轨道传播

class e2m2e.algorithms.CorrectionConfig(setup_type, symmetry_condition, fixed_parameters=<factory>, free_variables=<factory>, free_variable_indices=<factory>, target_conditions=<factory>, constraint_indices=<factory>, constraint_weights=<factory>, constraint_types=<factory>)[源代码]

基类:object

Immutable configuration for a differential correction strategy.

Encapsulates all correction parameters that were previously scattered across individual setup_* method bodies in DifferentialCorrection.

参数:
setup_type

Identifier string for the correction setup type.

Type:

str

symmetry_condition

Symmetry exploited by the correction (e.g. 'x_axis').

Type:

str

fixed_parameters

Parameter values held constant during correction.

Type:

dict[str, float]

free_variables

Names of variables the Newton solver adjusts.

Type:

list[str]

free_variable_indices

State-vector indices corresponding to free variables.

Type:

list[int]

target_conditions

Constraint names mapped to their target values.

Type:

dict[str, float]

constraint_indices

State-vector indices for constraint evaluation.

Type:

list[int]

constraint_weights

Per-constraint weighting factors for the Jacobian.

Type:

dict[str, float]

constraint_types

Per-constraint classification (e.g. 'equality').

Type:

dict[str, str]

setup_type: str
symmetry_condition: str
fixed_parameters: dict[str, float]
free_variables: list[str]
free_variable_indices: list[int]
target_conditions: dict[str, float]
constraint_indices: list[int]
constraint_weights: dict[str, float]
constraint_types: dict[str, str]

Submodules

e2m2e.algorithms.continuation module

轨道族延拓算法模块

提供自然参数延拓和伪弧长延拓方法,用于生成轨道族。 包括Halo轨道、Lyapunov轨道等的生成功能。

e2m2e.algorithms.continuation.compute_F_and_dF_symmetric_xz_plane(X, SV0, dynamics)[源代码]

计算XZ平面对称轨道的约束向量和雅可比矩阵

对应MATLAB: computeFdF_symPeriodicPlanes_CR3BP(X, SV0i, mu, plane=13)

X = [rx; rz; vy; tf2] - 自由变量向量 SV0 = [rx, ry, rz, vx, vy, vz] - 初始状态向量

F = [vx; vz; ry] - 约束向量(半周期终点状态) dF = ∂F/∂X - 约束雅可比矩阵 (3x4)

参数:
  • X (ndarray) -- 自由变量向量 [rx, rz, vy, tf2]

  • SV0 (ndarray) -- 初始状态向量 [x, y, z, vx, vy, vz]

  • dynamics (CR3BP_Dynamics) -- CR3BP_Dynamics实例,提供运动方程和雅可比矩阵

返回:

约束向量 [vx, vz, ry] dF: 雅可比矩阵 (3, 4)

返回类型:

F

e2m2e.algorithms.continuation.compute_tangent_vector(dF)[源代码]

计算切向量(约束雅可比矩阵的零空间)

对应MATLAB: Xdot = null(dF)

参数:

dF (ndarray) -- 约束雅可比矩阵 (3, 4)

返回:

切向量 (4,),单位化

返回类型:

Xdot

class e2m2e.algorithms.continuation.Continuation(corrector, step=None)[源代码]

基类:object

轨道族延拓

通过延拓算法生成一族周期轨道,支持自然参数延拓和伪弧长延拓。

参数:
corrector

DifferentialCorrection对象

continuation_parameter

延拓参数名称

Type:

str | None

step_size

延拓步长

family_orbits

轨道族列表

Type:

list[e2m2e.core.orbit.Orbit]

DEFAULT_STEP_SIZE = 0.01
MIN_STEP_SIZE = 1e-06
MAX_STEP_SIZE = 0.1
DEFAULT_PREDICTOR_ORDER = 1
__init__(corrector, step=None)[源代码]

初始化延拓器

Args: - corrector: DifferentialCorrection对象 - param: 延拓参数(如 "energy", "period", "amplitude", "x0", "z0") - step: 初始步长

参数:
返回类型:

None

continuation_parameter: str | None
family_orbits: list[Orbit]
family_parameters: list[float]
family_states: list[ndarray]
family_periods: list[float]
natural_continuation(seed_orbit, param_range, step_size, verbose=False)[源代码]

自然参数延拓

从种子轨道出发,逐步改变延拓参数,生成一族周期轨道。 支持双向延拓:如果param_range的最小值小于种子轨道参数值,则向小值方向延拓; 如果param_range的最大值大于种子轨道参数值,则向大值方向延拓。

参数:
  • seed_orbit -- Orbit, 种子轨道

  • param_range -- tuple, 延拓参数范围 (param_min, param_max)

  • step_size -- float, 步长(始终为正值,延拓方向由参数范围自动确定)

  • verbose -- 是否打印信息

返回:

包含轨道族的OrbitFamily对象

返回类型:

OrbitFamily

pseudo_arclength_continuation(seed_orbit, n_orbits=50, step_size=0.005, direction='positive', verbose=True, TolPAL=1e-06, TolDiffCorr=1e-06, IterMax=100, dc_scheme='adaptive', libration_point=1, directional_increment=False, target_vector=0, target_direction=1)[源代码]

伪弧长延拓(对应 MATLAB continuation_PAL_CR3BP,plane=13 / XZ 对称)

自由变量 X = [rx, rz, vy, T/2]Xdot = null(dF),PAL 约束 G = [F; (Xnew-X)·Xdot - DeltaS];内层用 Xnew 计算 F``(与 MATLAB 仅用固定 ``X 相比更一致)。每步后用微分修正闭合。

参数:
  • seed_orbit -- 种子轨道(Orbit对象)

  • n_orbits (int) -- 本支要生成的新轨道条数(与 MATLAB 中 N 一致)

  • step_size (float) -- 伪弧长步长 |DeltaS|``(正数;``direction 决定符号)

  • direction (str) -- positivenegative``(双侧延拓请调用两次或使用 ``halo_pseudo_arclength_continuation(..., direction='both')

  • dc_scheme (str) -- 见 Continuation 文档字符串

  • target_vector (int) -- 与 MATLAB TargetVector 对应的 0 基下标:0=rx,1=rz,2=vy,3=T/2

  • verbose (bool)

  • TolPAL (float)

  • TolDiffCorr (float)

  • IterMax (int)

  • libration_point (int)

  • directional_increment (bool)

  • target_direction (int)

返回:

仅含种子 + 本支新轨道(不重复添加种子)

返回类型:

OrbitFamily

generate_halo_seed_orbit(libration_point, amplitude_z, halo_class=0, verbose=False)[源代码]

生成Halo轨道初始猜测并修正

使用Richardson三阶近似生成初始猜测,然后通过微分修正 得到精确的周期轨道。

参数:
  • libration_point (int) -- 拉格朗日点 (1=L1, 2=L2)

  • amplitude_z (float) -- Z方向振幅

  • halo_class (int) -- 0=北Halo, 1=南Halo

  • verbose (bool) -- 是否打印详细信息

返回:

Halo周期轨道

返回类型:

Orbit

generate_halo_family(seed_orbit, n_orbits=50, direction='positive', step_size=0.001)[源代码]

生成Halo轨道族

使用自然参数延拓法生成Halo轨道族。

参数:
  • seed_orbit (Orbit) -- 种子轨道

  • n_orbits (int) -- 目标轨道数量

  • direction (str) -- 延拓方向 ("positive", "negative", "both")

  • step_size (float) -- 步长

返回:

Halo轨道族

返回类型:

List[Orbit]

halo_pseudo_arclength_continuation(seed_orbit, n_orbits=50, direction='both', step_size=0.0045, step_size_negative=None, verbose=True, TolPAL=1e-06, TolDiffCorr=1e-06, IterMax=100, dc_scheme='adaptive', directional_increment=True)[源代码]

Halo 轨道族伪弧长延拓(对齐 CR3BP_MATLAB_Library

对应 continuation_PAL_CR3BP + XZ 对称(X = [rx,rz,vy,T/2]),与 examples/FAMILY_L1Halo_North.m 的 PAL 步一致。

微分修正:MATLAB 脚本在 PAL 后使用 type=1``(固定 ``x0)。当前 Python 在 PAL 初值下 setup_halo_orbit_fixed_x0 易与 STM 牛顿耦合到非物理解,故默认 dc_scheme='adaptive'``(按 Δx/Δz fixed x / fixed z 间切换,与原 ``pseudo_arclength_continuation 一致)。若需对齐 MATLAB 的 fixedX,可设 matlab_halo_type1``(失败时会自动再试 ``fixed_z0)。 - DirectionalIncrementTargetVectorTargetDirection 与脚本一致:

正向支 TargetVector=2``(``rz)、TargetDirection=+1;负向支 TargetVector=1``(``rx)、TargetDirection=-1``(0 基下标见 ``pseudo_arclength_continuation)。

FAMILY_L1Halo_North.m 中正负支 DeltaS 不同(0.0045 与 0.009),可用 step_size_negative 单独指定负向步长模长。

参数:
  • seed_orbit (Orbit) -- 已收敛的 Halo 种子轨道

  • n_orbits (int) -- 每一支的新轨道条数 N``(与 MATLAB 一致;``direction=both 时 正向、负向各生成 n_orbits 条)

  • step_size (float) -- 正向伪弧长步长 ``|DeltaS|``(默认 0.0045)

  • step_size_negative (float | None) -- 负向步长模长,默认与 step_size 相同

  • dc_scheme (str) -- matlab_halo_type1 / matlab_halo_type2 / adaptive

  • direction (str)

  • verbose (bool)

  • TolPAL (float)

  • TolDiffCorr (float)

  • IterMax (int)

  • directional_increment (bool)

返回:

种子 + 各支新轨道(无重复种子)

返回类型:

OrbitFamily

e2m2e.algorithms.differential_correction module

微分修正算法模块

提供用于求解周期轨道的微分修正算法,支持多种对称性配置。 包括Halo轨道、Lyapunov轨道等的解析近似和微分修正功能。

e2m2e.algorithms.differential_correction.compute_halo_coefficients(mu, L)[源代码]

计算 Halo 轨道 Richardson 三阶近似所需的全部系数。

根据 Richardson (1980) 的三阶解析构造方法,在共线平动点附近将 CR3BP 运动方程 展开为非线性扰动级数。三阶近似将轨道位移分解为面内(u-v)和面外(w)分量, 用 Fourier 级数表示,包含基频 omega_p 的各阶谐波。

核心系数包括: - a_ij: u 方向(沿主天体连线)的振幅修正系数 - b_ij: v 方向(面内垂直方向)的振幅修正系数 - d_ij: w 方向(面外方向)的振幅修正系数 - k, delta: 与平动点位置相关的符号因子 - kappa1, kappa2: 频率修正系数(用于计算非线性周期)

参数:
  • mu (float) -- 质量比 mu = m2 / (m1 + m2)。

  • L (int) -- 拉格朗日点编号(1=L1, 2=L2)。

返回:

包含所有 Richardson 三阶近似系数的字典,键包括:
  • gamma: 次天体到平动点的距离(L2 时取负值)

  • omega_p: 面内振荡基频

  • c1, c2, c3: Legendre 系数(有效势展开的前三阶)

  • a21~a31, b21~b31, d21~d32: 各方向振幅修正系数

  • k, delta: 符号因子

  • l1~l3, kappa1, kappa2: 频率和周期修正系数

返回类型:

Dict[str, float]

抛出:

ValueError -- 当 L 不为 1 或 2 时。

引用

Richardson, D. L. (1980). Analytic construction of periodic orbits about the collinear points. Celestial Mechanics, 22(3), 303-320.

e2m2e.algorithms.differential_correction.halo_third_order_approximation(mu, Au, Aw, phi, L, tf, N, halo_class=0)[源代码]

计算Halo轨道三阶解析近似

参数:
  • mu (float) -- 质量比

  • Au (float) -- U方向振幅

  • Aw (float) -- W方向振幅

  • phi (float) -- 相位偏移

  • L (int) -- 拉格朗日点 (1=L1, 2=L2)

  • tf (float) -- 终止时间

  • N (int) -- 点数

  • halo_class (int) -- 0=Class I (北), 1=Class II (南)

返回:

状态向量序列 (N, 6),[u, v, w, u_dot, v_dot, w_dot] t: 时间序列 T: 周期

返回类型:

SV_uvw

Reference:

Richardson, D. L. (1980). Analytic construction of periodic orbits about the collinear points. Celestial Mechanics.

e2m2e.algorithms.differential_correction.compute_halo_initial_guess(mu, z_amplitude, L=1, halo_class=0)[源代码]

计算Halo轨道初始猜测参数

使用 Richardson 三阶近似系数生成初始猜测,配合微分修正器使用。 初始状态位于 XZ 平面穿越点(y=0),赤道面穿越处(z=0)。

参数:
  • mu (float) -- 质量比

  • z_amplitude (float) -- Z方向振幅

  • L (int) -- 拉格朗日点 (1=L1, 2=L2)

  • halo_class (int) -- 0=北Halo, 1=南Halo

返回:

  • x0: 初始x坐标

  • y0: 初始y坐标 (0)

  • z0: 初始z坐标 (0)

  • vx0: 初始vx (0)

  • vy0: 初始vy

  • vz0: 初始vz (0)

  • T_half: 半周期

  • Au: U方向振幅

  • Aw: W方向振幅

返回类型:

包含初始猜测参数的字典

class e2m2e.algorithms.differential_correction.DifferentialCorrection(dynamic, target=None, free_vars=None)[源代码]

基类:object

微分修正算法

通过迭代修正初始条件,使轨道满足指定的约束条件(如周期性、对称性等)。

支持的对称性配置: - 2D对称X固定X0: 平面对称周期轨道,固定初始x坐标 - 2D对称X固定T: 平面对称周期轨道,固定轨道周期 - 3D对称X固定X0: 空间对称周期轨道(Halo轨道等) - 3D对称XZ固定X0: 空间XZ对称周期轨道 - 3D对称XZ固定Z0: 空间XZ对称周期轨道,固定Z0

参数:
dynamic

CR3BP_Dynamics对象

target_conditions

目标约束条件字典

free_variables

自由变量列表

tolerance

收敛容差

max_iterations

最大迭代次数

convergence_history

收敛历史

Type:

list[float]

DEFAULT_TOLERANCE = 1e-12
DEFAULT_MAX_ITERATIONS = 50
DEFAULT_DAMPING_FACTOR = 1.0
VALID_SETUP_TYPES = ['2D_symmetric_x_fixed_x0', '2D_symmetric_x_fixed_t', '2D_symmetric_y_fixed_y0', '3D_symmetric_x_fixed_x0', '3D_symmetric_xz_fixed_x0', '3D_symmetric_xz_fixed_z0', 'halo_orbit_fixed_z0', 'halo_orbit_fixed_x0']
__init__(dynamic, target=None, free_vars=None)[源代码]

初始化修正器

参数:
  • dynamic (CR3BP_Dynamics) -- CR3BP_Dynamics对象

  • target (dict[str, Any] | None) -- 目标约束条件字典(可选)

  • free_vars (list[str] | None) -- 自由变量列表(可选)

返回类型:

None

convergence_history: list[float]
error_history: list[float]
correction_history: list[float]
constraint_indices: list[int]
constraint_weights: dict[str, float]
constraint_types: dict[str, str]
free_variable_indices: list[int]
setup_type: str | None
symmetry_condition: str | None
fixed_parameters: dict[str, float]
setup_2D_symmetric_x_fixed_x0(x0=0.0)[源代码]

配置平面问题中固定初始x坐标的对称周期轨道搜索

在平面圆形限制性三体问题(PCRTBP)模型中,动力学方程关于会合坐标系的x轴具有对称性。 利用这一性质,周期轨道的搜索可以简化为寻找合适的初始条件: 从x轴上一点垂直出发(y=0, x_dot=0),经过半周期T/2后再次垂直穿越x轴(y=0, x_dot=0)。

本函数针对这种对称性设置微分修正问题,固定初始x坐标x0,将初始y方向速度y_dot0 和半周期T/2作为自由变量进行调整,以满足终点处的垂直穿越条件。

参数:

x0 -- 固定的初始x坐标,轨道从点(x0, 0)垂直出发, 这里将值设置为0.0是因为在这个函数中只需要使用x0进行初始化

返回:

配置好的微分修正器实例

备注

  • 自由变量: [y_dot0, T_half] - 初始y方向速度和半周期时间

  • 目标约束: [y(T/2)=0, x_dot(T/2)=0] - 终点处再次垂直穿越x轴

  • 状态向量索引: [1, 3] 分别对应y坐标和x方向速度

此配置对应于Broucke(1968)等经典文献中寻找对称周期轨道的基本方法, 可用于生成围绕平动点或主天体的各类周期轨道家族。

Reference:

Broucke R A. Periodic orbits in the restricted three body problem with Earth-moon masses[R]. 1968.

setup_2D_symmetric_x_fixed_t(t_half)[源代码]

配置平面问题中固定半周期的对称周期轨道搜索

固定半周期T/2,调整初始条件x0和y_dot0满足约束。

参数:

t_half -- 固定的半周期

返回:

配置好的微分修正器实例

setup_2D_symmetric_y_fixed_y0(y0=0.0)[源代码]

配置平面问题中固定初始y坐标的y轴对称周期轨道搜索

适用于共振轨道(RO)等从y轴出发的周期轨道。 轨道从点(0, y0)出发(x=0, x_dot=0),经过半周期T/2后再次穿越y轴(x=0, x_dot=0)。

参数:

y0 -- 固定的初始y坐标

返回:

配置好的微分修正器实例

备注

  • 自由变量: [x_dot0, T_half] - 初始x方向速度和半周期时间

  • 目标约束: [x(T/2)=0, x_dot(T/2)=0] - 终点处再次穿越y轴

  • 状态向量索引: [0, 3] 分别对应x坐标和x方向速度

setup_3D_symmetric_x_fixed_x0(x0)[源代码]

配置空间问题中固定初始x坐标的对称周期轨道搜索(如Halo轨道)

参数:

x0 (float) -- 固定的初始x坐标

返回:

配置好的微分修正器实例

返回类型:

self

setup_3D_symmetric_xz_fixed_x0(x0)[源代码]

配置空间XZ对称周期轨道搜索,固定X0

参数:

x0 (float) -- 固定的初始x坐标

返回:

配置好的微分修正器实例

返回类型:

self

setup_3D_symmetric_xz_fixed_z0(z0)[源代码]

配置空间XZ对称周期轨道搜索,固定Z0

参数:

z0 (float) -- 固定的初始z坐标

返回:

配置好的微分修正器实例

返回类型:

self

setup_halo_orbit_fixed_z0(z0, libration_point=1)[源代码]

配置 Halo 轨道微分修正,固定初始 Z0(XZ 对称)

Halo 轨道具有 XZ 平面对称性,利用该对称性可以将问题简化为: 从 XZ 平面上一点 (x0, 0, z0) 出发,经过半周期 T/2 后再次到达 XZ 平面。

参数:
  • z0 (float) -- 固定的初始 z 坐标

  • libration_point (int) -- 平动点编号 (1=L1, 2=L2),默认 L1

返回:

配置好的微分修正器实例

返回类型:

self

备注

  • 自由变量: [x0, y_dot0, T_half] - 初始 x 坐标、y 方向速度和半周期时间

  • 目标约束: [y(T/2)=0, x_dot(T/2)=0, z_dot(T/2)=0] - 半周期处再次位于 XZ 平面

  • 状态向量索引: [0, 4, 6] 分别对应 x0、y_dot0 和时间 T_half

  • 注意: z 在半周期时会改变符号,不作为约束

setup_halo_orbit_fixed_x0(x0, libration_point=1)[源代码]

配置 Halo 轨道微分修正,固定初始 X0(XZ 对称)

Halo 轨道具有 XZ 平面对称性,利用该对称性可以将问题简化为: 从 XZ 平面上一点 (x0, 0, z0) 出发,经过半周期 T/2 后再次到达 XZ 平面。

参数:
  • x0 (float) -- 固定的初始 x 坐标

  • libration_point (int) -- 平动点编号 (1=L1, 2=L2),默认 L1

返回:

配置好的微分修正器实例

返回类型:

self

备注

  • 自由变量: [z0, y_dot0, T_half] - 初始 z 坐标、y 方向速度和半周期时间

  • 目标约束: [y(T/2)=0, x_dot(T/2)=0, z_dot(T/2)=0] - 半周期处再次位于 XZ 平面且垂直穿越

  • 状态向量索引: [2, 4, 6] 分别对应 z0、y_dot0 和时间 T_half

iterate_correction(initial_guess, verbose=False)[源代码]

迭代修正主算法(基于STM的牛顿法)

通过状态转移矩阵(STM)构建雅可比矩阵,使用牛顿迭代法修正自由变量, 使终点状态满足目标约束条件,从而找到精确的周期轨道。

参数:
  • initial_guess (Orbit) -- 初始猜测轨道,或初始状态向量

  • verbose (bool) -- 是否打印迭代过程信息

返回:

  • 返回修正后的 Orbit 对象

check_convergence()[源代码]

检查收敛性

返回:

是否收敛

返回类型:

bool

get_convergence_history()[源代码]

获取收敛历史

返回:

收敛历史数据

返回类型:

dict

e2m2e.algorithms.multiple_shooting module

class e2m2e.algorithms.multiple_shooting.MultipleShootingResult(t_patch, state_patch, converged, iterations, max_residual, residual_history)[源代码]

基类:object

多重打靶法迭代修正的结果。

t_patch

修正后的时间节点数组,形状 (N,)

state_patch

修正后的状态量数组,形状 (N, 6),每行依次为 [x, y, z, vx, vy, vz]

converged

是否在最大迭代次数内收敛

iterations

实际迭代次数

max_residual

最终迭代的最大残差

residual_history

每次迭代最大残差的历史记录

class e2m2e.algorithms.multiple_shooting.MultipleShooting(dynamics, n_workers=1, kernel_dir=None)[源代码]

基类:object

多重打靶法(Multiple Shooting)修正器。

将一条轨迹分为 N 个节点、n_seg = N-1 段弧段,对每段独立积分后, 通过匹配相邻段端点状态来构建残差向量,再利用雅可比矩阵(含 STM) 进行最小二乘修正,反复迭代直到残差满足容差。

当 var_time=True 时,时间节点也作为自由变量参与修正(适用于自由时间问题)。

并行策略

  • n_workers=1:串行(默认)。

  • n_workers>1, kernel_dir=None:多线程(ThreadPoolExecutor)。 适合 CR3BP 等纯 Python/NumPy 动力学,但受 GIL 限制,并发收益有限。

  • n_workers>1, kernel_dir=<路径>**多进程**(ProcessPoolExecutor)。 每个子进程重载 SPICE 内核,绕过 GIL,可充分利用多核 CPU。 仅适用于 ``EphemerisDynamics``(需 SPICE 内核)。

__init__(dynamics, n_workers=1, kernel_dir=None)[源代码]
参数:
  • dynamics --

    动力学模型对象,需提供以下接口: - propagate(state, time_span, with_stm=True):

    积分传播,返回含 "states" 和 "stm" 的字典

    • equations_of_motion(t, state): 计算状态导数(右端函数值)

  • n_workers (int) --

    并行工作进程/线程数,默认 1(串行)。 - n_workers>1kernel_dir 已设置:

    使用 ProcessPoolExecutor(多进程,推荐)。

    • n_workers>1kernel_dir=None: 使用 ThreadPoolExecutor(多线程,受 GIL 限制)。

  • kernel_dir (str | None) -- SPICE 内核目录路径 (含 de440.bspnaif0012.tls)。 仅在 n_workers>1 时需要。设置后自动启用多进程模式。 例:kernel_dir=os.environ.get("SPICE_KERNEL_DIR", "../e2m2e/kernels")

返回类型:

None

correct(t_patch, state_patch, var_time=False, max_iter=None, tolerance=None, verbose=False)[源代码]

执行多重打靶修正。

将整条轨迹分为若干弧段,对每段独立积分后检验节点处的状态连续性, 利用状态转移矩阵(STM)组装雅可比矩阵,通过最小二乘求解修正量并迭代。

参数:
  • t_patch (ndarray) -- 初始时间节点数组,长度 N

  • state_patch (ndarray) -- 初始状态量数组,形状 (N, 6),每行 [x, y, z, vx, vy, vz]

  • var_time (bool) -- 是否允许时间节点作为自由变量参与修正

  • max_iter (int | None) -- 最大迭代次数(默认使用 self.max_iter)

  • tolerance (float | None) -- 收敛容差(默认使用 self.tolerance)

  • verbose (bool) -- 是否显示进度条

返回:

包含修正后的时间/状态、收敛标志、迭代次数和残差历史

返回类型:

MultipleShootingResult

参数:
  • n_workers (int)

  • kernel_dir (str | None)

e2m2e.algorithms.multiple_shooting.sample_patch_points(orbit, n_points)[源代码]

沿周期轨道均匀采样 patch points(打靶节点)。

该方法用于多重打靶法(Multiple Shooting)的前处理,将一条周期轨道 在时间上均匀分割为 n_points 个节点,并通过线性插值获取每个节点处的状态。

参数:
  • orbit -- 轨道对象,需包含以下属性: - period: 轨道周期(归一化时间单位) - times: 时间数组,形状 (M,) - states: 状态数组,形状 (M, 6),每行 [x, y, z, vx, vy, vz]

  • n_points (int) -- 需要采样的节点数量

返回:

包含两个数组的元组:
  • t_patch: 采样时间节点数组,形状 (n_points,),归一化时间单位

  • states: 采样状态数组,形状 (n_points, 6),每行 [x, y, z, vx, vy, vz]

返回类型:

Tuple[np.ndarray, np.ndarray]

抛出:

ValueError -- 当轨道对象没有 period 属性时抛出

备注

  • 采样时间范围为 [0, period),不包含周期终点(endpoint=False)

  • 使用线性插值从原始轨道数据中获取节点状态

  • 适用于 CR3BP 归一化坐标系下的周期轨道采样

e2m2e.algorithms.multiple_shooting.convert_to_j2000(t_patch_syn, states_syn, syn_j2000, reference_et, tu_days=4.34811305)[源代码]

将 synodic 坐标系下的 patch points 转换到 J2000 惯性坐标系。

该方法用于将 CR3BP 归一化 synodic 坐标系中的轨道节点转换到 J2000 惯性坐标系,以便在星历模型(Ephemeris)中进行高精度轨道修正。

参数:
  • t_patch_syn (ArrayLike) -- synodic 坐标系下的时间节点数组,归一化时间单位(TU)

  • states_syn (ArrayLike) -- synodic 坐标系下的状态数组,形状 (N, 6), 每行 [x, y, z, vx, vy, vz],归一化单位(DU, DU/TU)

  • syn_j2000 -- SynodicJ2000Transformation 对象,提供坐标转换功能

  • reference_et (float) -- 参考历元的 SPICE ephemeris time(ET),单位秒

  • tu_days (float) -- 归一化时间单位(TU)对应的天数,默认值为 4.34811305 天

返回:

包含两个数组的元组:
  • t_patch_j2000: J2000 坐标系下的时间数组,SPICE ET(秒)

  • states_j2000: J2000 坐标系下的状态数组,形状 (N, 6),

    每行 [x, y, z, vx, vy, vz],单位(km, km/s)

返回类型:

Tuple[np.ndarray, np.ndarray]

备注

  • 时间转换公式:t_j2000 = reference_et + t_syn * (tu_days * 86400)

  • 状态转换使用 SynodicJ2000Transformation.batch_synodic_to_j2000() 方法

  • 适用于将 CR3BP 轨道转换到星历模型进行高精度修正的场景

  • 转换后的状态可用于 EphemerisDynamics 进行轨道传播

e2m2e.algorithms.stability module

稳定性分析模块

提供轨道稳定性分析功能,包括单值矩阵计算、Floquet乘子分析、分岔检测等。

class e2m2e.algorithms.stability.StabilityType(*values)[源代码]

基类:Enum

稳定性类型枚举

STABLE = 'stable'
UNSTABLE = 'unstable'
MARGINALLY_STABLE = 'marginally_stable'
HYPERBOLIC = 'hyperbolic'
ELLIPTIC = 'elliptic'
PARABOLIC = 'parabolic'
class e2m2e.algorithms.stability.BifurcationType(*values)[源代码]

基类:Enum

分岔类型枚举

NONE = 'none'
PERIOD_DOUBLING = 'period_doubling'
SADDLE_NODE = 'saddle_node'
TORUS = 'torus'
PITCHFORK = 'pitchfork'
TRANSCRITICAL = 'transcritical'
SECONDARY_HOPF = 'secondary_hopf'
class e2m2e.algorithms.stability.StabilityAnalysis(orbit, dynamics=None)[源代码]

基类:object

轨道稳定性分析

计算轨道的单值矩阵、Floquet乘子、稳定性指数等, 并进行稳定性分类和分岔检测。

参数:
orbit

Orbit对象

dynamic

CR3BP_Dynamics对象

monodromy_matrix

单值矩阵

eigenvalues

特征值

stability_indices

稳定性指数

STABILITY_THRESHOLD = 1e-06
BIFURCATION_TOLERANCE = 1e-08
__init__(orbit, dynamics=None)[源代码]

初始化分析器

参数:
  • orbit (Orbit) -- Orbit对象

  • dynamics (CR3BP_Dynamics | None) -- CR3BP_Dynamics对象(可选,如果orbit关联了system则自动创建)

返回类型:

None

stm_history: list[ndarray]
eigenvalue_pairs: list[tuple[complex, complex]]
lyapunov_exponents: list[float] | ndarray
compute_monodromy()[源代码]

计算单值矩阵

通过积分一个完整周期的状态转移矩阵获得单值矩阵。

返回:

6x6 单值矩阵

返回类型:

np.ndarray

compute_floquet_multipliers()[源代码]

计算Floquet乘子(特征值)

返回:

Floquet乘子

返回类型:

np.ndarray

compute_stability_index()[源代码]

计算稳定性指数

Broucke稳定性参数定义为:ν = λ + 1/λ, 其中λ是单值矩阵的特征值。

返回:

稳定性指数字典

返回类型:

dict

classify_orbit()[源代码]

对轨道进行稳定性分类

返回:

稳定性分类结果

返回类型:

dict

analyze_bifurcation()[源代码]

分析分岔类型

通过检查单值矩阵特征值的分布判断是否存在分岔。

返回:

分岔分析结果

返回类型:

dict

full_analysis()[源代码]

执行完整的稳定性分析

返回:

完整分析结果

返回类型:

dict

static detect_bifurcation_in_family(orbits, dynamics, tolerance=1e-08)[源代码]

检测轨道族中的分岔点

遍历轨道族中的每条轨道,计算其单值矩阵特征值, 检测是否有特征值接近 +1(切分岔/saddle-node bifurcation)。

参数:
  • orbits (list[Orbit]) -- Orbit对象列表(轨道族)

  • dynamics (CR3BP_Dynamics) -- CR3BP_Dynamics对象

  • tolerance (float) -- 特征值接近 +1 的容差,默认 1e-8

返回:

分岔点列表,每个元素包含:
  • orbit_index: 轨道在族中的索引

  • orbit: Orbit对象

  • eigenvalues: 特征值数组

  • eigenvalue_diff: |λ - 1| 的最小值

  • bifurcation_type: 分岔类型

返回类型:

List[Dict[str, Any]]

static find_nearest_bifurcation(orbits, dynamics, target_x0=None, tolerance=0.0001)[源代码]

在轨道族中找到最接近目标参数的分岔点

参数:
  • orbits (list[Orbit]) -- Orbit对象列表

  • dynamics (CR3BP_Dynamics) -- CR3BP_Dynamics对象

  • target_x0 (float | None) -- 目标x0坐标(可选)

  • tolerance (float) -- 搜索容差

返回:

分岔点字典,如果未找到则返回None

返回类型:

dict[str, Any] | None

e2m2e.algorithms.strategies package

Differential correction strategy functions.

Each strategy function returns an immutable CorrectionConfig that fully describes the correction setup (symmetry, free variables, constraints, etc.). The DifferentialCorrection class delegates to these functions so that configuration logic is separated from the iterative solver.

class e2m2e.algorithms.strategies.CorrectionConfig(setup_type, symmetry_condition, fixed_parameters=<factory>, free_variables=<factory>, free_variable_indices=<factory>, target_conditions=<factory>, constraint_indices=<factory>, constraint_weights=<factory>, constraint_types=<factory>)[源代码]

基类:object

Immutable configuration for a differential correction strategy.

Encapsulates all correction parameters that were previously scattered across individual setup_* method bodies in DifferentialCorrection.

参数:
setup_type

Identifier string for the correction setup type.

Type:

str

symmetry_condition

Symmetry exploited by the correction (e.g. 'x_axis').

Type:

str

fixed_parameters

Parameter values held constant during correction.

Type:

dict[str, float]

free_variables

Names of variables the Newton solver adjusts.

Type:

list[str]

free_variable_indices

State-vector indices corresponding to free variables.

Type:

list[int]

target_conditions

Constraint names mapped to their target values.

Type:

dict[str, float]

constraint_indices

State-vector indices for constraint evaluation.

Type:

list[int]

constraint_weights

Per-constraint weighting factors for the Jacobian.

Type:

dict[str, float]

constraint_types

Per-constraint classification (e.g. 'equality').

Type:

dict[str, str]

setup_type: str
symmetry_condition: str
fixed_parameters: dict[str, float]
free_variables: list[str]
free_variable_indices: list[int]
target_conditions: dict[str, float]
constraint_indices: list[int]
constraint_weights: dict[str, float]
constraint_types: dict[str, str]
e2m2e.algorithms.strategies.symmetric_2d_fixed_x0(x0=0.0)[源代码]

Fixed x0: free variables are y_dot0 and T_half.

Used for planar symmetric periodic orbits that cross the x-axis perpendicularly at both start and half-period.

参数:

x0 (float) -- Fixed initial x coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_2d_fixed_t(t_half)[源代码]

Fixed half-period: free variables are x0 and y_dot0.

参数:

t_half (float) -- Fixed half-period value.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_2d_fixed_y0(y0=0.0)[源代码]

Fixed y0 (y-axis symmetric): free variables are x_dot0 and T_half.

Suitable for resonant orbits that depart from the y-axis.

参数:

y0 (float) -- Fixed initial y coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_3d_fixed_x0(x0)[源代码]

Fixed x0 in 3D: free variables are z0, y_dot0, T_half.

Used for spatially symmetric periodic orbits such as Halo orbits with the x-axis symmetry condition.

参数:

x0 (float) -- Fixed initial x coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_xz_fixed_x0(x0)[源代码]

Fixed x0 with XZ-plane symmetry: free variables are z0, y_dot0, T_half.

参数:

x0 (float) -- Fixed initial x coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_xz_fixed_z0(z0)[源代码]

Fixed z0 with XZ-plane symmetry: free variables are x0, y_dot0, T_half.

参数:

z0 (float) -- Fixed initial z coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.halo_fixed_z0(z0, libration_point=1)[源代码]

Halo orbit correction with fixed z0 (XZ-plane symmetry).

Free variables are x0, y_dot0, and T_half. The libration point is recorded as a fixed parameter for bookkeeping.

参数:
  • z0 (float) -- Fixed initial z coordinate.

  • libration_point (int) -- Lagrange point number (1=L1, 2=L2).

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.halo_fixed_x0(x0, libration_point=1)[源代码]

Halo orbit correction with fixed x0 (XZ-plane symmetry).

Free variables are z0, y_dot0, and T_half.

参数:
  • x0 (float) -- Fixed initial x coordinate.

  • libration_point (int) -- Lagrange point number (1=L1, 2=L2).

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.base module

Base configuration dataclass for differential correction strategies.

class e2m2e.algorithms.strategies.base.CorrectionConfig(setup_type, symmetry_condition, fixed_parameters=<factory>, free_variables=<factory>, free_variable_indices=<factory>, target_conditions=<factory>, constraint_indices=<factory>, constraint_weights=<factory>, constraint_types=<factory>)[源代码]

基类:object

Immutable configuration for a differential correction strategy.

Encapsulates all correction parameters that were previously scattered across individual setup_* method bodies in DifferentialCorrection.

参数:
setup_type

Identifier string for the correction setup type.

Type:

str

symmetry_condition

Symmetry exploited by the correction (e.g. 'x_axis').

Type:

str

fixed_parameters

Parameter values held constant during correction.

Type:

dict[str, float]

free_variables

Names of variables the Newton solver adjusts.

Type:

list[str]

free_variable_indices

State-vector indices corresponding to free variables.

Type:

list[int]

target_conditions

Constraint names mapped to their target values.

Type:

dict[str, float]

constraint_indices

State-vector indices for constraint evaluation.

Type:

list[int]

constraint_weights

Per-constraint weighting factors for the Jacobian.

Type:

dict[str, float]

constraint_types

Per-constraint classification (e.g. 'equality').

Type:

dict[str, str]

setup_type: str
symmetry_condition: str
fixed_parameters: dict[str, float]
free_variables: list[str]
free_variable_indices: list[int]
target_conditions: dict[str, float]
constraint_indices: list[int]
constraint_weights: dict[str, float]
constraint_types: dict[str, str]

e2m2e.algorithms.strategies.halo module

Halo orbit correction strategies.

e2m2e.algorithms.strategies.halo.halo_fixed_z0(z0, libration_point=1)[源代码]

Halo orbit correction with fixed z0 (XZ-plane symmetry).

Free variables are x0, y_dot0, and T_half. The libration point is recorded as a fixed parameter for bookkeeping.

参数:
  • z0 (float) -- Fixed initial z coordinate.

  • libration_point (int) -- Lagrange point number (1=L1, 2=L2).

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.halo.halo_fixed_x0(x0, libration_point=1)[源代码]

Halo orbit correction with fixed x0 (XZ-plane symmetry).

Free variables are z0, y_dot0, and T_half.

参数:
  • x0 (float) -- Fixed initial x coordinate.

  • libration_point (int) -- Lagrange point number (1=L1, 2=L2).

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_2d module

2D symmetric correction strategies (planar CR3BP).

e2m2e.algorithms.strategies.symmetric_2d.symmetric_2d_fixed_x0(x0=0.0)[源代码]

Fixed x0: free variables are y_dot0 and T_half.

Used for planar symmetric periodic orbits that cross the x-axis perpendicularly at both start and half-period.

参数:

x0 (float) -- Fixed initial x coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_2d.symmetric_2d_fixed_t(t_half)[源代码]

Fixed half-period: free variables are x0 and y_dot0.

参数:

t_half (float) -- Fixed half-period value.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_2d.symmetric_2d_fixed_y0(y0=0.0)[源代码]

Fixed y0 (y-axis symmetric): free variables are x_dot0 and T_half.

Suitable for resonant orbits that depart from the y-axis.

参数:

y0 (float) -- Fixed initial y coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_3d module

3D symmetric correction strategies (spatial CR3BP).

e2m2e.algorithms.strategies.symmetric_3d.symmetric_3d_fixed_x0(x0)[源代码]

Fixed x0 in 3D: free variables are z0, y_dot0, T_half.

Used for spatially symmetric periodic orbits such as Halo orbits with the x-axis symmetry condition.

参数:

x0 (float) -- Fixed initial x coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_3d.symmetric_xz_fixed_x0(x0)[源代码]

Fixed x0 with XZ-plane symmetry: free variables are z0, y_dot0, T_half.

参数:

x0 (float) -- Fixed initial x coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig

e2m2e.algorithms.strategies.symmetric_3d.symmetric_xz_fixed_z0(z0)[源代码]

Fixed z0 with XZ-plane symmetry: free variables are x0, y_dot0, T_half.

参数:

z0 (float) -- Fixed initial z coordinate.

返回:

CorrectionConfig with the corresponding correction parameters.

返回类型:

CorrectionConfig