使用 GEKKO 的 MPC 中的变量范围

我正在尝试使用 MPC 和 GEKKO 实现恒温器控制。

状态变量（温度）应在预先指定的温度下限和上限范围内，temp_low and temp_upper在下面的代码中。

两者的界限每天都会变化：每小时一个值。

目标函数是使用供暖的成本。价格也会随着一天的变化而变化，TOU below. T_external是在微分方程中起作用的房间外部温度。

如何实施才能优化？

这是我的尝试：

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 = [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)

u = m.MV(value=[0]*24, lb=[0.0]*24, ub=[1.]*24)
u.STATUS = 1  # allow optimizer to change

# Controlled Variable
T = m.SV(value=[60]*24, lb=temp_low, ub=temp_upper)

m.Equation(T.dt() == k*(T_e-T) + b*u)

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

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

当我运行这个时，我得到：

@error: Model Expression
 *** Error in syntax of function string: Missing operator

Position: 4                   
 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0
    ?

Gekko 需要约束作为不等式表达式，其中变量T与上层相比TH或更低TL价值观。如果你有b=1.，这会导致一个不可行的解决方案，因为加热器的功率不足以将温度保持在上限和下限内。我将值更改为b=10以获得可行的解决方案。

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=10.)
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=60)

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)

一个可能更好的解决方案是通过将限制重新定义为错误来设置软约束。你可以最小化误差以保持在限制范围内 https://apmonitor.com/do/index.php/Main/ControllerTuning。即使它不能保持在限制范围内，优化器也会尽其所能来最大程度地减少不可行性。这也让您能够权衡多个目标 https://apmonitor.com/do/index.php/Main/MultiObjectiveOptimization同时考虑舒适性和成本之间的关系。

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=10.)
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=60)

# Soft constraints
eH = m.CV(value=0)
eL = m.CV(value=0)

eH.SPHI=0; eH.WSPHI=100; eH.WSPLO=0  ; eH.STATUS = 1
eL.SPLO=0; eL.WSPHI=0  ; eL.WSPLO=100; eL.STATUS = 1

m.Equations([eH==T-TH,eL==T-TL])

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()