基于Gekko的轨道推力与角度非线性优化不收敛问题求助

我开发了一个轨道推力与角度的非线性优化求解程序，但始终无法收敛。目标是通过优化各时间步的thrust与angle控制变量，最小化以推力时间积分定义的燃料消耗。尝试过Gekko库及scipy.optimize.minimize求解，均无法实现目标函数最小化，调整时间尺度也无效。

Python实现代码 #

from gekko import GEKKO
import numpy as np

# Initialize Model
m = GEKKO(remote=True)

m.time = np.linspace(0,140,801)

# State Variables
r = m.Var(value = 1.0)
theta = m.Var(value = 0.0)
vr = m.Var(value = 0.0)
vt = m.Var(value = 1.0)

# Control Variables
thrust = m.MV(value = 0.005, lb = 0.0, ub = 0.01)
thrust.STATUS = 1
angle = m.MV(value = 0.0, lb = -np.pi/2, ub = np.pi/2)
angle.STATUS = 1

# Optimise fuel consumption
fuel_consumption = m.Var(value = 0.0)

# Dynamics
m.Equation(r.dt() == vr)
m.Equation(theta.dt() == vt/r)
m.Equation(vr.dt() == (vt**2)/r - 1/(r**2) + thrust*m.sin(angle))
m.Equation(vt.dt() == -(vr*vt)/r + thrust*m.cos(angle))

m.Equation(fuel_consumption.dt() == thrust)  # Accumulate thrust over time

# Boundary Conditions
m.fix(r, pos=len(m.time)-1,val = 4.0)
m.fix(theta, pos=len(m.time)-1,val = 0.0)
m.fix(vr, pos=len(m.time)-1,val = 0)
m.fix(vt, pos=len(m.time)-1,val = 0.5)

# Objective function
m.Obj(fuel_consumption)

#Set global options
m.options.IMODE = 6  # Dynamic Optimization (Simultaneous ODEs)
m.options.NODES = 3  # Collocation nodes
m.options.SOLVER = 3  # IPOPT solver

#Solve simulation
m.solve(disp=True) # solve on public server

#Results
print('')
print('Results')
print('Objective: ',m.options.OBJFCNVAL)
print('Solution: ', thrust.value)

求解日志片段 #

----------------------------------------------------------------
 APMonitor, Version 1.0.1
 APMonitor Optimization Suite
 ----------------------------------------------------------------
 
 
 --------- APM Model Size ------------
 Each time step contains
   Objects      :            0
   Constants    :            0
   Variables    :            7
   Intermediates:            0
   Connections  :            8
   Equations    :            6
   Residuals    :            6
 
 Number of state variables:          25592
 Number of total equations: -        24000
 Number of slack variables: -            0
 ---------------------------------------
 Degrees of freedom       :           1592
 
 **********************************************
 Dynamic Control with Interior Point Solver
 **********************************************
  
  
 Info: Exact Hessian

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************

This is Ipopt version 3.12.10, running with linear solver ma57.

Number of nonzeros in equality constraint Jacobian...:    67966
Number of nonzeros in inequality constraint Jacobian.:     6400
Number of nonzeros in Lagrangian Hessian.............:    11195

Total number of variables............................:    25592
                     variables with only lower bounds:     3200
                variables with lower and upper bounds:     3200
                     variables with only upper bounds:        0
Total number of equality constraints.................:    20800
Total number of inequality constraints...............:     3200
        inequality constraints with only lower bounds:     3200
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.3103987e-02 3.00e+00 1.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
Reallocating memory for MA57: lfact (641403)
   1  3.8902526e+01 2.77e+00 4.18e+00  -1.3 1.47e+02    -  2.34e-01 7.51e-02f  1

......

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 250r 2.8607519e+02 1.73e+02 9.99e+02   2.2 0.00e+00   3.1 0.00e+00 2.58e-07R  5

Number of Iterations....: 250

                                   (scaled)                 (unscaled)
Objective...............:   2.8607519277584993e+02    2.8607519277584993e+02
Dual infeasibility......:   2.1703502621114161e+00    2.1703502621114161e+00
Constraint violation....:   1.7251865116637478e+02    1.7251865116637478e+02
Complementarity.........:   3.1072980675208779e+00    3.1072980675208779e+00
Overall NLP error.......:   1.7251865116637478e+02    1.7251865116637478e+02


Number of objective function evaluations             = 323
Number of objective gradient evaluations             = 252
Number of equality constraint evaluations            = 323
Number of inequality constraint evaluations          = 323
Number of equality constraint Jacobian evaluations   = 254
Number of inequality constraint Jacobian evaluations = 254
Number of Lagrangian Hessian evaluations             = 250
Total CPU secs in IPOPT (w/o function evaluations)   =     26.390
Total CPU secs in NLP function evaluations           =     45.994

EXIT: Maximum Number of Iterations Exceeded.

问题分析与解决建议 #

1. 建模潜在问题 #

• 终端约束过严：theta=0的终端约束不合理。轨道转移过程中极角必然会发生变化，强制终端回到初始值可能导致问题无解或收敛困难。可改为允许小范围偏差，或改用相对角度约束。
• 动力学方程数值不稳定：1/(r²)项在r较小时数值波动剧烈，建议对所有变量（r, vr, vt, thrust）进行无量纲化处理，缩小数值范围，提升求解稳定性。
• 初始值偏离可行域：thrust和angle的初始值可能远离最优解区域。先通过仿真模式（IMODE=4）验证从初始状态到终端状态的可行性，再以此为基础设置优化初始值。

2. 求解器参数调整 #

• 增加迭代次数：当前IPOPT迭代上限250，可通过m.options.MAX_ITER=500设置更高数值，给求解器更多收敛时间。
• 限制控制变量波动：给MV变量添加DCOST项（如thrust.DCOST=1e-6），惩罚控制变量的剧烈变化，避免求解震荡。
• 更换求解器：尝试使用SOLVER=1（APOPT），它在处理非线性问题时比IPOPT更稳健；或更换IPOPT的线性求解器为ma27。
• 减少时间节点：当前801个时间节点导致变量规模过大（25592个变量），可先减少到201个节点验证收敛性，再逐步提升精度。

3. 目标函数简化 #

• 直接使用m.Obj(m.integral(thrust))代替额外状态变量fuel_consumption，简化模型结构，降低求解复杂度。

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