python - 壁虎。 X值不超过一定点

Question

我需要解决一维飞机飞行最优控制问题。我有一架 1000 米高的飞机。我需要它沿 x 轴向前行驶一定距离 (x)，同时最大限度地减少油耗。当它达到那个距离 x 时，我需要程序停止。这个函数控制它:m.Equation(x*final<=1500) .

由于某些原因，在模拟期间我的 x 值不想高于 1310。

我怎样才能解决这个“障碍”？

我的gekko脚本:

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,1,nt)
m.time = tm

# Variables
Ro = m.Const(value=1.1)
g = m.Const(value=9.80665)
T = m.Var()
T0 = m.Const(value=273)
S = m.Const(value=122.6)
Cd = m.Const(value=0.1)
FuelFlow = m.Var()
D = m.Var()
Thrmax = value=200000
Thr = m.Var()
V = m.Var()
nu = m.Var(value=0)
nuu = nu.value
x = m.Var(value=0)
h = m.Var(value=1000)
mass = m.Var(value=60000)
    
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

m.options.MAX_ITER=40000 # iteration max number

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

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

# 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/60)))#
m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(x*final<=1500)
m.Equation(T==T0-h)

# Objective Function
m.Obj(-mass*tf*final)#
m.options.IMODE = 6
m.options.NODES = 3
m.options.MV_TYPE = 1
m.options.SOLVER = 3
#m.open_folder() # to search for infeasibilities
m.solve()


tm = tm * tf.value[0]

print('Final Time: ' + str(tf.value[-1]))
print('Final Speed: ' + str(V.value[-1]))
print('Final X: ' + str(x.value[-1]))

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(tm,Tcontr,'r-',LineWidth=2,label=r'$Tcontr$')
#plt.plot(m.time,x.value,'r--',LineWidth=2,label=r'$x$')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(tm,x.value,'r--',LineWidth=2,label=r'$x$')
#plt.plot(tm,mass.value,'g:',LineWidth=2,label=r'$mass$')
#plt.plot(tm,D.value,'g:',LineWidth=2,label=r'$D$')
#plt.plot(tm,V.value,'b-',LineWidth=2,label=r'$V$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()

Answer 1

Tcontr 的下界正在阻止 x值超过 1310。将较低的值设置为 0.1 会提高 x 的最终值。到 2469.76 如果约束 m.Equation(x*final<=1500)已删除。

Tcontr = m.MV(value=0.2,lb=0.1,ub=1)

当解决方案不是最优的或存在不可行的解决方案时，约束通常是罪魁祸首。检测有问题的约束的一种方法是创建绘图以验证解决方案未被人为约束。

制定约束的另一种方法是使用硬约束的组合 m.Equation((x-1500)*final==0)和一个软约束来指导解决方案为 m.Minimize(final*(x-1500)**2) .正确设置硬约束很重要。约束条件，例如 m.Equation(x*final==1500)意味着 x*0==1500 不可行当 final 不等于 1 时。不等式版本也更好地构成为 m.Equation((x-1500)*final<=0)但原始形式也有效，因为 x*0<1500 .

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,1,nt)
m.time = tm

# Variables
Ro = m.Const(value=1.1)
g = m.Const(value=9.80665)
T = m.Var()
T0 = m.Const(value=273)
S = m.Const(value=122.6)
Cd = m.Const(value=0.1)
FuelFlow = m.Var()
D = m.Var()
Thrmax = value=200000
Thr = m.Var()
V = m.Var()
nu = m.Var(value=0)
nuu = 0
x = m.Var(value=0)
h = m.Var(value=1000)
mass = m.Var(value=60000)
    
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)

m.options.MAX_ITER=1000 # iteration max number

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

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

# Equations
m.Equation(x.dt()==tf*(V*(math.cos(nuu))))#
m.Equation(Thr==Tcontr*Thrmax)
m.Equation(V.dt()==tf*((Thr-D)/mass))#
m.Equation(mass.dt()==tf*(-Thr*(FuelFlow/60)))#
m.Equation(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation((x-1500)*final==0)
m.Equation(T==T0-h)

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


tm = tm * tf.value[0]

print('Final Time: ' + str(tf.value[-1]))
print('Final Speed: ' + str(V.value[-1]))
print('Final X: ' + str(x.value[-1]))

plt.figure(1)
plt.subplot(2,1,1)
plt.plot(tm,Tcontr,'r-',lw=2,label=r'$Tcontr$')
#plt.plot(m.time,x.value,'r--',lw=2,label=r'$x$')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(tm,x.value,'r--',lw=2,label=r'$x$')
#plt.plot(tm,mass.value,'g:',lw=2,label=r'$mass$')
#plt.plot(tm,D.value,'g:',lw=2,label=r'$D$')
#plt.plot(tm,V.value,'b-',lw=2,label=r'$V$')
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Value')
plt.show()