Gekko 最优控制空气密度方程出错

如何解决Gekko 最优控制空气密度方程出错

所以我正在使用最佳控制优化飞机飞行。飞机飞行了一定距离（path 变量），然后模拟停止。求解器试图通过最大化质量值来最小化燃料消耗 (m.Maximize(masstffinal)。

我尝试实现气压和空气密度方程。压力方程工作正常，但密度不正常。

对于气压，我使用了这个：

它变成了这个m.Equation(pressure==101325*(1-(0.0065*h)/T0)**((g*0.0289652)/(8.31446*0.0065)))

对于空气密度，我使用了这个：

我把它变成了这样：m.Equation(Ro==(pressure*0.0289652)/(8.31446*T)) 但是出于某种原因，当我尝试在密度方程处于活动状态的情况下运行程序时，模拟会一直运行，直到达到最大循环计数。我对那个密度方程做错了什么？

我的代码：

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
import math
#Gekko model
m = GEKKO(remote=False)

#Time points
nt = 11
tm =  np.linspace(0,100,nt)
m.time = tm

# Variables
Ro = m.Var(value=1.1)#air density
g = m.Const(value=9.80665)
pressure = m.Var(value=101325)#
T = m.Var(value=281)#temp
T0 = m.Const(value=288)#temp at see level
S = m.Const(value=122.6)
Cd = m.Const(value=0.1)
FuelFlow = m.Var()
D = m.Var()#drag
Thrmax = m.Const(value=200000)#macimum throttle
Thr = m.Var()
V = m.Var(value=100,lb=0,ub=240)#veLocity
#Vmin = m.Var(value=100)
nu = m.Var(value=0)#angle
nuu = nu.value
x = m.Var(value=0,lb=0)#x position
h = m.Var(value=1000)# height
mass = m.Var(value=60000)
path = m.Const(value=5000000) #intended distance length
L = m.Var()

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

m.options.MAX_ITER=10000 # iteration number

#Fixed Variable
tf = m.FV(value=1,lb=0.0001,ub=1000.0)#
tf.STATUS = 1

# Parameters
Tcontr = m.MV(value=0.2,lb=0.2,ub=1)
Tcontr.STATUS = 1
Tcontr.dcosT = 0

# Equations
m.Equation(x.dt()==tf*(V*(math.cos(nuu.value))))#
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(V.dt()==tf*((Thr-D)/mass))#
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60000)))#

m.Equation(T==T0-(0.0065*h))
m.Equation(pressure==101325*(1-(0.0065*h)/T0)**((g*0.0289652)/(8.31446*0.0065)))# equation works
m.Equation(Ro==(pressure*0.0289652)/(8.31446*T))# equation does not work

m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(x*final<=path)

# Objective Function
m.Minimize(final*(x-path)**2)
m.Maximize(mass*tf*final)
m.options.IMODE = 6
m.options.NODES = 2 # it was 3 before
m.options.MV_TYPE = 1
m.options.soLVER = 3
#m.open_folder() # to search for infeasibilities
m.solve()

tm = tm * tf.value[0]
    
fig,axs = plt.subplots(6)
fig.suptitle('Results')
axs[0].plot(tm,Tcontr,'r-',linewidth=2,label=r'$Tcontr$')
axs[0].legend(loc='best')
axs[1].plot(tm,V.value,'b-',label=r'$V$')
axs[1].legend(loc='best')
axs[2].plot(tm,x.value,'r--',label=r'$x$')
axs[2].legend(loc='best')
axs[3].plot(tm,D.value,'g-',label=r'$D$')
axs[3].legend(loc='best')
axs[4].plot(tm,mass.value,'g:',label=r'$mass$')
axs[4].legend(loc='best')
axs[5].plot(tm,T.value,'p-',label=r'$T$')
axs[5].legend(loc='best')
#axs[6].plot(tm,Ro.value,label=r'$Ro$')
#axs[6].legend(loc='best')
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()

解决方法

一些可以帮助解决找不到解决方案的问题：

为变量添加约束以帮助求解器保持在合理范围内。

# Variables
T = m.Var(value=281,lb=100)#temp
h = m.Var(value=1000,lb=0)# height
mass = m.Var(value=60000,lb=1000)
tf = m.FV(value=1,lb=0.1,ub=100.0)#

如果解决方案受到约束，则考虑进一步打开约束条件并从先前的解决方案解决或重新开始。

通过重新排列方程避免被零除。

#m.Equation(Ro==(pressure*0.0289652)/(8.31446*T))# equation does not work
m.Equation(Ro*(8.31446*T)==(pressure*0.0289652))

通常，求解器在经过几百次迭代后要么成功要么失败。设置 MAX_ITER=100000 可能会导致求解不成功的时间过长。

import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
import math
#Gekko model
m = GEKKO(remote=False)

#Time points
nt = 101
tm =  np.linspace(0,100,nt)
m.time = tm

# Variables
Ro = m.Var(value=1.1)#air density
g = m.Const(value=9.80665)
pressure = m.Var(value=101325)#
T = m.Var(value=281,lb=100)#temp
T0 = m.Const(value=288)#temp at see level
S = m.Const(value=122.6)
Cd = m.Const(value=0.1)
FuelFlow = m.Var()
D = m.Var()#drag
Thrmax = m.Const(value=200000)#maximum throttle
Thr = m.Var()
V = m.Var(value=100,lb=0,ub=240)#velocity
#Vmin = m.Var(value=100)
nu = m.Var(value=0)#angle
nuu = nu.value
x = m.Var(value=0,lb=0)#x position
h = m.Var(value=1000,lb=1000)
path = m.Const(value=5000000) #intended distance length
L = m.Var()

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

m.options.MAX_ITER=200 # iteration number

#Fixed Variable
tf = m.FV(value=1,ub=100.0)#
tf.STATUS = 1

# Parameters
Tcontr = m.MV(value=0.2,lb=0.2,ub=1)
Tcontr.STATUS = 1
Tcontr.DCOST = 0

# Equations
m.Equation(x.dt()==tf*(V*(math.cos(nuu.value))))#
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(V.dt()==tf*((Thr-D)/mass))#
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60000)))#

m.Equation(T==T0-(0.0065*h))
m.Equation(pressure==101325*(1-(0.0065*h)/T0)**((g*0.0289652)/(8.31446*0.0065)))# equation works
m.Equation(Ro*(8.31446*T)==(pressure*0.0289652))# equation does not work

m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(x*final<=path)

# Objective Function
m.Minimize(final*(x-path)**2)
m.Maximize(mass*tf*final)
m.options.IMODE = 6
m.options.NODES = 2 # it was 3 before
m.options.MV_TYPE = 1
m.options.SOLVER = 3
#m.open_folder() # to search for infeasibilities
m.solve()

tm = tm * tf.value[0]
    
fig,axs = plt.subplots(6)
fig.suptitle('Results')
axs[0].plot(tm,Tcontr,'r-',lw=2,label=r'$Tcontr$')
axs[0].legend(loc='best')
axs[1].plot(tm,V.value,'b-',label=r'$V$')
axs[1].legend(loc='best')
axs[2].plot(tm,x.value,'r--',label=r'$x$')
axs[2].legend(loc='best')
axs[3].plot(tm,D.value,'g-',label=r'$D$')
axs[3].legend(loc='best')
axs[4].plot(tm,mass.value,'g:',label=r'$mass$')
axs[4].legend(loc='best')
axs[5].plot(tm,T.value,'p-',label=r'$T$')
axs[5].legend(loc='best')
#axs[6].plot(tm,Ro.value,label=r'$Ro$')
#axs[6].legend(loc='best')
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()