python - GEKKO如何解决一个最优的月球软着陆问题？

Question

我试图通过 GEKKO 解决一个二维最优月球软着陆问题。我假设月亮没有旋转。飞船应该在月球表面轻轻着陆，即最终垂直速度v和最终的水平速度 u应该为零，最终高度r应该是月球的半径。
这个问题可以说明如下:

状态和控制变量及方程如下(控制变量为推力 F 和姿态角 φ ):

我构建了 燃油优化问题如下(控制变量 φ 写成 f ):

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import time
import os

r0=1753000
v0=1675

time_start = time.time()

model = GEKKO()
nt = 501
model.time = np.linspace(0,1,nt)

# optimize final time
tf = model.FV(value=1.0,lb=0.1,ub=1000.0)
tf.STATUS = 1

# controls
F = model.MV(value=0,lb=0,ub=2000)
F.STATUS = 1
f = model.MV(value=-0.5*np.pi,lb=-0.5*np.pi,ub=0.5*np.pi)
f.STATUS = 1

# state variables
r = model.Var(value=r0,lb=1738000) # height
s = model.Var(value=0) # tru anomaly
v = model.Var(value=0) # vertical velocity
u = model.Var(value=v0) # horizional velocity
m = model.Var(value=600,lb=0,ub=600) # mass

# constants
mu = 4.90275*10**(12) # lunar gravitational constant
Isp = 300*9.8 # specific impulse

# Equations
model.Equation( r.dt() == tf * v )
model.Equation( s.dt() == tf * u/r )
model.Equation( v.dt() == tf * (F/m*model.cos(f)-mu/(r**2)+u**2/r) )
model.Equation( u.dt() == tf * (F/m*model.sin(f)-u*v/r) )
model.Equation( m.dt() == tf * (-F/Isp) )

# terminal constraints
_finalMask = np.zeros(nt)
_finalMask[-1] = 1.0
finalMask = model.Param(value=_finalMask)
model.Equation(v*finalMask>=0)
model.Equation(v*finalMask<=0.5)
model.fix_final(r,val=1738000)
model.fix_final(u,val=0)

model.Obj(-m*finalMask) # Objective function to be minimized

model.options.IMODE = 6
model.solver_options = ['max_iter 5000']
model.solve() # solve

time_end=time.time()
print('Calculation Time: ',time_end-time_start)

# scaled time
print('Landing Time: ' + str(tf.value[0]))
tm = np.linspace(0,tf.value[0],nt)
finaltime = tm[-1]

# PLOT
fig = plt.figure(1)
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 3)
ax4 = fig.add_subplot(2, 3, 4)
ax5 = fig.add_subplot(2, 3, 5)
ax6 = fig.add_subplot(2, 3, 6)
ax1.plot(tm,r.value,'k-',label=r'$r$')
ax1.set_xlabel('Time')
ax1.set_ylabel('r')
ax1.set_xlim(0,finaltime)
ax2.plot(tm,v.value,'b-',label=r'$v$')
ax2.set_xlabel('Time')
ax2.set_ylabel('v')
ax2.set_xlim(0,finaltime)
ax3.plot(tm,u.value,'g-',label=r'$w$')
ax3.set_xlabel('Time')
ax3.set_ylabel('u')
ax3.set_xlim(0,finaltime)
ax4.plot(tm,m.value,'y-',label=r'$m$')
ax4.set_xlabel('Time')
ax4.set_ylabel('m')
ax4.set_xlim(0,finaltime)
ax5.plot(tm,f.value,'c-',label=r'$f$')
ax5.set_xlabel('Time')
ax5.set_ylabel('f')
ax5.set_xlim(0,finaltime)
ax6.plot(tm,F.value,'r-',label=r'$F$')
ax6.set_xlabel('Time')
ax6.set_ylabel('F')
ax6.set_xlim(0,finaltime)
plt.tight_layout()
plt.show()

它可以找到一个与其他人的研究非常相似的解决方案。

但是 f出现了急转弯.这是 Not Acceptable ，因为角度 φ 应该不断地改变。
另外，我尝试缩放功能

scale = 1e-6
model.Equation( r.dt() == tf * v )
model.Equation( r*s.dt()*scale == tf * u*scale )
model.Equation( m*(r**2)*v.dt()*scale**2 == tf * ((r**2)*F*model.cos(f)-mu*m+(u**2)*r*m)*(scale**2) )
model.Equation( m*r*u.dt()*scale == tf * (F*r*model.sin(f)-u*v*m)*scale )
model.Equation( Isp*m.dt() == tf * (-F) )

但失败了

Solution Not Found
EXIT: Converged to a point of local infeasibility. Problem may be infeasible.

为了获得顺利 f ，我将第二个控制变量改为角加速度 a ，状态方程变为

代码变成了:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import time
import os

r0=1753000
v0=1675

time_start = time.time()

model = GEKKO() # remote=False
nt = 501
model.time = np.linspace(0,1,nt)

tf = model.FV(value=1.0,lb=0.1,ub=1000.0)
tf.STATUS = 1

