如何使用 Gekko 求解可重复使用再入飞行器的最优控制
How to use Gekko to solve for optimal control for a reusable reentry vehicle
我正在寻找最佳控制（aoa 和坡度角）以最大化使用 Gekko 的穿梭型再入飞行器的横向范围。下面是我目前的代码，我得到一个“未找到解决方案”和“退出：超过最大迭代次数”。模拟假设点质量具有非旋转地球框架。 EOMS 是 6 个耦合的非线性 ODE。我尝试过使用不同的求解器、implementing/removing 状态和控制约束、增加最大迭代次数等。我对 Gekko 中问题的设置和实现没有信心，希望得到一些关于我可以尝试的反馈下一个。我尝试遵循 APMonitor's Example 11. Optimal Control with Integral Objective, Inverted Pendulum Optimal Control, and Example 13. Optimal Control: Minimize Final Time 中的设置和布局。我正在寻找的解决方案如下。
非常感谢任何帮助！
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 和坡度角的解决方案应如下所示：Altitude Velocity AngleofAttack Bank
我通过减少最终时间（成功解的最大值=0.04）并重新排列方程以避免可能的被零除得到了一个成功的解：
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)))
从一个小时间段开始并观察状态有助于对模型进行故障排除。
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)。可能还需要额外的约束以将三角函数保持在 -2pi 到 2pi 区域中。有关初始化具有挑战性的问题的更多信息：
- Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, 计算机与化学工程, 2015,卷。 78，第 39-50 页，DOI：10.1016/j.compchemeng.2015.04.016.
成功！！感谢@John Hedengren 的反馈和指点。结合我在这里学到的知识并恢复状态的界限后，我能够成功地使用 Gekko 生成与以前的解决方案一致的答案。对 Gekko 和这里的潜力感到非常兴奋！
我正在寻找最佳控制（aoa 和坡度角）以最大化使用 Gekko 的穿梭型再入飞行器的横向范围。下面是我目前的代码，我得到一个“未找到解决方案”和“退出：超过最大迭代次数”。模拟假设点质量具有非旋转地球框架。 EOMS 是 6 个耦合的非线性 ODE。我尝试过使用不同的求解器、implementing/removing 状态和控制约束、增加最大迭代次数等。我对 Gekko 中问题的设置和实现没有信心，希望得到一些关于我可以尝试的反馈下一个。我尝试遵循 APMonitor's Example 11. Optimal Control with Integral Objective, Inverted Pendulum Optimal Control, and Example 13. Optimal Control: Minimize Final Time 中的设置和布局。我正在寻找的解决方案如下。
非常感谢任何帮助！
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 和坡度角的解决方案应如下所示：Altitude Velocity AngleofAttack Bank
我通过减少最终时间（成功解的最大值=0.04）并重新排列方程以避免可能的被零除得到了一个成功的解：
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)))
从一个小时间段开始并观察状态有助于对模型进行故障排除。
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)。可能还需要额外的约束以将三角函数保持在 -2pi 到 2pi 区域中。有关初始化具有挑战性的问题的更多信息：
- Safdarnejad, S.M., Hedengren, J.D., Lewis, N.R., Haseltine, E., Initialization Strategies for Optimization of Dynamic Systems, 计算机与化学工程, 2015,卷。 78，第 39-50 页，DOI：10.1016/j.compchemeng.2015.04.016.
成功！！感谢@John Hedengren 的反馈和指点。结合我在这里学到的知识并恢复状态的界限后，我能够成功地使用 Gekko 生成与以前的解决方案一致的答案。对 Gekko 和这里的潜力感到非常兴奋！