GEKKO RTO与MPC模式

Question

这是从这个一派生出来的一个问题。在发布了我的问题后，我找到了一个解决方案(更像是一个补丁来强制优化器进行优化)。有件事让我很困惑。约翰·赫登格伦正确地指出，ODE中的b=1.0会导致IMODE=6不可行的解决方案。然而，在我与IMODE=3的杂乱无章的工作中，我确实得到了一个解决方案。

我试图了解这里发生了什么，阅读盖科氏文档中的IMODE=3和6，但对我来说还不清楚

IMODE=3

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

IMODE=6

MPC模型预测控制(MPC)是以IMODE=6作为同步解决方案，IMODE=9作为顺序射击方法实现的。

为什么b=1.在一种模式下工作，而不是在另一种模式下工作？

这是我与IMODE=3和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()

结果：

​

​

​

​

​

​

Answer 1

不可行IMODE=6与可行IMODE=3的区别在于，IMODE=3情况允许优化器调整温度初始条件。优化器认识到初始条件可以改变，因此它将其修改为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的默认离散化是NODES=2，相当于欧拉有限差分法。设置NODES=3-6可以提高有限元上的正交配置的准确性。