python - GEKKO 最佳控制。如何强制值高于某些阈值?
问题描述
所以我在模拟飞机飞行。飞机飞行一定距离(path
变量),然后模拟停止。求解器试图(m.Maximize(mass*tf*final)
通过最大化质量值来最小化燃料消耗。求解器控制油门踏板位置 ( Tcontr
)。我需要它来尝试保持速度值 V 高于某个阈值(例如 >=100)。我需要添加什么样的方程式才能发生这种情况?
我的代码:
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.Const(value=1.1)#
g = m.Const(value=9.80665)
pressure = m.Var(value=101325)#
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 = m.Const(value=200000)
Thr = m.Var()
V = m.Var(value=100,lb=0,ub=240)
Vmin = m.Var(value=100)
nu = m.Var(value=0)
nuu = nu.value
x = m.Var(value=0,lb=0)
h = m.Var(value=1000)
mass = m.Var(value=60000)
path = m.Const(value=5000000) #intended distance travelled
L = m.Var()
p = np.zeros(nt)
p[-1] = 1.0
final = m.Param(value=p)
m.options.MAX_ITER=50000 # 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(D==0.5*Ro*(V**2)*Cd*S)
m.Equation(FuelFlow==0.75882*(1+(V/2938.5)))
m.Equation(x*final<=path)
m.Equation(T==T0-(0.0065*h))
# 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]
print('Final Time: ' + str((tf.value[-1])*100))
print('Final Speed: ' + str(V.value[-1]))
print('Final X: ' + str(x.value[-1]))
print('Final Mass: ' + str(mass.value[-1]))
fig, axs = plt.subplots(5)
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,mass.value,'g:',LineWidth=2,label=r'$mass$')
axs[4].legend(loc='best')
plt.xlabel('Time')
#plt.ylabel('Value')
plt.show()
解决方案
在 Gekko 中有两种实现方式。
硬约束
第一种方法是将 Velocity 的下限设置为 100:
V = m.Var(value=100,lb=100,ub=240)
这是一个硬约束——除非满足这一点,否则求解器不会返回成功的解决方案。然而,这可能导致不可行的解决方案。
软约束
如果需要,但不需要 Velocity 保持在 100 以上,或者存在不可行解决方案的问题,那么软约束可能是更好的选择。设置以下以实现 100 到 240 之间的死区:
V = m.CV(value=100,lb=0,ub=240)
V.SPHI = 240; V.SPLO = 100
V.WSPHI = 1; V.WSPLO = 10 # penalty for violating bound
V.STATUS = 1; V.TR_INIT = 0
机器学习和动态优化课程中提供了有关此类变量的更多信息,并详细描述了 l1-norm 目标。
推荐阅读
- angular-cli - 为 angular-cli v6 全局关闭进度显示
- curl - 如何使用 cURL 获取 FULL cookie?
- html - 我想要每个节点 d3 的特定路径颜色
- java - 矩阵格式的控制台输出
- python - 使用 python sql 连接器执行简单的存储过程
- python - 替换字符串中的每个第二个字符
- ubuntu - XAMPP 另一个 Web 服务器已经在运行。尝试
- tensorflow - 无法将 Tensorflow 标量摘要写入事件日志
- wpf - WPF DataGrid:获取单元格行索引的有效方法
- javascript - 在 Node.js 中的 Array 中使对象数组成为对象的值