# controls
F = model.MV(value=0,lb=0,ub=2000)
F.STATUS = 1
a = model.MV(value=0,lb=-0.5*np.pi/180,ub=0.5*np.pi/180)
a.STATUS = 1

# state variables
r = model.Var(value=r0,lb=1738000) # height
s = model.Var(value=0) # tru anomaly
v = model.Var(value=0) # vertical velocity
u = model.Var(value=v0) # horizional velocity
f = model.Var(value=-0.5*np.pi) # angle
w = model.Var(value=0, lb=-10*np.pi/180,ub=10*np.pi/180) # angular velocity
m = model.Var(value=600,lb=0,ub=600) # mass

# constants
mu = 4.90275*10**(12) # lunar gravitational constant
Isp = 300*9.8 # specific impulse

# Equations
model.Equation( r.dt() == tf * v )
model.Equation( s.dt() == tf * u/r )
model.Equation( v.dt() == tf * (F/m*model.cos(f)-mu/(r**2)+u**2/r) )
model.Equation( u.dt() == tf * (F/m*model.sin(f)-u*v/r) )
model.Equation( f.dt() == tf * (w - u/r) ) # ---newly added
model.Equation( w.dt() == tf * a ) # ---newly added
model.Equation( m.dt() == tf * (-F/Isp) )

# terminal constraints
_finalMask = np.zeros(nt)
_finalMask[-1] = 1.0
finalMask = model.Param(value=_finalMask)
model.Equation(v*finalMask>=0)
model.Equation(v*finalMask<=0.5)
model.Equation(f*finalMask>=-5*np.pi/180) # ***newly added
model.Equation(f*finalMask<=5*np.pi/180) # ***newly added
model.fix_final(r,val=1738000)
model.fix_final(u,val=0)

model.Obj(-m*finalMask) # Objective function to be minimized

model.options.IMODE = 6
model.solver_options = ['max_iter 99999']
model.solve() # solve

time_end=time.time()
print('Calculation Time: ',time_end-time_start)

# scaled time
print('Landing Time: ' + str(tf.value[0]))
tm = np.linspace(0,tf.value[0],nt)
finaltime = tm[-1]

# PLOT
fig = plt.figure(1)
ax1 = fig.add_subplot(2, 4, 1)
ax2 = fig.add_subplot(2, 4, 2)
ax3 = fig.add_subplot(2, 4, 3)
ax4 = fig.add_subplot(2, 4, 4)
ax5 = fig.add_subplot(2, 4, 5)
ax6 = fig.add_subplot(2, 4, 6)
ax7 = fig.add_subplot(2, 4, 7)
ax8 = fig.add_subplot(2, 4, 8)
ax1.plot(tm,r.value,'k-',label=r'$r$')
ax1.set_xlabel('Time')
ax1.set_ylabel('r')
ax1.set_xlim(0,finaltime)
ax2.plot(tm,v.value,'b-',label=r'$v$')
ax2.set_xlabel('Time')
ax2.set_ylabel('v')
ax2.set_xlim(0,finaltime)
ax3.plot(tm,u.value,'g-',label=r'$u$')
ax3.set_xlabel('Time')
ax3.set_ylabel('u')
ax3.set_xlim(0,finaltime)
ax4.plot(tm,m.value,'y-',label=r'$m$')
ax4.set_xlabel('Time')
ax4.set_ylabel('m')
ax4.set_xlim(0,finaltime)
ax5.plot(tm,F.value,'r-',label=r'$F$')
ax5.set_xlabel('Time')
ax5.set_ylabel('F')
ax5.set_xlim(0,finaltime)
ax6.plot(tm,f.value,'k-',label=r'$f$')
ax6.set_xlabel('Time')
ax6.set_ylabel('f')
ax6.set_xlim(0,finaltime)
ax7.plot(tm,w.value,'b-',label=r'$w$')
ax7.set_xlabel('Time')
ax7.set_ylabel('w')
ax7.set_xlim(0,finaltime)
ax8.plot(tm,a.value,'r-',label=r'$a$')
ax8.set_xlabel('Time')
ax8.set_ylabel('a')
ax8.set_xlim(0,finaltime)
plt.tight_layout()
plt.show()

但是，我已经运行了很多小时的代码，但没有得到解决方案或结果。更具体地说，它在某些迭代时停止:

我尝试使用软终端约束但也没有得到任何结果:

model.Equation(r*finalMask>=1738000)
model.Equation(r*finalMask<=1738001)
model.Equation(v*finalMask>=0)
model.Equation(v*finalMask<=0.5)
model.Equation(u*finalMask>=-0.5)
model.Equation(u*finalMask<=0.5)
model.Equation(f*finalMask>=-5*np.pi/180)
model.Equation(f*finalMask<=5*np.pi/180)

Answer 1

我无法运行您的代码，因为 w未定义但可能已更改为 u .此外，有很多地方需要分号来分隔同一行的语句。我试图获得一个可以运行但尚未成功的版本。您能否验证发布的代码是否运行以显示您描述的问题，或者发布成功运行的更新代码？
我修改了方程式以与您的问题陈述保持一致，并避免除以零可能出现的潜在问题。我还包括了一个比例因子，因为一些方程有 r**2这是一个非常大的数字。

