文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
社区首页 >问答首页 >如何利用Gekko求解可再入飞行器的最优控制

问如何利用Gekko求解可再入飞行器的最优控制
EN

Stack Overflow用户
提问于 2022-01-13 17:20:01
回答 2查看 111关注 0票数 3

我正在寻求最优控制(aoa和岸角)，以最大限度地跨越距离的穿梭机型再入飞行器使用Gekko。下面是我当前的代码，我得到了一个“未找到的解决方案”和"EXIT:迭代次数超过的最大值“。模拟假设点质量与一个非旋转的地球框架.EOMS为6种耦合的非线性ODEs。我尝试过使用不同的求解器，实现/删除状态和控制约束，增加最大迭代次数等。我对我在Gekko中的问题的设置和实现没有信心，我希望得到一些关于我下一步可以尝试什么的反馈。我试图遵循APMonitor示例11.带积分目标的最优控制、倒立摆最优控制和示例13.最优控制:最大限度地减少最后时间中的设置和布局。我正在寻找的解决方案如下。

任何帮助都是非常感谢的！

代码语言：javascript
复制
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math

pi = math.pi

########################
######FRONT MATTER######
########################
m = GEKKO() # initialize GEKKO
nt = 2501 #simulation time is 2500 seconds
tfin = 2500
m.time = np.linspace(0,tfin,nt) #time array

#==================#
#PARAMS
#==================#
Re  = m.Param(value = 6371203.92)           # radius of the earth, m
S   = m.Param(value = 249.9091776)          # vehicle surface area, m^2
cl0 = m.Param(value = -0.2070)              # coeff lift param 1
cl1 = m.Param(value = 1.6756)               # coeff lift param 2
cd0 = m.Param(value = 0.0785)               # coeff drag param 1
cd1 = m.Param(value = -0.3529)              # coeff drag param 2
cd2 = m.Param(value = 2.0400)               # coeff drag param 3
H   = m.Param(value = 7254.24)              # density scale height, m
rho0= m.Param(value = 1.225570827014494)    # sea level atmospheric density, kg/m^3
mu  = m.Param(value = 3.986031954093051e14) #earth gravitational param, m^3/s^2
mass= m.Param(value = 92079.2525560557)     #vehicle mass, kg

#===============================#
#BOUNDARY CONDITIONS
#===============================#
t0     = 0
alt0   = 79248
rad0   = alt0+Re
altf   = 24384
radf   = altf+Re
lon0   = 0
lat0   = 0
speed0 = +7802.88
speedf = +762
fpa0   = -1*pi/180
fpaf   = -5*pi/180
azi0   = +90*pi/180
azif   = -90*pi/180

#===============================#
#LIMITS ON VARIABLES
#===============================#
tfMin = 0;              tfMax = 3000;
radMin = Re;            radMax = rad0;
lonMin = -pi;           lonMax = -lonMin;
latMin = -70*pi/180;    latMax = -latMin;
speedMin = 10;          speedMax = 45000;
fpaMin = -80*pi/180;    fpaMax =  80*pi/180;
aziMin = -180*pi/180;   aziMax =  180*pi/180;
aoaMin = -90*pi/180;    aoaMax = -aoaMin;
bankMin = -90*pi/180;   bankMax =   1*pi/180;

#===============================#
#VARIABLES
#===============================#

#state variables and bounds
rad = m.Var(value=rad0, lb=radMin, ub=radMax)       # radius, m
lon = m.Var(value=lon0, lb=lonMin, ub=lonMax)       # longitude, rad
lat = m.Var(value=lat0, lb=latMin, ub=latMax)       # latitude, rad
vel = m.Var(value=speed0, lb=speedMin, ub=speedMax) # velocity, m/sec
fpa = m.Var(value=fpa0, lb=fpaMin, ub=fpaMax)       # flight path angle, rad
azi = m.Var(value=azi0, lb=aziMin, ub=aziMax)       # azimuth angle, rad

#control variables
aoa     = m.MV(value=-20, lb=aoaMin, ub=aoaMax)    # angle of attack, rad
aoa.STATUS = 1
aoa.DCOST = 1e-2
bank    = m.MV(value=0, lb=bankMin, ub=bankMax)    # bank angle, rad
bank.STATUS = 1
bank.DCOST = 1e-2

#===============================#
#INTERMEDIATE VARIABLES
#===============================#
altitude = m.Intermediate(rad - Re)
CD       = m.Intermediate(cd0+cd1*aoa+cd2*aoa**2)
rho      = m.Intermediate(rho0*m.exp(-altitude/H))
CL       = m.Intermediate(cl0+cl1*aoa)
q        = m.Intermediate(0.5*rho*vel**2)
D        = m.Intermediate(q*S*CD/mass)
L        = m.Intermediate(q*S*CL/mass)
gravity  = m.Intermediate(mu/rad**2)

