GEKKO+IPOPT航天器轨道传播优化飞行时长问题咨询

用户问题描述 #

尝试用GEKKO/Python的IPOPT求解器优化航天器轨道飞行时长，核心需求是从初始轨道oe_i以delta_t=3600s步长传播到目标轨道oe_f，但轨道传播部分存在逻辑错误，具体问题：

• 无循环时模型可运行但无有效求解结果（未实现长时间传播）
• 用for循环时出现tf非整数错误
• 基于tf计算值循环时，mass变量无法参与计算

原代码如下：

m = GEKKO()

#manipulating variables and initial guesses
al_a = m.MV(value = -1, lb = -2, ub = 2)
al_a.STATUS = 1
l_e = m.MV(value = 0.001, lb = 0, ub = 10**6)
l_e.STATUS = 1
l_i = m.MV(value = 1, lb = 0, ub = 10**6)
l_i.STATUS = 1

#variables and initial guesses
a = m.Var(value = oe_i[0], lb = oe_i[0] - 6378000, ub = oe_f[0] + 6378000, name = 'sma')
e = m.Var(value = oe_i[1], lb = 0, ub = 1, name = 'ecc')
i = m.Var(value = oe_i[2], lb = 0, ub = math.radians(90), name = 'inc')
Om = m.Var(value = oe_i[3], lb = 0, ub = math.radians(360), name = 'raan')
om = m.Var(value = oe_i[4], lb = 0, ub = math.radians(360), name = 'ap')
nu = m.Var(value = oe_i[5], lb = 0, ub = math.radians(360), name = 'ta')
mass = m.Var(value = m0, lb = 0, ub = m0, name = 'mass')

#objective function
tf = m.FV(1.2 * ((m0 - mass)/dm), lb = 0, ub = t_max)
tf.STATUS = 1

#propagate
t = 0

while t <= tf:
    deltas, Tp = propagate(a, e, i, Om, om, nu, mass)
    m.Equation(Tp * a.dt() == (deltas[0] * delta_t * deltas[7]))
    m.Equation(Tp * e.dt() == (deltas[1] * delta_t * deltas[7]))
    m.Equation(Tp * i.dt() == (deltas[2] * delta_t * deltas[7]))
    m.Equation(Tp * Om.dt() == (deltas[3] * delta_t * deltas[7]))
    m.Equation(Tp * om.dt() == (deltas[4] * delta_t * deltas[7]))
    m.Equation(nu.dt() == deltas[5] * delta_t)
    m.Equation(Tp * mass.dt() == (deltas[6] * delta_t * deltas[7]))
    
    t = t + delta_t

#starting constraints
m.fix(a, pos = 0, val = oe_i[0])
m.fix(e, pos = 0, val = oe_i[1])
m.fix(i, pos = 0, val = oe_i[2])
m.fix(Om, pos = 0, val = oe_i[3])
m.fix(om, pos = 0, val = oe_i[4])
m.fix(nu, pos = 0, val = oe_i[5])
m.fix(mass, pos = 0, val = m0)

#boundary constraints
m.fix(a, pos = len(m.time) - 1, val = oe_f[0])
m.fix(e, pos = len(m.time) - 1, val = oe_f[1])
m.fix(i, pos = len(m.time) - 1, val = oe_f[2])
m.fix(Om, pos = len(m.time) - 1, val = oe_f[3])
m.fix(om, pos = len(m.time) - 1, val = oe_f[4])
m.fix(nu, pos = len(m.time) - 1, val = oe_f[5])
m.fix(mass, pos = len(m.time) - 1, val = 0)

m.time = np.linspace(0,0.2,100)
m.Obj(tf)
m.options.IMODE = 6 # non-linear model
m.options.SOLVER = 3 # solver (IPOPT)
m.options.MAX_ITER = 15000
m.options.RTOL = 1e-7
m.options.OTOL = 1e-7
m.open_folder()
m.solve(display=false) # Solve
print('Optimal time: ' + str(tf.value[0]))
m.solve()
m.open_folder(infeasibilities.txt)

核心错误分析 #