# equations
scale = 1e-6
model.Equation( r.dt() == tf * v )
model.Equation( r * s.dt()*scale == tf * u*scale )
model.Equation( m * (r**2) * v.dt() * scale**2 \
                == tf * ((r**2)*F*model.cos(f)-mu+(u**2)*(r**3)) *(scale**2))
model.Equation( m * r * u.dt()*scale == tf * (-u*v*m+F*model.cos(f))*scale )
model.Equation( Isp * m.dt() == -tf * F )

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

model = GEKKO() # initialize gekko
nt = 501
model.time = np.linspace(0,1,nt)

# optimize final time
tf = model.FV(value=100.0,lb=1.0,ub=1000.0)
tf.STATUS = 1

# controls
F = model.MV(value=0,lb=0,ub=1500)
F.STATUS = 1
f = model.MV(value=0,lb=-0.5*np.pi,ub=0.5*np.pi)
f.STATUS = 1

# state variables
# height = 1753 km
r = model.Var(value=1753000,lb=1738000) 
# vertical velocity = 0 m/s
v = model.Var(value=0)
# azimuth = 0 rad
s = model.Var(value=0)
# azimuth angular velocity rad/s
u = model.Var(value=9.65e-4)
# mass = 600 kg
m = model.Var(value=600,lb=0,ub=600) 

# constants
mu = 4.90275e12 # lunar gravitational constant
Isp = 300*9.8 # specific impulse

# equations
scale = 1e-6
model.Equation( r.dt() == tf * v )
model.Equation( r * s.dt()*scale == tf * u*scale )
model.Equation( m * (r**2) * v.dt() * scale**2 \
                == tf * ((r**2)*F*model.cos(f)-mu+(u**2)*(r**3)) *(scale**2))
model.Equation( m * r * u.dt()*scale == tf * (-u*v*m+F*model.cos(f))*scale )
model.Equation( Isp * m.dt() == -tf * F )

# terminal conditions
_finalMask = np.zeros(nt)
_finalMask[-1] = 1.0
finalMask = model.Param(value=_finalMask)

# Terminal Constraint 1
if False:
    model.Equation((r-1738000)*finalMask==0)
    model.Equation(v*finalMask==0)
    model.Equation(u*finalMask==0)

# Terminal Constraint 2
if False:
    model.Equation((r-1738000)*finalMask<=0)
    model.Equation(v*finalMask<=0)
    model.Equation(u*finalMask<=0)

if True:
    # terminal constraints
    _finalMask = np.zeros(nt)
    _finalMask[-1] = 1.0
    finalMask = model.Param(value=_finalMask)
    model.Equation(v*finalMask>=0)
    model.Equation(v*finalMask<=0.2)
    model.fix_final(r,val=1738000)
    model.fix_final(u,val=0)

model.Minimize(tf)
model.options.IMODE = 6
model.solve()

print('Final Time: ' + str(tf.value[0]))
tm = np.linspace(0,tf.value[0],nt)

# PLOT
fig = plt.figure(1)
ax1 = fig.add_subplot(2, 3, 1)
ax2 = fig.add_subplot(2, 3, 2)
ax3 = fig.add_subplot(2, 3, 3)
ax4 = fig.add_subplot(2, 3, 4)
ax5 = fig.add_subplot(2, 3, 5)
ax6 = fig.add_subplot(2, 3, 6)
ax1.plot(tm,r.value,'k-',label=r'$r$')
ax1.set_xlabel('Time'); ax1.set_ylabel('r')
ax2.plot(tm,v.value,'b-',label=r'$v$')
ax2.set_xlabel('Time'); ax2.set_ylabel('v')
ax3.plot(tm,s.value,'g-',label=r'$s$')
ax3.set_xlabel('Time'); ax3.set_ylabel('s')
ax4.plot(tm,f.value,'y-',label=r'$f$')
ax4.set_xlabel('Time'); ax4.set_ylabel('f')
ax5.plot(tm,F.value,'c-',label=r'$F$')
ax5.set_xlabel('Time'); ax5.set_ylabel('F')
ax6.plot(tm,m.value,'r-',label=r'$m$')
ax6.set_xlabel('Time'); ax6.set_ylabel('m')
plt.tight_layout()
plt.show()

如果所有变量都在 1 附近并且方程残差从大约相同的水平开始，则数值解更容易求解。 Gekko 根据初始条件自动缩放变量。
正如您所观察到的，与精确地达到精确的最终条件相比，在小范围内找到解决方案会更容易。另一件有帮助的事情是包含一个目标(软约束)和平等约束(硬约束)。这种组合通常更适合这些类型的终端约束问题，最终时间可调整。
您可能还想使用指定的最终时间进行初始化，然后再次求解。您可以关闭自由度，例如:

# solve first time without adjustable tf
tf.STATUS = 0
model.solve()

# solve again with tf adjustable
tf.STATUS = 1
# don't advance the initial condition
model.options.TIME_SHIFT=0
model.solve()

另一种选择是用 F.UPPER=1500 解决然后尝试 F.UPPER=2000为另一个解决。 Jennings optimal control problem是另一个可能有帮助的例子。