#===============================#
#EOMS
#===============================#

p = np.zeros(nt) # mark final time point
p[-1] = 1.0
final = m.Param(value=p)

m.Equation(rad.dt() == vel*m.sin(fpa))
m.Equation((rad*m.cos(lat))*lon.dt() == vel*m.cos(fpa)*m.sin(azi))
m.Equation(rad*lat.dt() == vel*m.cos(fpa)*m.cos(azi))
m.Equation(vel.dt() == -D-gravity*m.sin(fpa))
m.Equation(vel*fpa.dt() == (L*m.cos(bank)-m.cos(fpa)*(gravity-vel**2/rad)))
m.Equation(vel*azi.dt() == (L*m.sin(bank)/m.cos(fpa)+vel**2*m.cos(fpa)*m.sin(azi)*m.tan(lat)/rad))

#===============================#
#OPTIMIZATION SOLVER
#===============================#
m.Obj(-lat*final)


# m.options.SOLVER = 3
# m.options.IMODE = 6
# m.solve(disp=True)

m.options.MAX_ITER = 500
m.options.IMODE = 6
# m.options.NODES = 3
# m.options.MV_TYPE = 1
m.options.SOLVER = 3
# m.open_folder()
m.solve()

对于高度、速度、AOA和岸角，解决方案应该如下所示：海拔高度 速度 AngleofAttack 银行

ode
nonlinear-optimization
gekko
python
EN

回答 2

Stack Overflow用户

回答已采纳

发布于 2022-01-14 04:17:31

通过减少最后的时间(max=0.04表示成功的解)并重新排列方程以避免可能的零除法，我得到了一个成功的解：

代码语言：javascript
复制
m.Equation(rad.dt() == vel*m.sin(fpa))
m.Equation((rad*m.cos(lat))*lon.dt() == vel*m.cos(fpa)*m.sin(azi))
m.Equation(rad*lat.dt() == vel*m.cos(fpa)*m.cos(azi))
m.Equation(vel.dt() == -D-gravity*m.sin(fpa))
m.Equation(vel*fpa.dt() == (L*m.cos(bank)-m.cos(fpa)*(gravity-vel**2/rad)))
m.Equation(m.cos(fpa)*rad*vel*azi.dt() == \
           (L*m.sin(bank)*rad+vel**2*(m.cos(fpa))**2*m.sin(azi)*m.tan(lat)))

从一小段时间开始，观察状态可以帮助解决模型的疑难问题。

​

​

代码语言：javascript
复制
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math

pi = math.pi

########################
######FRONT MATTER######
########################
m = GEKKO() # initialize GEKKO
nt = 101 # 2501 #simulation time is 2500 seconds
tfin = 0.04
m.time = np.linspace(0,tfin,nt) #time array

#==================#
#PARAMS
#==================#
Re  = m.Param(value = 6371203.92)           # radius of the earth, m
S   = m.Param(value = 249.9091776)          # vehicle surface area, m^2
cl0 = m.Param(value = -0.2070)              # coeff lift param 1
cl1 = m.Param(value = 1.6756)               # coeff lift param 2
cd0 = m.Param(value = 0.0785)               # coeff drag param 1
cd1 = m.Param(value = -0.3529)              # coeff drag param 2
cd2 = m.Param(value = 2.0400)               # coeff drag param 3
H   = m.Param(value = 7254.24)              # density scale height, m
rho0= m.Param(value = 1.225570827014494)    # sea level atmospheric density, kg/m^3
mu  = m.Param(value = 3.986031954093051e14) #earth gravitational param, m^3/s^2
mass= m.Param(value = 92079.2525560557)     #vehicle mass, kg

#===============================#
#BOUNDARY CONDITIONS
#===============================#
t0     = 0
alt0   = 79248
rad0   = alt0+Re
altf   = 24384
radf   = altf+Re
lon0   = 0
lat0   = 0
speed0 = +7802.88
speedf = +762
fpa0   = -1*pi/180
fpaf   = -5*pi/180
azi0   = +90*pi/180
azif   = -90*pi/180

