使用 GEKKO 的 MPC 中的可变边界

Question

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

状态变量（温度）应该在预先指定的温度值的上限和上限内，并在下面的代码中。temp_lowtemp_upper

两者绑定在一天中都有所不同：每小时一个值。

目标函数是使用加热的成本。价格也随天变化，如下所示。 是房间的外部温度，在微分方程中起作用。TOUT_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

    ?

Answer 1

Gekko 需要将 constaint 作为不等式表达式，其中变量与上限或下限值进行比较。如果你有，它会导致一个不可行的解决方案，因为加热器的功率不足以将温度保持在上限和下限。我将值更改为 以获得可行的解决方案。TTHTLb=1.b=10from gekko import GEKKOimport numpy as npm = GEKKO(remote=False)m.time = np.linspace(0,23,24)#initialize variablesT_external = [50.,50.,50.,50.,45.,45.,45.,60.,60.,63.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 64.,45.,45.,50.,52.,53.,53.,54.,54.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 53.,52.,51.,50.,45.]temp_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 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.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 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.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;200.,200.,200.,200.,200.,200.,200.,200.,200.,200.,0.05,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;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&nbsp; # allow optimizer to change# Controlled VariableT = 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 = 6m.solve(disp=True,debug=True)一个可能更好的解决方案是通过将限制重新定义为错误来设置软约束。您可以将误差降至最低以保持在限制范围内。即使它不能保持在限制范围内，优化程序也会尽最大努力将不可行性降至最低。这还允许您同时权衡多个目标，例如在舒适度和成本之间from gekko import GEKKOimport numpy as npm = GEKKO(remote=False)m.time = np.linspace(0,23,24)#initialize variablesT_external = [50.,50.,50.,50.,45.,45.,45.,60.,60.,63.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 64.,45.,45.,50.,52.,53.,53.,54.,54.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 53.,52.,51.,50.,45.]temp_low = [55.,55.,55.,55.,55.,55.,55.,68.,68.,68.,68.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 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.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; 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.,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;200.,200.,200.,200.,200.,200.,200.,200.,200.,200.,0.05,\&nbsp; &nbsp; &nbsp; &nbsp; &nbsp;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&nbsp; # allow optimizer to change# Controlled VariableT = m.SV(value=60)# Soft constraintseH = m.CV(value=0)eL = m.CV(value=0)eH.SPHI=0; eH.WSPHI=100; eH.WSPLO=0&nbsp; ; eH.STATUS = 1eL.SPLO=0; eL.WSPHI=0&nbsp; ; eL.WSPLO=100; eL.STATUS = 1m.Equations([eH==T-TH,eL==T-TL])m.Equation(T.dt() == k*(T_e-T) + b*u)m.Minimize(TOU*u)m.options.IMODE = 6m.solve(disp=True,debug=True)import matplotlib.pyplot as pltplt.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()