使用 Gekko 的混合整数模型预测控制

Question

问题陈述：

尝试利用 GEKKO 优化 select 可用容量内的电器数量，同时最大限度地降低能源成本。

问题模型：

我不能 post 图片，但 here 是 link 来查看公式，您也可以复制并粘贴到下面以在视觉上更好地找到它：

\min \sum _{t=1}^{24} n x_{i} P_{i}

受制于：

\sum _{i=1}^{24} C_{i} - x_{i} = 0 
C_{min} < C_{i}< C_{max}
x_{min} < x_{i}< x_{max}
n_{min} < n_{i}< n_{max}
\sum_{t=1}^{T} \tau_{i}^{t} = I_{i}

其中， x和n是表示赋值的状态变量一个站内电器的功率和数量。 C是容量，是工作状态设备（ON/OFF 条件）在时间段 t 和 I_{i}是每个电器允许的间隔时间（即 5 小时）

目前完成的工作：

感谢 post 我整理了以下内容：

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
cap_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.]
cap_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,50.,50.,50.,50.,50.,\
       50.,50.,50.,50.,50.,50.,50.,50.,50.,0.05,0.05,0.05]

k = 1.1 #Safety factor
CL = m.Param(value=cap_low)
CH = m.Param(value=cap_upper)
TOU = m.Param(value=TOU)

x = m.MV(lb=6, ub=32)
x.STATUS = 1  # allow optimizer to change

n = m.MV(lb=1, ub=10, integer=True)
n.STATUS = 1  # allow optimizer to change

# Controlled Variable
C = m.SV(value=55)
m.Equations([C>=CL, C<=CH])

m.Equation(k*C - n* x == 0)
m.Minimize(TOU*n*x)

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

plt.subplot(3,1,1)
plt.plot(m.time,cap_low,'k--')
plt.plot(m.time,cap_upper,'k--')
plt.plot(m.time,C.value,'r-')
plt.ylabel('Demand')
plt.subplot(3,1,2)
plt.step(m.time,x.value,'b:')
plt.ylabel('Load')
plt.xlabel('Time (hr)')
plt.subplot(3,1,3)
plt.step(m.time,n.value,'r:')
plt.ylabel('No. of Appliances')
plt.xlabel('Time (hr)')
plt.show()

结果

上面的代码产生了一个有希望的结果，但是缺少限制每个设备运行时间的最后一个约束\sum_{t=1}^{T} \tau_{i}^{t} = I_{i}。

问题：

这个问题看起来像是模型预测控制（跟踪可用负载能力）和混合整数线性规划（具有多个线性运算方程的多个设备）的组合。所以问题是：

如何在上面的代码中定义每个设备的允许间隔以限制它们在一天内的运行。
在 Gekko 中集成 MPC 和 MILP 的最佳方式是什么

Answer 1

问题 2 的答案：切换到 APOPT 求解器来解决 MILP 或 MINLP 问题。

问题 1 的答案：将 n 分隔为单独的设备状态：

# decide when to turn on each appliance
na = 20 # number of appliances
a = m.Array(m.MV,na,lb=0, ub=1, integer=True)
for ai in a:
    ai.STATUS = 1
    ai.DCOST = 0

# number of appliances running
n = m.Var(lb=1, ub=5)
m.Equation(n==m.sum(a))

底部的子图显示各个电器何时打开或关闭。将设备限制为每天 5 小时会导致收敛问题：

ax = [m.Intermediate(m.integral(ai)) for ai in a]
for ai in ax:
    # maximum of 5 hours each day
    m.Minimize(1e-4*ax**2)
    m.Equation(ax<=5)

也许您可以从这段代码开始并尝试其他策略来提高收敛性，例如不同的 APOPT solver options 或约束形式。创建一个软约束（用经济激励而不是硬限制来惩罚）也有助于提高收敛性。

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

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

#initialize variables
cap_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.]
cap_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,50.,50.,50.,50.,50.,\
       50.,50.,50.,50.,50.,50.,50.,50.,50.,0.05,0.05,0.05]

k = 1.1 #Safety factor
CL = m.Param(value=cap_low)
CH = m.Param(value=cap_upper)
TOU = m.Param(value=TOU)

x = m.MV(lb=6, ub=32)
x.STATUS = 1  # allow optimizer to change

# decide when to turn on each appliance
na = 20 # number of appliances
a = m.Array(m.MV,na,lb=0, ub=1, integer=True)
for ai in a:
    ai.STATUS = 1
    ai.DCOST = 0

ax = [m.Intermediate(m.integral(ai)) for ai in a]
for ai in ax:
    # maximum of 5 hours each day
for ai in ax:
    # maximum of 5 hours each day
    m.Minimize(1e-4*ai**2)
    m.Equation(ai<=5)

# number of appliances running
n = m.Var(lb=1, ub=5)
m.Equation(n==m.sum(a))

# Controlled Variable
C = m.SV(value=55)
m.Equations([C>=CL, C<=CH])

m.Equation(k*C - n* x == 0)
m.Minimize(TOU*n*x)

m.options.IMODE = 6
m.options.NODES = 2
m.options.SOLVER = 1

m.solver_options = ['minlp_gap_tol 1.0e-2',\
                    'minlp_maximum_iterations 10000',\
                    'minlp_max_iter_with_int_sol 500',\
                    'minlp_branch_method 3']

m.solve(disp=True)

plt.subplot(4,1,1)
plt.plot(m.time,cap_low,'k--')
plt.plot(m.time,cap_upper,'k--')
plt.plot(m.time,C.value,'r-')
plt.ylabel('Demand')
plt.subplot(4,1,2)
plt.step(m.time,x.value,'b:')
plt.ylabel('Load')
plt.xlabel('Time (hr)')
plt.subplot(4,1,3)
plt.step(m.time,n.value,'r:')
plt.ylabel('No. of Appliances')
plt.subplot(4,1,4)
for ai in a:
    plt.plot(m.time,ai.value)
plt.xlabel('Time (hr)')
plt.show()

使用 Gekko 的混合整数模型预测控制

Mixed-Integer Model Predictive Control using Gekko

python

mathematical-optimization

gekko

mixed-integer-programming