Границы переменных в MPC с GEKKO

Я пытаюсь реализовать управление термостатом с помощью MPC и GEKKO.

Переменная состояния (температура) должна находиться в пределах нижнего и верхнего предустановленных значений температуры, temp_low и temp_upper в приведенном ниже коде.

Обе границы меняются в течение дня: одно значение в час.

Целевая функция - это затраты на отопление. Цена также меняется в течение дня, TOU ниже. 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
    ?

person Arraval    schedule 19.03.2020    source источник


Ответы (1)


Гекко нужны константы как выражения неравенства, в которых переменная 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)

Потенциально лучшим решением является установка мягких ограничений путем переопределения пределов как ошибки. Вы можете минимизировать ошибку, чтобы оставаться в установленных пределах. Даже если он не может оставаться в установленных пределах, оптимизатор сделает все возможное, чтобы минимизировать неосуществимость. Это также позволяет вам найти компромисс между несколькими целями одновременно, например, между комфортом и Стоимость.

Не выходите за пределы допустимой температуры

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()
person John Hedengren    schedule 19.03.2020
comment
Большое спасибо за помощь, это действительно впечатляет. Есть кое-что, что меня заинтриговало. Вчера я продолжал пытаться заставить его работать, и я, наконец, сделал, но с IMODE=3 и b=1.0. Чем отличается то, что в одном режиме он работает, а в другом приводит к недопустимому решению. Я разместил вопрос здесь: stackoverflow.com/questions/60772432 / gekko-rto-vs-mpc-Mode. С уважением. - person Arraval; 20.03.2020
comment
Я буду рад добавить более подробный ответ на этот пост. IMODE=3 устанавливает все производные значения равными нулю, но ваше решение с этим другим режимом также позволяет начальному условию изменяться на более высокое значение, поэтому нагреватель не должен быть таким мощным, чтобы оставаться выполнимым. - person John Hedengren; 20.03.2020
comment
Как мне написать такие мягкие ограничения в целевой функции? Я хочу документировать свой алгоритм. ? Могу ли я просто добавить вес e_ {H} * и вес e_ {L} * в целевую функцию? @ Джон Хеденгрен. - person Tej Tizaoui; 13.10.2020
comment
Дополнительная информация о целевой функции CV с целью l1-norm: apmonitor.com/ сделать / index.php / Main / ControllerObjective - person John Hedengren; 13.10.2020