Gekko 最优控制。我无法将一个值保持在其范围内
Gekko optimal control. I can't keep one value within its bounds
我正在尝试模拟 3D 飞机飞行。我对伽马值（飞行路径角度）有疑问。它超出了范围，然后模拟停止。伽马值是通过以下等式计算的：
我把它变成了这个：m.Equation(gamma.dt()==tf*((L*m.cos(Mu)-mass*g*m.cos(gamma))/mass*V))
模拟的目标是让飞机达到特定的 X 和 Y 值（m.Minimize(w*final*(x-pathx)**2) 和 m.Minimize(w*final*(pathy-y)**2)），同时最大限度地减少燃料消耗 m.Maximize(0.2*mass*tf*final)。
求解器通过控制升力系数 Cl 来控制 gamma 值，这会影响升力值 L，进而影响 gamma 值。计算提升 L 值的方程如下所示：m.Equation(L==0.5*Ro*(V**2)*(Cl)*S)。但最终求解器无法控制 gamma 值，直到飞机到达目的地。
有什么问题？
我的代码：
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,lb=100)#temperature
T0 = m.Const(value=288)#temperature at see level
S = m.Const(value=122.6)
Cd = m.Var(value=0.025)#drag coef 0.06 works
#Cl = m.Const(value=0.3)#lift couef
FuelFlow = m.Var()
D = m.Var(value=10000,lb=0)#drag
Thrmax = m.Const(value=200000)#maximum throttle
Thr = m.Var()
V = m.Var(value=200,lb=45,ub=240)#velocity
gamma = m.Var(value=0,lb=-0.6,ub=1.2)# Flight-path angle
#gammaa = gamma.value
Xi = m.Var(value=0, lb=-2, ub=2.0)# Heading angle
x = m.Var(value=0,lb=-1000,ub=1015000)#x position
y = m.Var(value=0,lb=-1000,ub=1011000)#y position
h = m.Var(value=1000,lb=-20000,ub=50000)# height
targetH = m.Param(value=10000) #target flight altitude
mass = m.Var(value=60000,lb=10000)
pathx = m.Const(value=1000000) #intended distance in x direction
pathy = m.Const(value=1000000) #intended distance in y direction
L = m.Var(value=250000)#lift
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=100.0)#
tf.STATUS = 1
# Controlled parameters
Tcontr = m.MV(value=0.6,lb=0.0,ub=1)# solver controls throttle pedal position
Tcontr.STATUS = 1
Tcontr.DCOST = 0
Mu = m.MV(value=0,lb=-1,ub=1)# solver controls bank angle
Mu.STATUS = 1
Mu.DCOST = 1e-3
Cl = m.MV(value=0.1,lb=-0.3,ub=0.9)# solver controls lift couef
Cl.STATUS = 1
Cl.DCOST = 1e-3
# Equations
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60000)))#
m.Equation(V.dt()==tf*((Thr-D-mass*g*m.cos(gamma))/mass)) #
m.Equation(x.dt()==tf*(V*(m.cos(gamma))*(m.cos(Xi))))#
m.Equation(x*final<=pathx)
#pressure and density part
m.Equation(T==T0-(0.0065*h))
m.Equation(pressure==101325*(1-(0.0065*h)/T0)**((g*0.0289652)/(8.31446*0.0065)))#
m.Equation(Ro*(8.31446*T)==(pressure*0.0289652))
#2D addition part
m.Equation(L==0.5*Ro*(V**2)*(Cl)*S)# Problem here or no problem, idk
m.Equation(mass*m.cos(gamma)*V*Xi.dt()==tf*((L*m.sin(Mu)))) #
m.Equation(y.dt()==tf*(V*(m.cos(gamma))*(m.sin(Xi))))#
#m.Equation((y-pathy)*final==0)
m.Equation(y*final<=pathy)
#3D addition part
m.Equation(h.dt()==tf*(V*(m.sin(gamma))))#
m.Equation(h*final<=targetH)
m.Equation(gamma.dt()==tf*((L*m.cos(Mu)-mass*g*m.cos(gamma))/mass*V))#
#Cd equation
m.Equation(Cd==((Cl)**2)/10)
# Objective Function
w = 1e4
m.Minimize(w*final*(x-pathx)**2) #1D part (x)
m.Minimize(w*final*(pathy-y)**2) #2D part (y)
m.Maximize(0.2*mass*tf*final) #objective function
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(11)
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-',LineWidth=2,label=r'$V$')
axs[1].legend(loc='best')
axs[2].plot(tm,x.value,'r--',LineWidth=2,label=r'$x$')
axs[2].legend(loc='best')
axs[3].plot(tm,D.value,'g-',LineWidth=2,label=r'$D$')
axs[3].legend(loc='best')
axs[4].plot(tm,L.value,'p-',LineWidth=2,label=r'$L$')
axs[4].legend(loc='best')
axs[5].plot(tm,y.value,'p-',LineWidth=2,label=r'$y$')
axs[5].legend(loc='best')
axs[6].plot(tm,Ro.value,'p-',LineWidth=2,label=r'$Ro$')
axs[6].legend(loc='best')
axs[7].plot(tm,h.value,'p-',LineWidth=2,label=r'$h$')
axs[7].legend(loc='best')
axs[8].plot(tm,gamma.value,'p-',LineWidth=2,label=r'$gamma$')
axs[8].legend(loc='best')
axs[9].plot(tm,Cl.value,'p-',LineWidth=2,label=r'$Cl$')
axs[9].legend(loc='best')
axs[10].plot(tm,Cd.value,'p-',LineWidth=2,label=r'$Cd$')
axs[10].legend(loc='best')
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()
重要更新
我希望求解器能够将 gamma 值保持在其边界内，直到获得必要的 x 和 y 值。求解器有问题。
图表 A.
在这种情况下，模拟停止了，因为它无法阻止 gamma 超出其下限。 gamma 通常在 L (lift) * cos(mu) 大于 mass*g*cos(gamma) 时变大，反之则变小： mass*g*cos(gamma) 大于 L (lift) * cos(mu)。那么问题就变成了：为什么突然 mass*g*cos(gamma) 变得比 L (lift) * cos(mu) 大这么多？ g 是不变的。在模拟的最后时刻，质量没有太大变化。 L（提升）并没有变得特别小。在这部分模拟过程中，cos(mu) 通常等于 1。
图表 B.
在这种情况下，模拟再次停止，因为它无法将 gamma 值保持在其上限内。可以看出，在模拟的最后时刻，Cl 值正在上升，有必要防止伽玛值超过其下限，但在那之后，gamma 值飙升，并且由于某种原因, 求解器没有降低 Cl 值，这迫使 gamma 值超过其上限，从而迫使模拟在达到目标之前停止。
在这种情况下，问题是：为什么求解器不降低 Cl 值以阻止 gamma 超过其上限？
我不确定到底是什么问题。看起来优化器正在整个范围内控制伽玛，并且它始终保持在 -0.6 到 1.2 的范围内。您能否提供有关问题所在的更多信息？
max(gamma.value)
>>>> 1.2
min(gamma.value)
>>>> -0.6
这是我通过 运行 您的代码生成的解决方案：
请注意，如果您使用更高的时间分辨率解决问题，您会得到一个完全不同的解决方案：
我正在尝试模拟 3D 飞机飞行。我对伽马值（飞行路径角度）有疑问。它超出了范围，然后模拟停止。伽马值是通过以下等式计算的：
我把它变成了这个：m.Equation(gamma.dt()==tf*((L*m.cos(Mu)-mass*g*m.cos(gamma))/mass*V))
模拟的目标是让飞机达到特定的 X 和 Y 值（m.Minimize(w*final*(x-pathx)**2) 和 m.Minimize(w*final*(pathy-y)**2)），同时最大限度地减少燃料消耗 m.Maximize(0.2*mass*tf*final)。
求解器通过控制升力系数 Cl 来控制 gamma 值，这会影响升力值 L，进而影响 gamma 值。计算提升 L 值的方程如下所示：m.Equation(L==0.5*Ro*(V**2)*(Cl)*S)。但最终求解器无法控制 gamma 值，直到飞机到达目的地。
有什么问题？
我的代码：
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,lb=100)#temperature
T0 = m.Const(value=288)#temperature at see level
S = m.Const(value=122.6)
Cd = m.Var(value=0.025)#drag coef 0.06 works
#Cl = m.Const(value=0.3)#lift couef
FuelFlow = m.Var()
D = m.Var(value=10000,lb=0)#drag
Thrmax = m.Const(value=200000)#maximum throttle
Thr = m.Var()
V = m.Var(value=200,lb=45,ub=240)#velocity
gamma = m.Var(value=0,lb=-0.6,ub=1.2)# Flight-path angle
#gammaa = gamma.value
Xi = m.Var(value=0, lb=-2, ub=2.0)# Heading angle
x = m.Var(value=0,lb=-1000,ub=1015000)#x position
y = m.Var(value=0,lb=-1000,ub=1011000)#y position
h = m.Var(value=1000,lb=-20000,ub=50000)# height
targetH = m.Param(value=10000) #target flight altitude
mass = m.Var(value=60000,lb=10000)
pathx = m.Const(value=1000000) #intended distance in x direction
pathy = m.Const(value=1000000) #intended distance in y direction
L = m.Var(value=250000)#lift
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=100.0)#
tf.STATUS = 1
# Controlled parameters
Tcontr = m.MV(value=0.6,lb=0.0,ub=1)# solver controls throttle pedal position
Tcontr.STATUS = 1
Tcontr.DCOST = 0
Mu = m.MV(value=0,lb=-1,ub=1)# solver controls bank angle
Mu.STATUS = 1
Mu.DCOST = 1e-3
Cl = m.MV(value=0.1,lb=-0.3,ub=0.9)# solver controls lift couef
Cl.STATUS = 1
Cl.DCOST = 1e-3
# Equations
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60000)))#
m.Equation(V.dt()==tf*((Thr-D-mass*g*m.cos(gamma))/mass)) #
m.Equation(x.dt()==tf*(V*(m.cos(gamma))*(m.cos(Xi))))#
m.Equation(x*final<=pathx)
#pressure and density part
m.Equation(T==T0-(0.0065*h))
m.Equation(pressure==101325*(1-(0.0065*h)/T0)**((g*0.0289652)/(8.31446*0.0065)))#
m.Equation(Ro*(8.31446*T)==(pressure*0.0289652))
#2D addition part
m.Equation(L==0.5*Ro*(V**2)*(Cl)*S)# Problem here or no problem, idk
m.Equation(mass*m.cos(gamma)*V*Xi.dt()==tf*((L*m.sin(Mu)))) #
m.Equation(y.dt()==tf*(V*(m.cos(gamma))*(m.sin(Xi))))#
#m.Equation((y-pathy)*final==0)
m.Equation(y*final<=pathy)
#3D addition part
m.Equation(h.dt()==tf*(V*(m.sin(gamma))))#
m.Equation(h*final<=targetH)
m.Equation(gamma.dt()==tf*((L*m.cos(Mu)-mass*g*m.cos(gamma))/mass*V))#
#Cd equation
m.Equation(Cd==((Cl)**2)/10)
# Objective Function
w = 1e4
m.Minimize(w*final*(x-pathx)**2) #1D part (x)
m.Minimize(w*final*(pathy-y)**2) #2D part (y)
m.Maximize(0.2*mass*tf*final) #objective function
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(11)
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-',LineWidth=2,label=r'$V$')
axs[1].legend(loc='best')
axs[2].plot(tm,x.value,'r--',LineWidth=2,label=r'$x$')
axs[2].legend(loc='best')
axs[3].plot(tm,D.value,'g-',LineWidth=2,label=r'$D$')
axs[3].legend(loc='best')
axs[4].plot(tm,L.value,'p-',LineWidth=2,label=r'$L$')
axs[4].legend(loc='best')
axs[5].plot(tm,y.value,'p-',LineWidth=2,label=r'$y$')
axs[5].legend(loc='best')
axs[6].plot(tm,Ro.value,'p-',LineWidth=2,label=r'$Ro$')
axs[6].legend(loc='best')
axs[7].plot(tm,h.value,'p-',LineWidth=2,label=r'$h$')
axs[7].legend(loc='best')
axs[8].plot(tm,gamma.value,'p-',LineWidth=2,label=r'$gamma$')
axs[8].legend(loc='best')
axs[9].plot(tm,Cl.value,'p-',LineWidth=2,label=r'$Cl$')
axs[9].legend(loc='best')
axs[10].plot(tm,Cd.value,'p-',LineWidth=2,label=r'$Cd$')
axs[10].legend(loc='best')
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()
重要更新
我希望求解器能够将 gamma 值保持在其边界内，直到获得必要的 x 和 y 值。求解器有问题。
图表 A.
gamma 超出其下限。 gamma 通常在 L (lift) * cos(mu) 大于 mass*g*cos(gamma) 时变大，反之则变小： mass*g*cos(gamma) 大于 L (lift) * cos(mu)。那么问题就变成了：为什么突然 mass*g*cos(gamma) 变得比 L (lift) * cos(mu) 大这么多？ g 是不变的。在模拟的最后时刻，质量没有太大变化。 L（提升）并没有变得特别小。在这部分模拟过程中，cos(mu) 通常等于 1。
图表 B.
gamma 值保持在其上限内。可以看出，在模拟的最后时刻，Cl 值正在上升，有必要防止伽玛值超过其下限，但在那之后，gamma 值飙升，并且由于某种原因, 求解器没有降低 Cl 值，这迫使 gamma 值超过其上限，从而迫使模拟在达到目标之前停止。
在这种情况下，问题是：为什么求解器不降低 Cl 值以阻止 gamma 超过其上限？
我不确定到底是什么问题。看起来优化器正在整个范围内控制伽玛，并且它始终保持在 -0.6 到 1.2 的范围内。您能否提供有关问题所在的更多信息？
max(gamma.value)
>>>> 1.2
min(gamma.value)
>>>> -0.6
这是我通过 运行 您的代码生成的解决方案：
请注意，如果您使用更高的时间分辨率解决问题，您会得到一个完全不同的解决方案：