Python GEKKO-IPOPT轨道优化代码收敛故障求助

问题 #

使用GEKKO编写轨道优化代码，调用IPOPT求解时触发错误：

Restoration phase is called at point that is almost feasible, with constraint violation 0.000000e+00. Abort.

已尝试调整变量类型、MV初始值、目标函数形式及优化方向，但均未解决问题，需实现求解收敛。

完整代码 #

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math

#constants
mu = 3.98574405096E14
g = 9.81
R_E = 6.3781E6
J2 = 1.08262668E-3

P = 10E3
eta = 0.65
Isp = 3300.0
m0 = 1200.0

t_max = 86400.0 * 365 * 2

E_en = math.pi
E_ex = -math.pi

oe_i = np.array([R_E + 200000.0, 0, math.radians(28.5), math.radians(0.0), math.radians(0.0), math.radians(0.0)])
oe_f = np.array([R_E + 210000.0, 0, math.radians(28.5), math.radians(0.0), math.radians(0.0), math.radians(0.0)])

dm = (2 * eta * P)/((g * Isp)**2)
delta_t = 3600.0

#initialize model
traj = GEKKO()
nt = 200
traj.time = np.linspace(0, 1, nt)

#state variables
a = traj.SV(value = oe_i[0], lb = R_E, ub = oe_f[0] + 5000.0, name = 'sma')
e = traj.SV(value = oe_i[1], lb = 0, ub = 1, name = 'ecc')
i = traj.SV(value = oe_i[2], lb = math.radians(-90), ub = math.radians(90), name = 'inc')
Om = traj.SV(value = oe_i[3], lb = math.radians(-180), ub = math.radians(180), name = 'raan')
om = traj.SV(value = oe_i[4], lb = math.radians(-180), ub = math.radians(180), name = 'ap')
nu = traj.SV(value = oe_i[5], lb = math.radians(-180), ub = math.radians(180), name = 'ta')
m = traj.SV(value = m0, lb = 1000.0, ub = m0, name = 'mass')
t = traj.SV(value = 0.0, lb = 0.0, ub = t_max, name = 'time')

traj.periodic(i)
traj.periodic(Om)
traj.periodic(om)
traj.periodic(nu)

q = np.zeros(nt)
q[-1] = 1.0
final = traj.Param(value = q)

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

#manipulating variables and initial guesses
al_a = traj.MV(value = -1.0, lb = -2.0, ub = 2.0, name = 'al_a')
al_a.STATUS = 1
l_e = traj.MV(value = 0.001, lb = 0.0, ub = 1.0E6, name = 'l_e')
l_e.STATUS = 1
l_i = traj.MV(value = 1.0, lb = 0.0, ub = 1.0E6, name = 'l_i')
l_i.STATUS = 1

#equations
p = a * (1 - e**2)
r = p/(1 + e * traj.cos(nu))
rv = (np.array([r * traj.cos(nu), r * traj.sin(nu), 0]))
vv = (np.array([-traj.sin(nu) * traj.sqrt(mu/p), (e + traj.cos(nu)) * traj.sqrt(mu/p), 0]))
cO = traj.cos(Om)
sO = traj.sin(Om)
co = traj.cos(om)
so = traj.sin(om)
ci = traj.cos(i)
si = traj.sin(i)
R = np.array([[(cO * co - sO * so * ci), (-cO * so - sO * co * ci), (sO * si)], [(sO * co + cO * so * ci), (-sO * so + cO * co * ci), (-cO * si)], [(so * si), (co * si), ci]])
ri = np.array([R[0][0] * rv[0] + R[0][1] * rv[1] + R[0][2] * rv[2], R[1][0] * rv[0] + R[1][1] * rv[1] + R[1][2] * rv[2], R[2][0] * rv[0] + R[2][1] * rv[1] + R[2][2] * rv[2]])
vi = np.array([R[0][0] * vv[0] + R[0][1] * vv[1] + R[0][2] * vv[2], R[1][0] * vv[0] + R[1][1] * vv[1] + R[1][2] * vv[2], R[2][0] * vv[0] + R[2][1] * vv[1] + R[2][2] * vv[2]])
r = traj.sqrt(ri[0]**2 + ri[1]**2 + ri[2]**2)
v = traj.sqrt(vi[0]**2 + vi[1]**2 + vi[2]**2)
hi = np.cross(ri, vi)
h = traj.sqrt(hi[0]**2 + hi[1]**2 + hi[2]**2)
s_a = traj.Intermediate((-l_e * (r/a) * traj.sin(nu))/(traj.sqrt(4 * ((al_a * (a * v**2)/mu) + l_e * (e + traj.cos(nu)))**2 + l_e**2 * (r**2/a**2) * (traj.sin(nu))**2)))
c_a = traj.Intermediate((-2 * (((al_a * a * v**2)/mu) + l_e * (e + traj.cos(nu))))/(traj.sqrt(4 * ((al_a * a * v**2)/mu + l_e * (e + traj.cos(nu)))**2 + l_e**2 * (r**2/a**2) * (traj.sin(nu))**2)))
s_b = traj.Intermediate((-l_i * ((r * v)/h) * traj.cos(om + nu))/(traj.sqrt(l_i**2 * ((r**2 * v**2)/h**2) * (traj.cos(om + nu))**2 + ((4 * al_a**2 * a**2 * v**4)/mu**2) * c_a**2 + l_e**2 * ((2 * (e + traj.cos(nu)) * c_a + (r/a) * traj.sin(nu) * s_a)**2))))
c_b = traj.Intermediate((((-al_a * 2 * a * v**2)/mu) * c_a - l_e * (2 * (e + traj.cos(nu)) * c_a + (r/a) * traj.sin(nu) * s_a))/(traj.sqrt(l_i**2 * ((r**2 * v**2)/h**2) * (traj.cos(om + nu))**2 + ((4 * al_a**2 * a**2 * v**4)/mu**2) * c_a**2 + l_e**2 * ((2 * (e + traj.cos(nu)) * c_a + (r/a) * traj.sin(nu) * s_a)**2))))
a_T = (2 * eta * P)/(m * g * Isp)
a_n = a_T * s_a * c_b
a_t = a_T * c_a * c_b
a_h = a_T * s_b
n = traj.sqrt(mu/a**3)
Om_J2 = ((-3 * n * R_E**2 * J2)/(2 * a**2 * (1 - e**2)**2)) * traj.cos(i)
om_J2 = ((3 * n * R_E**2 * J2)/(4 * a**2 * (1 - e**2)**2)) * (4 - 5 * (traj.sin(i))**2)
dt_dE = r/(n * a)
Tp = (2 * math.pi/traj.sqrt(mu)) * a**(3/2)