• 循环逻辑误用GEKKO机制：GEKKO动态优化（IMODE=6）基于m.time数组自动离散化时间步，无需手动编写while/for循环。手动循环会导致微分方程无法与时间维度绑定，且tf是符号优化变量，不能作为Python循环的数值判断条件。
• tf定义错误：tf = m.FV(1.2 * ((m0 - mass)/dm))中，mass是随时间变化的状态变量，不能用于初始化FV的数值，FV初始值必须为常数。
• 时间尺度不匹配：m.time = np.linspace(0,0.2,100)的时间范围与实际需求的3600s步长完全脱节，无法覆盖轨道传播的实际时长。

修正方案 #

1. 重构时间与微分方程定义 #

去掉手动循环，利用GEKKO自动离散化机制，通过m.time生成固定步长的时间数组：

delta_t = 3600  # 步长3600s
num_steps = int(t_max / delta_t) + 1
# 先预设足够大的时间范围，tf会自动优化到最优值
m.time = np.linspace(0, t_max, num_steps)

2. 修正tf的定义 #

将tf设为独立FV变量，通过约束关联质量消耗：

# 初始值设为初始轨道周期的合理倍数
tf = m.FV(value=2*np.pi*np.sqrt(oe_i[0]**3/mu), lb=0, ub=t_max)
tf.STATUS = 1

# 质量消耗与tf的关联约束（dm为单位时间质量消耗率）
m.Equation(mass[-1] == m0 - dm * tf)

3. 正确定义轨道传播微分方程 #

直接将propagate返回的速率项与变量导数绑定，无需手动处理步长：

deltas, Tp = propagate(a, e, i, Om, om, nu, mass)
# 若propagate返回的是步长增量，需除以delta_t转换为速率
m.Equation(a.dt() == deltas[0] * deltas[7] / Tp)
m.Equation(e.dt() == deltas[1] * deltas[7] / Tp)
m.Equation(i.dt() == deltas[2] * deltas[7] / Tp)
m.Equation(Om.dt() == deltas[3] * deltas[7] / Tp)
m.Equation(om.dt() == deltas[4] * deltas[7] / Tp)
m.Equation(nu.dt() == deltas[5])
m.Equation(mass.dt() == deltas[6] * deltas[7] / Tp)

4. 修正边界约束 #

用[-1]直接引用终端时间点的变量值，避免手动计算索引：

# 初始约束
m.fix(a, pos=0, val=oe_i[0])
m.fix(e, pos=0, val=oe_i[1])
m.fix(i, pos=0, val=oe_i[2])
m.fix(Om, pos=0, val=oe_i[3])
m.fix(om, pos=0, val=oe_i[4])
m.fix(nu, pos=0, val=oe_i[5])
m.fix(mass, pos=0, val=m0)

# 终端约束（放宽质量约束避免不可行）
m.fix(a, pos=-1, val=oe_f[0])
m.fix(e, pos=-1, val=oe_f[1])
m.fix(i, pos=-1, val=oe_f[2])
m.fix(Om, pos=-1, val=oe_f[3])
m.fix(om, pos=-1, val=oe_f[4])
m.fix(nu, pos=-1, val=oe_f[5])
m.Equation(mass[-1] >= 0)

5. 修正求解调用 #

去掉重复求解，修正display参数大小写：

m.Obj(tf)
m.options.IMODE = 6
m.options.SOLVER = 3
m.options.MAX_ITER = 15000
m.options.RTOL = 1e-7
m.options.OTOL = 1e-7
m.solve(display=False)
print('Optimal time: ' + str(tf.value[0]))
m.open_folder()

同类参考示例 #

GEKKO最优控制的典型结构示例：

import numpy as np
from gekko import GEKKO

m = GEKKO()
# 定义时间数组
m.time = np.linspace(0, 10, 100)
# 定义变量
u = m.MV(value=0, lb=-1, ub=1)
u.STATUS = 1
x = m.Var(value=0)
tf = m.FV(value=10, lb=5, ub=15)
tf.STATUS = 1
# 时间缩放（适配tf优化）
m.time = m.time * tf
# 微分方程
m.Equation(x.dt() == u)
# 边界约束
m.fix(x, pos=0, val=0)
m.fix(x, pos=-1, val=1)
# 目标函数
m.Obj(tf)
# 求解
m.options.IMODE = 6
m.solve(display=False)

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