GEKKO RTO 与 MPC 模式

这是一个由此衍生的问题one https://stackoverflow.com/questions/60761440/variable-bounds-in-mpc-with-gekko。在发布我的问题后，我找到了一个解决方案（更像是强制优化器优化的补丁）。有件事让我困惑。约翰·赫登格伦 https://stackoverflow.com/users/2366941/john-hedengren正确地指出ab=1.0在 ODE 中会导致不可行的解决方案IMODE=6。然而，在我的零散工作中IMODE=3我确实得到了解决方案。

我试图了解这里发生的事情阅读GEKKO's https://gekko.readthedocs.io/en/latest/imode.html的文档IMODE=3 and 6但我不清楚

IMODE=3

RTO 实时优化 (RTO) 是一种稳态模式，允许 决策变量（STATUS=1 的 FV 或 MV 类型）或其他 变量超过方程的数量。目标函数 指导附加变量的选择以选择最优的 可行的解决方案。 RTO 是 Gekko 的默认模式，如果 未指定 m.options.IMODE。

IMODE=6

MPC 模型预测控制 (MPC) 是通过 IMODE=6 实现的 同时解决方案或使用 IMODE=9 作为连续拍摄方法。

为什么b=1。在一种模式下有效但在另一种模式下无效？

这是我的零碎工作IMODE=3 and b=1.0:

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

m = GEKKO(remote=False)
m.time = np.linspace(0,23,24)

#initialize variables
T_e = [50.,50.,50.,50.,45.,45.,45.,60.,60.,63.,\
       64.,45.,45.,50.,52.,53.,53.,54.,54.,53.,52.,51.,50.,45.]
temp_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,55.,55.,68.,\
            68.,68.,68.,55.,55.,55.,55.,55.,55.,55.]
temp_upper = [75.,75.,75.,75.,75.,75.,75.,70.,70.,70.,70.,75.,75.,\
              70.,70.,70.,70.,75.,75.,75.,75.,75.,75.,75.]
TOU = [0.05,0.05,0.05,0.05,0.05,0.05,0.05,200.,200.,\
       200.,200.,200.,200.,200.,200.,200.,200.,200.,\
       200.,200.,200.,0.05,0.05,0.05]

b = m.Param(value=1.)
k = m.Param(value=0.05)

u = [m.MV(0.,lb=0.,ub=1.) for i in range(24)]

# Controlled Variable
T = [m.SV(60.,lb=temp_low[i],ub=temp_upper[i]) for i in range(24)]

for i in range(24):
    u[i].STATUS = 1

for i in range(23):
    m.Equation( T[i+1]-T[i]-k*(T_e[i]-T[i])-b*u[i]==0.0 )

m.Obj(np.dot(TOU,u))

m.options.IMODE = 3
m.solve(debug=True)
myu =[u[0:][i][0] for i in range(24)]
print myu
myt =[T[0:][i][0] for i in range(24)]
plt.plot(myt)
plt.plot(temp_low)
plt.plot(temp_upper)
plt.show()
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(myu,color='b')
ax2.plot(TOU,color='k')
plt.show()

Results:

不可行与不可行的区别IMODE=6以及可行的IMODE=3那是IMODE=3case 允许优化器调整温度初始条件。优化器认识到初始条件可以更改，因此将其修改为 75，以保持可行并最大限度地减少未来的能耗。

from gekko import GEKKO
import numpy as np

m = GEKKO(remote=False)
m.time = np.linspace(0,23,24)

#initialize variables
T_external = [50.,50.,50.,50.,45.,45.,45.,60.,60.,63.,\
              64.,45.,45.,50.,52.,53.,53.,54.,54.,\
              53.,52.,51.,50.,45.]
temp_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,\
            55.,55.,68.,68.,68.,68.,55.,55.,55.,55.,55.,55.,55.]
temp_upper = [75.,75.,75.,75.,75.,75.,75.,70.,70.,70.,70.,75.,\
              75.,70.,70.,70.,70.,75.,75.,75.,75.,75.,75.,75.]
TOU_v = [0.05,0.05,0.05,0.05,0.05,0.05,0.05,200.,200.,200.,200.,\
         200.,200.,200.,200.,200.,200.,200.,200.,200.,200.,0.05,\
         0.05,0.05]

b = m.Param(value=1.)
k = m.Param(value=0.05)
T_e = m.Param(value=T_external)
TL = m.Param(value=temp_low)
TH = m.Param(value=temp_upper)
TOU = m.Param(value=TOU_v)

u = m.MV(lb=0, ub=1)
u.STATUS = 1  # allow optimizer to change

# Controlled Variable
T = m.SV(value=75)

m.Equations([T>=TL,T<=TH])
m.Equation(T.dt() == k*(T_e-T) + b*u)

m.Minimize(TOU*u)

m.options.IMODE = 6
m.solve(disp=True,debug=True)

import matplotlib.pyplot as plt
plt.subplot(2,1,1)
plt.plot(m.time,temp_low,'k--')
plt.plot(m.time,temp_upper,'k--')
plt.plot(m.time,T.value,'r-')
plt.ylabel('Temperature')
plt.subplot(2,1,2)
plt.step(m.time,u.value,'b:')
plt.ylabel('Heater')
plt.xlabel('Time (hr)')
plt.show()

如果你再去一天（48 小时），你可能会发现这个问题最终将不可行，因为较小的加热器b=1无法满足温度较低的限制。

使用的优点之一IMODE=6是你可以写出微分方程，而不是自己进行离散化。和IMODE=3，您可以使用欧拉方法来求解微分方程。默认离散化为IMODE>=4 is NODES=2，相当于欧拉有限差分法。环境NODES=3-6提高准确度有限元正交配置 https://apmonitor.com/do/index.php/Main/OrthogonalCollocation.