基于GEKKO+Python调整日期变量最小化delta-V impulse

问题核心 #

这是Gekko+IPOPT轨道优化问题的扩展，目标是通过调整发射日期、飞掠日期、抵达日期（已定义为Gekko操作变量MV）最小化delta-V冲量。当前实现的slingshot函数仅返回初始猜测日期对应的结果，未触发优化迭代，需要调整代码让Gekko能通过该函数优化日期参数。

问题根源 #

直接调用slingshot()时，代码访问了launch.value、flyby.value、arrive.value，这只会获取MV变量的初始值，而非让Gekko将这些变量作为优化参数传入函数迭代计算。Gekko无法识别直接访问.value的操作，必须通过包装函数或建立约束的方式，让优化器能动态调整MV变量并计算对应冲量。

解决方案 #

1. 重构轨道计算函数 #

修改slingshot函数，让它接受日期数值作为输入参数，而非直接访问Gekko变量：

def slingshot(launch_jd, flyby_jd, arrive_jd):
    # 日期转换
    date_launchE = Time(launch_jd, format="jd", scale="utc").tdb
    date_flyby = Time(flyby_jd, format="jd", scale="utc").tdb
    date_arrivalE = Time(arrive_jd, format="jd", scale="utc").tdb

    # 天体轨道数据获取
    venus = Ephem.from_body(Venus, time_range(date_launchE, end=date_arrivalE, periods=500))
    earth = Ephem.from_body(Earth, time_range(date_launchE, end=date_arrivalE, periods=500))

    r_e0, v_e0 = earth.rv(date_launchE)

    # 初始轨道与飞掠轨道构建
    ss_e0 = Orbit.from_ephem(Sun, earth, date_launchE)
    ss_fly = Orbit.from_ephem(Sun, venus, date_flyby)

    # Lambert转移轨道计算（地球到飞掠天体）
    man_launch = Maneuver.lambert(ss_e0, ss_fly)
    cruise1 = ss_e0.apply_maneuver(man_launch)
    cruise1_end = cruise1.propagate(date_flyby)

    # 飞掠到地球的Lambert转移轨道
    ss_earrive = Orbit.from_ephem(Sun, earth, date_arrivalE)
    man_flyby = Maneuver.lambert(cruise1_end, ss_earrive)
    imp_a, imp_b = man_flyby.impulses

    # 轨道传播与最终状态获取
    cruise2, _ = cruise1_end.apply_maneuver(man_flyby, intermediate=True)
    cruise2_end = cruise2.propagate(date_arrivalE)
    ss_target = Orbit.from_ephem(Sun, earth, date_arrivalE)

    # 提取各节点的位置、速度与冲量
    r1, v1 = Ephem.from_orbit(ss_e0, date_launchE).rv(date_launchE)
    r2, v2 = Ephem.from_orbit(ss_fly, date_flyby).rv(date_flyby)
    r3, v3 = Ephem.from_orbit(ss_target, date_arrivalE).rv(date_arrivalE)
    
    # 转换为数值列表并计算delta-V模长
    r1_vals = [x.value for x in r1]
    v1_vals = [x.value for x in v1]
    r2_vals = [x.value for x in r2]
    v2_vals = [x.value for x in v2]
    r3_vals = [x.value for x in r3]
    v3_vals = [x.value for x in v3]
    imp_a_vals = [x.value for x in imp_a]
    delta_v = np.linalg.norm(imp_a_vals)
    
    return delta_v, r1_vals, v1_vals, r2_vals, v2_vals, r3_vals, v3_vals

2. 用Gekko自定义函数包装实现优化 #

由于轨道计算依赖外部库（无法自动求导），需用Gekko的user_defined功能将函数包装为优化器可调用的黑箱函数，建立目标与变量约束：

from gekko import GEKKO
import numpy as np

# 初始化Gekko模型
m = GEKKO(remote=False)

# 定义操作变量（MV）
launch = m.MV(value=2460159.5, lb=2460159.5, ub=2460525.5)
launch.STATUS = 1
flyby = m.MV(value=2460846.5, lb=2460704.5, ub=2460908.5)
flyby.STATUS = 1
arrive = m.MV(value=2461534.5, lb=2461250.5, ub=2461658.5)
arrive.STATUS = 1

# 定义需跟踪的结果变量
delta_v = m.Var(value=0)
r1 = m.Array(m.Var, 3, value=0)
v1 = m.Array(m.Var, 3, value=0)
r2 = m.Array(m.Var, 3, value=0)
v2 = m.Array(m.Var, 3, value=0)
r3 = m.Array(m.Var, 3, value=0)
v3 = m.Array(m.Var, 3, value=0)

# 定义Gekko兼容的包装函数
def slingshot_gekko():
    # 获取当前迭代的日期数值
    lj = launch.value[0]
    fj = flyby.value[0]
    aj = arrive.value[0]
    # 调用轨道计算函数
    dv, r1v, v1v, r2v, v2v, r3v, v3v = slingshot(lj, fj, aj)
    # 返回所有结果（顺序需与user_defined输出列表对应）
    return [dv] + r1v + v1v + r2v + v2v + r3v + v3v

# 注册自定义函数：输入为3个MV，输出为delta-V+所有位置速度变量
m.user_defined(slingshot_gekko, 
               inputs=[launch, flyby, arrive], 
               outputs=[delta_v] + r1 + v1 + r2 + v2 + r3 + v3)

# 设置优化目标：最小化delta-V
m.Minimize(delta_v)

# 设置静态优化模式
m.IMODE = 2

# 运行优化（disp=True显示迭代过程）
m.solve(disp=True)

# 输出优化结果
print("=== 优化结果 ===")
print(f"发射日期: {launch.value[0]:.2f} JD")
print(f"飞掠日期: {flyby.value[0]:.2f} JD")
print(f"抵达日期: {arrive.value[0]:.2f} JD")
print(f"最小delta-V: {delta_v.value[0]:.4f} km/s")
print(f"发射点位置r1: {[round(x.value[0], 2) for x in r1]}")
print(f"发射点速度v1: {[round(x.value[0], 4) for x in v1]}")
# 其他变量可按同样方式打印

3. 关键注意事项 #

• 禁止直接在优化逻辑中访问Gekko变量的.value，需通过user_defined让Gekko在每次迭代时自动传递当前MV数值。
• 由于是黑箱优化，Gekko会用有限差分法计算梯度，可调整m.options.RTOL、m.options.OTOL参数平衡精度与速度。
• 轨道计算的periods参数可适当减小，加快函数运行速度，提升优化迭代效率。

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