利用IPOPT GEKKO对航天器轨迹进行优化

Question

因此，我试图解决IPOPT (python)中的一个非线性微分优化问题。我需要最小化时间的航天器轨道，其中的运动方程是由圆形限制的三体动力学描述。轨道是从地球到L1之间的太阳-地球系统.我有固定的初始和端点位置和速度坐标，我需要积分微分方程，这样从地球到最终位置的时间是最小的。在CRTBP的动力学中，我在x，y，z，方向加上无量纲推力，是这个优化问题的控制变量。这导致路径约束，在每个推力段，总推力大小必须小于或等于每个推力分量的平方之和(见下面的问题公式图)。

​

​

我对最优控制相关的问题和软件，如IPOPT和GEKKO都很陌生。我试着用IPOPT作为优化软件在GEKKO中编写这个问题，但它似乎从未收敛。我尝试遵循这个轨迹优化示例，特别是2D部分：optimisation.html。

显然，在我的问题上，我有更复杂的动力学以及三维轨迹优化。但考虑到我有固定的起点和终点，并且给出了每个变量的初始猜测(接近解)，最后，这只是一个在A和B之间寻找最佳路径的问题，给出了三体动力学。虽然我一直在犯这个错误，

出口:收敛到局部不可行点。问题可能是不可行的。

这是我的代码：

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
#initialize model
m = GEKKO()

# optional solver settings with APOPT
Nsim = 100 #steps when thrust is constant
# tf = 12 days /tau_star nondim

m.time = np.linspace(0, 0.206, Nsim)
#constants
mu = 3.003481e-6
l_star = 1.49597871e11
G_c = 6.67408 * 10**-11
m_star = 1.9885e30
tau_st = (l_star**3/G_c/m_star)**0.5 # circa 58 days
Isp = 3000 #  impulse per second
g0 = 9.80665* 10**-3 # km/s^2 => g = 9.8
T_max = 0.1 # {N, 100 mN}
m0 = 500 #kg
f = (T_max*(tau_st**2)/(l_star * m0)) #maximal nondim thrust
sub_L1 = 0.0158
x_final = sub_L1 - (1-mu)

#manipulating variables thrust components in x,y,z direction
#initial values represent full-thrust in x-direction as well as downwards
ux = m.MV(0.99*f, lb=-f, ub=f)
ux.STATUS = 1
uy = m.MV(0, lb=-f, ub=f)
uy.STATUS = 1
uz = m.MV(m.sqrt(1 - 0.99**2)*f, lb=-f, ub=f)
uz.STATUS = 1

#variables + initial guesses??
x1 = m.Var(value=x_final, lb=0, ub=1)
x2 = m.Var(value= -0.0023, lb=0, ub=1)
x3 = m.Var(value= 0, lb=0, ub=1)
x4 = m.Var(value= 0.3, lb=-1, ub=1)
x5 = m.Var(value= 0.3, lb=-0.5, ub=0.5)
x6 = m.Var(value= 0.3, lb=-0.5, ub=0.5)
#constraint on objetive function
# guess value for tf is 9 days meaning for an optimal path less time is needed to get to x_final
tf = m.FV(value = 0.15, lb=0.03, ub=0.5) # time upper and lower constraints
tf.STATUS = 1


#defining r1 and r2 as equations to be solved implicity together
a = m.Var(value=0.99)
b = m.Var(value=0.004)
m.Equation(((x1 - mu)**2 + x2**2 + x3**2)**(3/2) == a) #r1
m.Equation(((x1 + 1 - mu)**2 + x2**2 + x3**2)**(3/2) == b) #r2
#dynamics of three-body-problem
m.Equation( x1.dt() == x4*tf)
m.Equation( x2.dt() == x5*tf)
m.Equation( x3.dt() == x6*tf)
m.Equation( x4.dt() == tf*(2*x5 + x1 - ((1-mu)*(x1-mu)/a) - (mu*(x1+1-mu))/b + ux))
m.Equation( x5.dt() == tf*(-2*x4 + x2 - ((1-mu)*x2/a) - (mu*x2)/b + uy))
m.Equation( x6.dt() == tf*(-(1-mu) * x3 / a - (mu*x3)/b + uz))

#path constraints
m.Equation(x1 <= -0.1)
eq = m.Param(value=f)
m.Equation(ux**2 + uy**2 + uz**2 <= eq**2)
#starting constraints, starting x-value is from earth's escape
m.fix(x1, pos=0,val=-0.9952)
m.fix(x2, pos=0,val=-0.0023)
m.fix(x3, pos=0,val=-0.0010)
#m.fix(x4, pos=0,val=1.3890)
m.fix(x4, pos=0,val=0.002)
#m.fix(x5, pos=0,val=1.0585)
m.fix(x5, pos=0,val=0.001)
m.fix(x6, pos=0,val=0.0273)
#boundary constraints
m.fix(x1, pos=len(m.time)-1,val=x_final) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x2, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x3, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x4, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x5, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x6, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.Obj(tf) # minimize final time
m.options.IMODE = 6 # non-linar model
#m.options.NODES = 3 # collocation nodes
m.options.SOLVER = 3 # solver (IPOPT)
m.options.MAX_ITER = 15000
m.options.RTOL = 1e-3
m.options.OTOL = 1e-3
m.solve() # Solve
print('Optimal time: ' + str(tf.value[0]))
m.solve(disp=True)
m.open_folder(infeasibilities.txt)

Answer 1

一个不可行的条件是方程m.Equation(x1 <= -0.1)，但变量定义为0的下界，1的上界与x1 = m.Var(value=x_final, lb=0, ub=1)。尝试删除状态变量上的不等式约束，直到问题变得可行为止。依次添加约束，直到再次变得不可行为止，以确定导致不可行解决方案的约束(S)。

有时问题是可行的，但求解者无法找到解决方案或陷入不可行区域。对于求解者来说，端点硬约束可能是一种挑战：

m.fix(x1, pos=len(m.time)-1,val=x_final) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x2, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x3, pos=len(m.time)-1,val=0) # final destination is sub-L1 in the Sun-Earth-system
m.fix(x4, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x5, pos=len(m.time)-1,val=0.0) # stationary in sub-L1
m.fix(x6, pos=len(m.time)-1,val=0.0) # stationary in sub-L1

一种策略是将软(客观)约束与硬约束相结合，以帮助指导求解者。看看倒立摆问题所采取的策略。

​

​

它通过以下组合更好地解决了问题：

# soft constraints
m.Minimize(final*ya**2)
m.Minimize(final*va**2)
m.Minimize(final*theta_a**2)
m.Minimize(final*qa**2)

# hard constraints
m.fix(ya,pos=end_loc,val=0.0)
m.fix(va,pos=end_loc,val=0.0)
m.fix(theta_a,pos=end_loc,val=0.0)
m.fix(qa,pos=end_loc,val=0.0)

另一个可能有帮助的问题是黑尔飞机优化问题。

• 马丁，R.A.，盖茨，N.，宁，A.，赫登仁，J.D.，“在站台约束下高空太阳能飞机飞行轨迹的动态优化”，“制导、控制和动力学杂志”，2018年，doi: 10.2514/1.G003737。预印 _ 文章 _ 源代码

​

​