#deltas
tp = traj.if3(t - tf, 1, 0)
traj.Equation(tp * Tp * a.dt() == (a_t * (2 * v * a**2)/mu) * delta_t * dt_dE)
traj.Equation(tp * Tp * e.dt() == ((1/v) * (2 * (e + traj.cos(nu)) * a_t + (r/a) * a_n * traj.sin(nu))) * delta_t * dt_dE)
traj.Equation(tp * Tp * i.dt() == ((r/h) * a_h * traj.cos(om + nu)) * delta_t * dt_dE)
traj.Equation(tp * Tp * Om.dt() == ((r/(h * traj.sin(i))) * a_h * traj.sin(om + nu) + Om_J2) * delta_t * dt_dE)
traj.Equation(tp * Tp * om.dt() == ((1/(e * v)) * (2 * a_t * traj.sin(nu) - (2 * e + (r/a) * traj.cos(nu)) * a_n) - (r/(h * traj.sin(i))) * a_h * traj.sin(om + nu) * traj.cos(i) + om_J2) * delta_t * dt_dE)
traj.Equation(tp * nu.dt() == (traj.acos((traj.cos((1/dt_dE) * delta_t) - e)/(1 - e * traj.cos((1/dt_dE) * delta_t))) - nu))
traj.Equation(tp * Tp * m.dt() == ((-2 * eta * P)/((g * Isp)**2)) * delta_t * dt_dE)
traj.Equation(t.dt() == delta_t)
traj.Equation(a * final == oe_f[0])
traj.Equation(e * final == oe_f[1])
traj.Equation(i * final == oe_f[2])
traj.Equation(Om * final == oe_f[3])
traj.Equation(om * final == oe_f[4])
traj.Equation(nu * final == oe_f[5])

#solve
traj.Obj(tf)
traj.options.IMODE = 6
traj.options.SOLVER = 3
traj.options.MAX_ITER = 15000
traj.options.RTOL = 1e-6
traj.options.OTOL = 1e-6
#traj.open_folder()
traj.solve()
print('Optimal time: ' + str(tf.value[0]))
traj.solve()
#traj.open_folder(infeasibilities.txt)

修复方案 #

1. 修正目标函数 #

FV变量初始值不能依赖状态变量，将tf改为固定初始值：

tf = traj.FV(value=3600*24*7, lb=0.0, ub=t_max)
tf.STATUS = 1

2. 删除冲突的周期性约束 #

移除traj.periodic(i)、traj.periodic(Om)、traj.periodic(om)、traj.periodic(nu)，因为最终轨道参数是固定目标值，无需与初始值周期一致。

3. 替换真近点角(nu)的微分方程 #

原方程数值不稳定，改用开普勒方程推导的关系：

M = traj.Intermediate(traj.sqrt(mu/a**3)*t)
E = traj.Intermediate(M + e*traj.sin(M))
traj.Equation(nu == 2*traj.atan(traj.sqrt((1+e)/(1-e))*traj.tan(E/2)))

删除原有的traj.Equation(tp * nu.dt() == ...)。

4. 调整MV的上下界与初始值 #

缩小l_e、l_i的上界，调整初始值为无推力状态：

al_a = traj.MV(value=0.0, lb=-2.0, ub=2.0, name='al_a')
al_a.STATUS = 1
l_e = traj.MV(value=0.0, lb=0.0, ub=1000.0, name='l_e')
l_e.STATUS = 1
l_i = traj.MV(value=0.0, lb=0.0, ub=1000.0, name='l_i')
l_i.STATUS = 1

5. 修正质量微分方程 #

使用正确的火箭质量消耗公式（删除原方程，替换为）：

a_T = (eta * P)/(m * g * Isp)  # 移除原公式中的多余因子2
traj.Equation(m.dt() == -a_T / g)

6. 优化时间缩放逻辑 #

将traj.Equation(t.dt() == delta_t)改为：

traj.Equation(t.dt() == tf)

让求解器自动处理时间缩放，避免固定delta_t带来的数值刚性。

7. 配置IPOPT求解器参数 #

添加参数提升数值稳定性：

traj.solver_options = ['max_iter 10000',
                       'tol 1e-6',
                       'constr_viol_tol 1e-6',
                       'compl_inf_tol 1e-6',
                       'mu_strategy adaptive',
                       'linear_solver ma57']

说明 #

以上修改针对数值病态、约束冲突、物理公式错误三个核心问题，修复后可大幅提升求解收敛概率。

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