Gekko 最优控制我不能在其范围内保持一个值

问题描述

我正在尝试模拟 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)。

求解器通过控制升力系数 gamma 来控制 Cl 值，该系数影响升力值 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,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-',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,L.value,'p-',label=r'$L$')
axs[4].legend(loc='best')
axs[5].plot(tm,y.value,label=r'$y$')
axs[5].legend(loc='best')
axs[6].plot(tm,Ro.value,label=r'$Ro$')
axs[6].legend(loc='best')
axs[7].plot(tm,h.value,label=r'$h$')
axs[7].legend(loc='best')
axs[8].plot(tm,gamma.value,label=r'$gamma$')
axs[8].legend(loc='best')
axs[9].plot(tm,Cl.value,label=r'$Cl$')
axs[9].legend(loc='best')
axs[10].plot(tm,Cd.value,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 超过其下限，但在那之后，gamma 值飙升，并且由于某种原因，求解器没有降低 Cl 值，这迫使 gamma 值超过其上限，这迫使模拟在达到目标之前停止。

在这种情况下，问题是：为什么求解器不降低 Cl 值以阻止 gamma 超出其上限？

解决方法

我不确定到底是什么问题。看起来优化器在整个范围内控制伽马，并且它始终保持在 -0.6 到 1.2 的范围内。您能否提供更多有关问题的信息？

max(gamma.value)
>>>> 1.2

min(gamma.value)
>>>> -0.6

这是我通过运行您的代码生成的解决方案：

请注意，如果您使用更高的时间分辨率解决问题，则会得到完全不同的解决方案：

gekko python