python - GEKKO Python + IPOPT 最小化 delta-V

Question

我的 Python 中的 GEKKO 模型如下:

# 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 = 2460887.5, ub = 2460908.5) # 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])) # Venus/Mars
# arrive = m.MV(value = 2461534.5, lb = 2461250.5, ub = 2461658.5) # 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' 是速度变化。

Answer 1

m.time 向量需要与用于提供value 的 numpy 数组具有相同的维度。 value 可以是单个值，也可以是 m.time 的大小。在这种情况下，需要定义3个时间点。要使用 img[0]...img[2]，需要定义一个数组来创建该变量列表。每个变量都有 3 个值，这些值在相应的 m.time 点处求解。 slingshot() 函数未定义。即使已定义，也不应使用它来重新定义 Gekko 变量。

下面是成功解决的修改，但是，这可能不是预期的问题。

import numpy as np
from gekko import GEKKO

# Initialize model
m = GEKKO()

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

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

# Objective function
dV = m.Minimize(m.sqrt(imp[0]**2 + imp[1]**2 + imp[2]**2))

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

m.time = [0,1,2]
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=True) # Solve

解决方案:

  36  5.8037359e-06 9.31e-10 4.02e-08  -6.6 2.38e+03 -12.1 1.00e+00 9.28e-03f  1
  37  4.8012697e-09 0.00e+00 5.99e-11 -11.0 1.33e-01 -12.6 9.96e-01 9.99e-01h  1

Number of Iterations....: 37

                                   (scaled)                 (unscaled)
Objective...............:   4.8012696909311766e-09    4.8012696909311766e-09
Dual infeasibility......:   5.9928342286739218e-11    5.9928342286739218e-11
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.1687338056472621e-09    1.1687338056472621e-09
Overall NLP error.......:   1.1687338056472621e-09    1.1687338056472621e-09


Number of objective function evaluations             = 38
Number of objective gradient evaluations             = 38
Number of equality constraint evaluations            = 38
Number of inequality constraint evaluations          = 38
Number of equality constraint Jacobian evaluations   = 38
Number of inequality constraint Jacobian evaluations = 38
Number of Lagrangian Hessian evaluations             = 37
Total CPU secs in IPOPT (w/o function evaluations)   =      0.021
Total CPU secs in NLP function evaluations           =      0.004

EXIT: Optimal Solution Found.
 
 The solution was found.
 
 The final value of the objective function is   4.801269690931177E-009
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   3.239999999641441E-002 sec
 Objective      :   4.801269690931177E-009
 Successful solution
 ---------------------------------------------------

如果这不是正确的解决方案，请使用此信息并可能创建一个新问题，并尽力解决该问题。有Gekko tutorials可以更完整地概述如何配置动态优化问题。