轨道保持示例（main_control.py）#

设计 Halo 标称轨道后做轨道保持蒙特卡洛仿真，绘制标称与受控轨迹对比（会合系 x-z 投影，含五个平动点标注）。样本量刻意压小（2 个控制周期、 1 个蒙特卡洛样本），是快速演示版；工程评估把样本提到惯例值 100。

前置：SPICE 内核在仓库根 kernels/。

uv run --no-sync python examples/main_control.py --save

#!/usr/bin/env python3
"""main_control —— 轨道保持示例

先用 ``design_orbit`` 设计一条 Halo 标称轨道，再用 ``control_orbit``
做轨道保持蒙特卡洛仿真，绘制受控轨迹与标称轨迹对比。

用法：
    python examples/main_control.py            # 交互式出图
    python examples/main_control.py --save     # 存成 PNG（无头服务器可用）

前置条件：SPICE 内核位于仓库根 ``kernels/``（或设 ``$SPICE_KERNEL_DIR``）。
"""

from __future__ import annotations

import argparse
import pathlib

import numpy as np

# 输出图片保存到脚本所在目录（无论从哪运行）
_OUT_DIR = pathlib.Path(__file__).resolve().parent


def main() -> None:
    parser = argparse.ArgumentParser(description="轨道保持示例（Halo 蒙特卡洛）")
    parser.add_argument("--save", action="store_true", help="存为 PNG 而非交互式显示")
    parser.add_argument(
        "--log-level", default="WARNING", help="日志级别（DEBUG/INFO/WARNING/ERROR）"
    )
    args = parser.parse_args()

    from e2m2e.tools.logging import configure_logging

    configure_logging(level=args.log_level)

    print("=" * 60)
    print("e2m2e 轨道保持示例（Halo）")
    print("=" * 60)

    if args.save:
        import matplotlib

        matplotlib.use("Agg")

    from _plot_setup import setup_cjk_font

    setup_cjk_font()

    from e2m2e.algorithm.design import design_orbit
    from e2m2e.algorithm.station_keeping import control_orbit
    from e2m2e.api.models import DesignOrbitRequest

    # 1. 设计一条短弧 Halo 标称轨道（供轨道保持）
    print("\n1. 设计 L2 Halo 标称轨道")
    # 摄动开关：太阳第三体引力 + 地月非球形（10 阶）+ 炮弹模型光压
    perturbation = {
        "sun_body": 1,
        "planets": 0,
        "earth_nonspherical": 1,
        "moon_nonspherical": 1,
        "solar_radiation": 1,
        "atmosphere": 0,
        "relativity": 0,
        "tide": 0,
        "coupling": 0,
    }
    print("   摄动开关：")
    for k, v in perturbation.items():
        print(f"     {k} = {v}")
    result = design_orbit(
        DesignOrbitRequest(
            orbit_type="HALO",
            collinear_point=2,
            amplitude=30000.0,
            phase=0.0,
            duration=0.1 * 365.25 * 86400,
            output_step=3600.0,
            perturbation=perturbation,
        )
    )
    print(f"   星历行数 = {len(result.ephemeris)}")

    # 2. 轨道保持蒙特卡洛仿真（少量控制/样本，快速演示）
    print("\n2. 轨道保持蒙特卡洛仿真")
    print("   控制模式 1（目标点宽松），2 个控制周期，1 个蒙特卡洛样本")
    ctl = control_orbit(
        result.ephemeris,
        control_mode=1,
        num_controls=2,
        num_monte_carlo=1,
        control_interval=10.0,
        output_step=3600.0,
        perturbation=perturbation,
    )
    print(f"   失败样本数 = {ctl.num_failed}")
    rows = np.asarray(ctl.sk_statistic.rows)
    if rows.size:
        print(f"   总 Δv = {rows[0, 0]:.3f} m/s，最大 Δv = {rows[0, 1]:.3f} m/s")
    else:
        print("   统计为空（无有效样本）")

    # 3. 绘制标称 vs 受控轨迹（会合系 x-z）
    print("\n3. 绘制标称与受控轨迹对比")
    from e2m2e.algorithm.family.cr3bp_orbits import earth_moon_system

    system = earth_moon_system()

    # 标称：会合系无量纲状态
    nominal = result.ephemeris.synodic_position

    # 直接用 matplotlib 绘 x-z 投影：标称 vs 受控 + 地月天体 + 平动点
    import matplotlib.pyplot as plt

    from e2m2e.algorithm.dynamics import LibrationPoint

    fig, ax1 = plt.subplots(figsize=(12, 10), dpi=100)
    ax1.plot(
        nominal[:, 0],
        nominal[:, 2],
        label="标称轨道",
        linewidth=1.5,
        alpha=0.8,
    )

    # 受控：最后一次样本的受控星历（若可用）
    if ctl.controlled_ephemeris is not None:
        controlled = ctl.controlled_ephemeris.synodic_position
        ax1.plot(
            controlled[:, 0],
            controlled[:, 2],
            color="orange",
            label="受控轨道",
            linewidth=1.5,
            alpha=0.8,
        )

    # 天体标记（质心归一坐标：主天体在 -mu，次天体在 1-mu）
    mu = system.mu
    ax1.scatter(
        -mu,
        0,
        color="#2E86AB",
        s=200,
        edgecolors="#1A5276",
        linewidth=1.5,
        zorder=10,
        label="Earth",
    )
    ax1.scatter(
        1 - mu,
        0,
        color="#95A5A6",
        s=100,
        edgecolors="#566573",
        linewidth=1.5,
        zorder=10,
        label="Moon",
    )

    # 五个平动点（灰色三角 + 标签）
    for i, lp in enumerate(LibrationPoint):
        coord = system.L_points[lp]
        ax1.scatter(coord[0], coord[1], color="gray", marker="^", s=60, zorder=5)
        ax1.annotate(
            f"L{i + 1}",
            (coord[0], coord[1]),
            textcoords="offset points",
            xytext=(5, 5),
            fontsize=16,
        )

    ax1.set_xlabel("X（无量纲）")
    ax1.set_ylabel("Z（无量纲）")
    ax1.set_title("L2 Halo 轨道保持：标称 vs 受控（会合系 x-z）")
    ax1.legend(loc="upper right")

    if args.save:
        fig.savefig(
            str(_OUT_DIR / "main_control_halo.png"),
            dpi=150,
            bbox_inches="tight",
            pad_inches=0.1,
        )
        print(f"   已保存 {_OUT_DIR / 'main_control_halo.png'}")
    else:
        plt.show()

    print("\n" + "=" * 60)
    print("示例完成！")
    print("=" * 60)


if __name__ == "__main__":
    main()

产出与解读#

图：examples/main_control_halo.png——标称轨道与受控轨道的 x-z 投影对比，地月天体与 L1–L5 标注同坐标系。
结果字段：ctl.num_failed 是蒙特卡洛失败样本数； ctl.sk_statistic.rows 首行前两列是总 Δv 与最大单次 Δv（m/s）； ctl.controlled_ephemeris.synodic_position 是受控星历（全部样本失败时为 None），画图前判空。
control_orbit(result.ephemeris, control_mode=1, num_controls=2, num_monte_carlo=1, control_interval=10.0, ...)：控制模式 1 = 目标点宽松；控制力模型与真实力模型的摄动开关可分别指定。

对照：接口层等价调用见轨道保持。