#===============================#
#LIMITS ON VARIABLES
#===============================#
tfMin = 0;              tfMax = 3000;
radMin = Re;            radMax = rad0;
lonMin = -pi;           lonMax = -lonMin;
latMin = -70*pi/180;    latMax = -latMin;
speedMin = 10;          speedMax = 45000;
fpaMin = -80*pi/180;    fpaMax =  80*pi/180;
aziMin = -180*pi/180;   aziMax =  180*pi/180;
aoaMin = -90*pi/180;    aoaMax = -aoaMin;
bankMin = -90*pi/180;   bankMax =   1*pi/180;

#===============================#
#VARIABLES
#===============================#

#state variables and bounds
rad = m.Var(value=rad0, lb=radMin, ub=radMax)       # radius, m
lon = m.Var(value=lon0, lb=lonMin, ub=lonMax)       # longitude, rad
lat = m.Var(value=lat0, lb=latMin, ub=latMax)       # latitude, rad
vel = m.Var(value=speed0, lb=speedMin, ub=speedMax) # velocity, m/sec
fpa = m.Var(value=fpa0, lb=fpaMin, ub=fpaMax)       # flight path angle, rad
azi = m.Var(value=azi0, lb=aziMin, ub=aziMax)       # azimuth angle, rad

#control variables
aoa     = m.MV(value=-20, lb=aoaMin, ub=aoaMax)    # angle of attack, rad
aoa.STATUS = 1
bank    = m.MV(value=0, lb=bankMin, ub=bankMax)    # bank angle, rad
bank.STATUS = 1

#===============================#
#INTERMEDIATE VARIABLES
#===============================#
altitude = rad - Re
CD       = cd0+cd1*aoa+cd2*aoa**2
rho      = rho0*m.exp(-altitude/H)
CL       = cl0+cl1*aoa
q        = 0.5*rho*vel**2
D        = q*S*CD/mass
L        = q*S*CL/mass
gravity  = mu/rad**2

#===============================#
#EOMS
#===============================#

p = np.zeros(nt) # mark final time point
p[-1] = 1.0
final = m.Param(value=p)

m.Equation(rad.dt() == vel*m.sin(fpa))
m.Equation((rad*m.cos(lat))*lon.dt() == vel*m.cos(fpa)*m.sin(azi))
m.Equation(rad*lat.dt() == vel*m.cos(fpa)*m.cos(azi))
m.Equation(vel.dt() == -D-gravity*m.sin(fpa))
m.Equation(vel*fpa.dt() == (L*m.cos(bank)-m.cos(fpa)*(gravity-vel**2/rad)))
m.Equation(m.cos(fpa)*rad*vel*azi.dt() == \
           (L*m.sin(bank)*rad+vel**2*(m.cos(fpa))**2*m.sin(azi)*m.tan(lat)))

#===============================#
#OPTIMIZATION SOLVER
#===============================#
m.Maximize(lat*final)
m.options.SOLVER = 3
m.options.IMODE = 6
m.solve(disp=True)

plt.subplot(4,2,1)
plt.plot(m.time,rad.value,label='rad')
plt.legend()
plt.subplot(4,2,2)
plt.plot(m.time,lon.value,label='lon')
plt.legend()
plt.subplot(4,2,3)
plt.plot(m.time,lat.value,label='lat')
plt.legend()
plt.subplot(4,2,4)
plt.plot(m.time,vel.value,label='vel')
plt.legend()
plt.subplot(4,2,5)
plt.plot(m.time,fpa.value,label='fpa')
plt.legend()
plt.subplot(4,2,6)
plt.plot(m.time,azi.value,label='azi')
plt.legend()
plt.subplot(4,2,7)
plt.plot(m.time,aoa.value,label='aoa')
plt.xlabel('Time')
plt.legend()
plt.subplot(4,2,8)
plt.plot(m.time,bank.value,label='bank')
plt.xlabel('Time')
plt.legend()
plt.show()

一个建议是关闭自由度，以验证解决方案(STATUS=0)与较小的时间范围。还可能需要额外的约束来将-2_pi中的三角函数保持为2_pi区域。还有关于初始化具有挑战性的问题的其他信息：

  • Safdarnejad，S.M.，Hedengren，J.D.，Lewis，N.R.，Haseltine，E.，动态系统优化的初始化策略，Computers and Chemical，2015年，第78卷，第39至50页，DOI: 10.1016/j.compchemeng.2015.04.016。
票数 2
EN

Stack Overflow用户

发布于 2022-01-16 05:50:28

成功！！谢谢您的反馈和指点@。在吸收了我在这里学到的知识，并恢复了状态的界限之后，我成功地使用了Gekko来产生一个与以前的解决方案一致的答案。非常兴奋的Gekko和这里的潜力！

​

​

票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：

https://stackoverflow.com/questions/70700585

复制
相关文章

相似问题

领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 粤公网安备44030502008569号

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档