python - GEKKO如何解决一个最优的月球软着陆问题?
问题描述
我试图通过 GEKKO 解决 2D 最优月球软着陆问题。我假设月亮没有旋转。航天器应在月球表面软着陆,即最终垂直速度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
。这是不可接受的,因为角度φ应该不断变化。
另外,我尝试扩展功能
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)
解决方案
我无法让您的代码运行,因为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 最优控制问题是另一个可能有帮助的例子。
推荐阅读
- azure-active-directory - msonline 模块 - where 语句在 azure 函数 powershell 中不起作用
- javascript - 使用输入文本框值搜索电影评论 API AJAX 方法 - JS
- solr - Solr中具有相同名称的多个字段
- python - z3 power mod 约束抛出 Z3Exception
- visual-studio-code - 使用 **watch** 在 VSCode 中调试独立的 Blazor WebAssembly 应用程序
- python - 使用理解遍历字典以将所有日期时间值转换为 MM/DD/YYYY 字符串
- python-3.x - Pandas DataFrame to_dict 和 from_dict 会丢失列顺序吗?
- java - 更新操作但没有更新行时的 HttpStatus 代码
- google-cloud-platform - 我的 GCP 项目仪表板中的未知项目
- javascript - 使用javascript变量html将链接放入href