GEKKO+IPOPT最小化Δ-V时数组长度匹配错误求助

我用Python构建了含引力弹弓的轨道转移优化GEKKO模型，目标是最小化Δ-V，采用IPOPT求解器，代码如下：

# Initialize model
m = GEKKO()

# Manipulating variables and initial guesses
launch = m.MV(value = np.array([2460310.5, 0, 0]), lb = np.array([2460310.5, 0, 0]), ub = np.array([2460340.5, 0, 0]))
launch.STATUS = 1
flyby = m.MV(value = np.array([2460575.5, 0, 0]), lb = np.array([2460493.5, 0, 0]), ub = np.array([2460340.5, 0, 0])) # Venus/Mars
# flyby = m.MV(value = 2460997.5, lb = np.array([2460887.5, 0, 0]), ub = np.array([2460908.5, 0, 0])) # Jupiter
flyby.STATUS = 1
arrive = m.MV(value = np.array([2460845.5, 0, 0]), lb = np.array([2460631.5, 0, 0]), ub = np.array([2460660.5, 0, 0])) # Venus/Mars
# arrive = m.MV(value = 2461534.5, lb = np.array([2461250.5, 0, 0]), ub = np.array([2461658.5, 0, 0])) # Jupiter
arrive.STATUS = 1

# Variables
r1 = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e10, -1e10, -1e10]), ub = np.array([1e10, 1e10, 1e10]), name = "r1")
v1 = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e5, -1e5, -1e5]), ub = np.array([1e5, 1e5, 1e5]), name = "v1")
r2 = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e10, -1e10, -1e10]), ub = np.array([1e10, 1e10, 1e10]), name = "r2")
v2 = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e5, -1e5, -1e5]), ub = np.array([1e5, 1e5, 1e5]), name = "v2")
r3 = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e10, -1e10, -1e10]), ub = np.array([1e10, 1e10, 1e10]), name = "r3")
v3 = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e5, -1e5, -1e5]), ub = np.array([1e5, 1e5, 1e5]), name = "v3")
l = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e5, -1e5, -1e5]), ub = np.array([1e5, 1e5, 1e5]), name = "launch")
imp = m.Var(value = np.array([0, 0, 0]), lb = np.array([-1e5, -1e5, -1e5]), ub = np.array([1e5, 1e5, 1e5]), name = "impulse")

# Objective function
dV = m.FV(value = m.sqrt(imp[0]**2 + imp[1]**2 + imp[2]**2), lb = 0, ub = 10000)
dV.STATUS = 1

# Slingshot maneuver
r1, v1, r2, v2, r3, v3, l_mag, imp_mag, v_final = slingshot()

m.Obj(dV) #minimize delta-V
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.solve(disp=False) # Solve

运行时抛出错误：

Exception: Data arrays must have the same length, and match time discretization in dynamic problems

我尝试过修改m.Var变量的数据类型、将m.MV改为数组、移除/保留m.time等方法，但均未解决问题。其中r为位置矢量，v为速度矢量，l和imp为速度增量。

问题根源与修正方案 #

核心问题分析

报错的本质是GEKKO变量/参数的维度与模型模式不匹配，代码里存在几个关键错误：

1. MV变量错误使用数组初始化：launch、flyby、arrive是时间点（标量），但你用了3维数组赋值，在IMODE=6（稳态优化）下会引发维度冲突。
2. FV变量误用：dV是目标函数，不需要用m.FV（FV用于参数估计），直接定义目标即可。
3. slingshot()函数变量覆盖：你重新赋值了已定义的GEKKO变量，导致模型无法识别变量关联；如果函数返回普通数组而非GEKKO表达式，也会引发维度问题。
4. 矢量变量定义不规范：位置、速度这类3维矢量，直接用m.Var赋值数组会导致GEKKO无法正确识别分量维度。

具体修正步骤

1. 修正MV变量为标量：

   launch = m.MV(value=2460310.5, lb=2460310.5, ub=2460340.5)
   flyby = m.MV(value=2460575.5, lb=2460493.5, ub=2460340.5)
   arrive = m.MV(value=2460845.5, lb=2460631.5, ub=2460660.5)
   

2. 修正目标函数定义：

   # 直接将速度增量的模作为优化目标，无需FV
   m.Obj(m.sqrt(imp[0]**2 + imp[1]**2 + imp[2]**2))
   

3. 规范矢量变量定义：用m.Array创建3维矢量变量，方便GEKKO识别分量：

   r1 = m.Array(m.Var, 3, value=0, lb=-1e10, ub=1e10, name="r1")
   v1 = m.Array(m.Var, 3, value=0, lb=-1e5, ub=1e5, name="v1")
   r2 = m.Array(m.Var, 3, value=0, lb=-1e10, ub=1e10, name="r2")
   v2 = m.Array(m.Var, 3, value=0, lb=-1e5, ub=1e5, name="v2")
   r3 = m.Array(m.Var, 3, value=0, lb=-1e10, ub=1e10, name="r3")
   v3 = m.Array(m.Var, 3, value=0, lb=-1e5, ub=1e5, name="v3")
   l = m.Array(m.Var, 3, value=0, lb=-1e5, ub=1e5, name="launch")
   imp = m.Array(m.Var, 3, value=0, lb=-1e5, ub=1e5, name="impulse")
   

4. 修复slingshot()函数调用：让函数接收GEKKO变量作为输入，通过m.Equation()绑定约束关系，而非覆盖变量：

   def slingshot(m, launch, flyby, arrive, r1, v1, r2, v2, r3, v3, l, imp):
       # 在这里添加引力弹弓的约束方程，例如轨道动力学、位置速度匹配等
       # m.Equation(r1[0] == ...)
       # m.Equation(v1[1] == ...)
       # 返回GEKKO表达式而非普通数组
       l_mag = m.sqrt(l[0]**2 + l[1]**2 + l[2]**2)
       imp_mag = m.sqrt(imp[0]**2 + imp[1]**2 + imp[2]**2)
       v_final = v3 + imp
       return r1, v1, r2, v2, r3, v3, l_mag, imp_mag, v_final
   
   # 调用函数并绑定约束
   r1, v1, r2, v2, r3, v3, l_mag, imp_mag, v_final = slingshot(m, launch, flyby, arrive, r1, v1, r2, v2, r3, v3, l, imp)
   

5. 确认模型模式：如果是稳态优化（仅优化关键时间点的状态），保持IMODE=6；如果是动态轨迹优化，需要定义m.time并切换到动态模式（如IMODE=9），同时确保所有变量的初始值数组长度与m.time一致。

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