Python 能量系统的 GEKKO MINLP 优化：如何构建二维数组的中间体
Python GEKKO MINLP optimization of energy system: How to build intermediates that are 2D arrays
我目前正在 Python GEKKO 中实施 MINLP 优化问题，以确定三联产能源系统的最佳运行策略。由于我将不同代表性日期的所有时段的能源需求视为输入数据，因此基本上我所有的决策变量、中间值等都是二维数组。
我怀疑二维中间体的声明是我的问题。现在我使用列表理解来声明 2D 中间体，但似乎 python 不能使用这些中间体。此外，错误 This steady-state IMODE only allows scalar values. 发生。
每当我像这样使用 GEKKO m.Array 函数时：
e_GT = m.Array(m.Intermediate(E_GT[z][p]/E_max_GT) for z in range(Z) for p in range(P), (Z,P))
它说，无法调用 GEKKO 对象 m.Intermediate。
如果有人能给我提示，我将不胜感激。
完整代码如下：
"""
Created on Fri Nov 22 10:18:33 2019
@author: julia
"""
# __Get GEKKO & numpy___
from gekko import GEKKO
import numpy as np
# ___Initialize model___
m = GEKKO()
# ___Global options_____
m.options.SOLVER = 1 # APOPT is MINLP Solver
# ______Constants_______
i = m.Const(value=0.05)
n = m.Const(value=10)
C_GT = m.Const(value=100000)
C_RB = m.Const(value=10000)
C_HB = m.Const(value=10000)
C_RS = m.Const(value=10000)
C_RE = m.Const(value=10000)
Z = 12
P = 24
E_min_GT = m.Const(value=1000)
E_max_GT = m.Const(value=50000)
F_max_GT = m.Const(value=100000)
Q_max_GT = m.Const(value=100000)
a = m.Const(value=1)
b = m.Const(value=1)
c = m.Const(value=1)
d = m.Const(value=1)
eta_RB = m.Const(value=0.01)
eta_HB = m.Const(value=0.01)
eta_RS = m.Const(value=0.01)
eta_RE = m.Const(value=0.01)
alpha = m.Const(value=0.01)
# ______Parameters______
T_z = m.Param([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
C_p_Gas = m.Param(np.ones([P]))
C_p_Elec = m.Param(np.ones([P]))
E_d = np.ones([Z,P])
H_d = np.ones([Z,P])
K_d = np.ones([Z,P])
# _______Variables______
E_purch = m.Array(m.Var, (Z,P), lb=0)
E_GT = m.Array(m.Var, (Z,P), lb=0)
F_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT_RB = m.Array(m.Var, (Z,P), lb=0)
Q_disp = m.Array(m.Var, (Z,P), lb=0)
Q_HB = m.Array(m.Var, (Z,P), lb=0)
K_RS = m.Array(m.Var, (Z,P), lb=0)
K_RE = m.Array(m.Var, (Z,P), lb=0)
delta_GT = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_HB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RS = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RE = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
# ____Intermediates_____
R = m.Intermediate((i*(1+i)**n)/((1+i)**n-1))
e_min_GT = m.Intermediate(E_min_GT/E_max_GT)
e_GT = [m.Intermediate(E_GT[z][p]/E_max_GT) for z in range(Z) for p in range(P)]
f_GT = [m.Intermediate(F_GT[z][p]/F_max_GT) for z in range(Z) for p in range(P)]
q_GT = [m.Intermediate(Q_GT[z][p]/Q_max_GT) for z in range(Z) for p in range(P)]
Q_RB = [m.Intermediate(eta_RB*Q_GT_RB[z][p]*delta_RB[z][p]) for z in range(Z) for p in range(P)]
F_HB = [m.Intermediate(eta_HB*Q_HB[z][p]*delta_HB[z][p]) for z in range(Z) for p in range(P)]
Q_RS = [m.Intermediate(eta_RS*K_RS[z][p]*delta_RS[z][p]) for z in range(Z) for p in range(P)]
E_RE = [m.Intermediate(eta_RE*K_RE[z][p]*delta_RE[z][p]) for z in range(Z) for p in range(P)]
F_Gas = [m.Intermediate(F_GT[z][p] + eta_HB*Q_HB[z][p]*delta_HB[z][p]) for z in range(Z) for p in range(P)]
Cc = m.Intermediate(R*(C_GT + C_RB + C_HB + C_RS + C_RE))
Cr_z = m.Intermediate((sum(C_p_Gas[p]*F_Gas[z][p] + C_p_Elec[p]*E_purch[z][p]) for p in range(P)) for z in range(Z))
Cr = m.Intermediate(sum(Cr_z[z]*T_z[z]) for z in range(Z))
# ______Equations_______
m.Equation(e_min_GT[z][p]*delta_GT[z][p] <= e_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(e_GT[z][p] <= 1*delta_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(f_GT [z][p]== a*delta_GT[z][p] + b*e_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(q_GT [z][p]== c*delta_GT[z][p] + d*e_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(E_purch[z][p] + E_GT[z][p] == E_RE[z][p] + E_d[z][p] for z in range(Z) for p in range(P))
m.Equation(Q_GT[z][p] == Q_disp[z][p] + Q_GT_RB[z][p] for z in range(Z) for p in range(P))
m.Equation(Q_RB[z][p] + Q_HB[z][p] == Q_RS[z][p] + H_d[z][p] for z in range(Z) for p in range(P))
m.Equation(K_RS[z][p] + K_RE[z][p] == K_d[z][p] for z in range(Z) for p in range(P))
m.Equation(Q_disp[z][p] <= alpha*Q_GT[z][p] for z in range(Z) for p in range(P))
# ______Objective_______
m.Obj(Cc + Cr)
#_____Solve Problem_____
m.solve()
为了诊断问题，我在解决命令之前添加了打开 运行 文件夹的调用。
#_____Solve Problem_____
m.open_folder()
m.solve()
我用文本编辑器打开了 gk_model0.apm 模型文件来查看模型的文本版本。最下面显示最后两个中间体和9个方程有问题
i2327=<generator object <genexpr> at 0x0E51BC30>
i2328=<generator object <genexpr> at 0x0E51BC30>
End Intermediates
Equations
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
minimize (i2326+i2328)
End Equations
End Model
当 Param 或 Var 的 .value 属性 是列表或 numpy 数组而不是标量值时，就会发生这种情况。
# ______Parameters______
#T_z = m.Param([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
T_z = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
还有一些其他问题，例如
e_min_GT 在第一个等式中被引用为列表的元素而不是标量值- 中间列表理解中缺少方括号，使列表成为二维的，例如
[[2 for p in range(P)] for z in range(Z)]
- 尺寸不匹配
for z in range(Z)] for p in range(P)] 而不是 for p in range(P)] for z in range(Z)]
- 为了提高效率，请使用
m.sum 而不是 sum
我还进行了其他一些更改。差分程序应该显示它们。
"""
Created on Fri Nov 22 10:18:33 2019
@author: julia
"""
# __Get GEKKO & numpy___
from gekko import GEKKO
import numpy as np
# ___Initialize model___
m = GEKKO()
# ___Global options_____
m.options.SOLVER = 1 # APOPT is MINLP Solver
# ______Constants_______
i = m.Const(value=0.05)
n = m.Const(value=10)
C_GT = m.Const(value=100000)
C_RB = m.Const(value=10000)
C_HB = m.Const(value=10000)
C_RS = m.Const(value=10000)
C_RE = m.Const(value=10000)
Z = 12
P = 24
E_min_GT = m.Const(value=1000)
E_max_GT = m.Const(value=50000)
F_max_GT = m.Const(value=100000)
Q_max_GT = m.Const(value=100000)
a = m.Const(value=1)
b = m.Const(value=1)
c = m.Const(value=1)
d = m.Const(value=1)
eta_RB = m.Const(value=0.01)
eta_HB = m.Const(value=0.01)
eta_RS = m.Const(value=0.01)
eta_RE = m.Const(value=0.01)
alpha = m.Const(value=0.01)
# ______Parameters______
T_z = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
C_p_Gas = np.ones(P)
C_p_Elec = np.ones(P)
E_d = np.ones([Z,P])
H_d = np.ones([Z,P])
K_d = np.ones([Z,P])
# _______Variables______
E_purch = m.Array(m.Var, (Z,P), lb=0)
E_GT = m.Array(m.Var, (Z,P), lb=0)
F_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT_RB = m.Array(m.Var, (Z,P), lb=0)
Q_disp = m.Array(m.Var, (Z,P), lb=0)
Q_HB = m.Array(m.Var, (Z,P), lb=0)
K_RS = m.Array(m.Var, (Z,P), lb=0)
K_RE = m.Array(m.Var, (Z,P), lb=0)
delta_GT = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_HB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RS = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RE = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
# ____Intermediates_____
R = m.Intermediate((i*(1+i)**n)/((1+i)**n-1))
e_min_GT = m.Intermediate(E_min_GT/E_max_GT)
e_GT = [[m.Intermediate(E_GT[z][p]/E_max_GT) for p in range(P)] for z in range(Z)]
f_GT = [[m.Intermediate(F_GT[z][p]/F_max_GT) for p in range(P)] for z in range(Z)]
q_GT = [[m.Intermediate(Q_GT[z][p]/Q_max_GT) for p in range(P)] for z in range(Z)]
Q_RB = [[m.Intermediate(eta_RB*Q_GT_RB[z][p]*delta_RB[z][p]) for p in range(P)] for z in range(Z)]
F_HB = [[m.Intermediate(eta_HB*Q_HB[z][p]*delta_HB[z][p]) for p in range(P)] for z in range(Z)]
Q_RS = [[m.Intermediate(eta_RS*K_RS[z][p]*delta_RS[z][p]) for p in range(P)] for z in range(Z)]
E_RE = [[m.Intermediate(eta_RE*K_RE[z][p]*delta_RE[z][p]) for p in range(P)] for z in range(Z)]
F_Gas = [[m.Intermediate(F_GT[z][p] + eta_HB*Q_HB[z][p]*delta_HB[z][p]) for p in range(P)] for z in range(Z)]
Cc = m.Intermediate(R*(C_GT + C_RB + C_HB + C_RS + C_RE))
Cr_z = [m.Intermediate(m.sum([C_p_Gas[p]*F_Gas[z][p] + C_p_Elec[p]*E_purch[z][p] for p in range(P)])) for z in range(Z)]
Cr = m.Intermediate(m.sum([Cr_z[z]*T_z[z] for z in range(Z)]))
# ______Equations_______
m.Equation([e_min_GT*delta_GT[z][p] <= e_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([e_GT[z][p] <= 1*delta_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([f_GT [z][p]== a*delta_GT[z][p] + b*e_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([q_GT [z][p]== c*delta_GT[z][p] + d*e_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([E_purch[z][p] + E_GT[z][p] == E_RE[z][p] + E_d[z][p] for z in range(Z) for p in range(P)])
m.Equation([Q_GT[z][p] == Q_disp[z][p] + Q_GT_RB[z][p] for z in range(Z) for p in range(P)])
m.Equation([Q_RB[z][p] + Q_HB[z][p] == Q_RS[z][p] + H_d[z][p] for z in range(Z) for p in range(P)])
m.Equation([K_RS[z][p] + K_RE[z][p] == K_d[z][p] for z in range(Z) for p in range(P)])
m.Equation([Q_disp[z][p] <= alpha*Q_GT[z][p] for z in range(Z) for p in range(P)])
# ______Objective_______
m.Obj(Cc + Cr)
#_____Solve Problem_____
#m.open_folder()
m.solve()
问题求解时间为 1.6 sec。
--------- APM Model Size ------------
Each time step contains
Objects : 13
Constants : 20
Variables : 5209
Intermediates: 2320
Connections : 313
Equations : 5213
Residuals : 2893
Number of state variables: 5209
Number of total equations: - 2905
Number of slack variables: - 864
---------------------------------------
Degrees of freedom : 1440
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 1.53 NLPi: 3 Dpth: 0 Lvs: 0 Obj: 2.69E+04 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 1.59730000000854 sec
Objective : 26890.6404951639
Successful solution
---------------------------------------------------
二维列表定义需要额外的方括号。这给出了一个 3 行和 4 列的二维列表。
[[p+10*z for p in range(3)] for z in range(4)]
# Result: [[0, 1, 2], [10, 11, 12], [20, 21, 22], [30, 31, 32]]
如果省略内部括号，它是一个长度为 12 的一维列表。
[p+10*z for p in range(3) for z in range(4)]
# Result: [0, 10, 20, 30, 1, 11, 21, 31, 2, 12, 22, 32]
当列表的每个元素都是 Gekko 时它也有效 Intermediate。
[[m.Intermediate(p+10*z) for p in range(3)] for z in range(4)]
我目前正在 Python GEKKO 中实施 MINLP 优化问题，以确定三联产能源系统的最佳运行策略。由于我将不同代表性日期的所有时段的能源需求视为输入数据，因此基本上我所有的决策变量、中间值等都是二维数组。 我怀疑二维中间体的声明是我的问题。现在我使用列表理解来声明 2D 中间体，但似乎 python 不能使用这些中间体。此外，错误 This steady-state IMODE only allows scalar values. 发生。
每当我像这样使用 GEKKO m.Array 函数时：
e_GT = m.Array(m.Intermediate(E_GT[z][p]/E_max_GT) for z in range(Z) for p in range(P), (Z,P))
它说，无法调用 GEKKO 对象 m.Intermediate。
如果有人能给我提示，我将不胜感激。
完整代码如下：
"""
Created on Fri Nov 22 10:18:33 2019
@author: julia
"""
# __Get GEKKO & numpy___
from gekko import GEKKO
import numpy as np
# ___Initialize model___
m = GEKKO()
# ___Global options_____
m.options.SOLVER = 1 # APOPT is MINLP Solver
# ______Constants_______
i = m.Const(value=0.05)
n = m.Const(value=10)
C_GT = m.Const(value=100000)
C_RB = m.Const(value=10000)
C_HB = m.Const(value=10000)
C_RS = m.Const(value=10000)
C_RE = m.Const(value=10000)
Z = 12
P = 24
E_min_GT = m.Const(value=1000)
E_max_GT = m.Const(value=50000)
F_max_GT = m.Const(value=100000)
Q_max_GT = m.Const(value=100000)
a = m.Const(value=1)
b = m.Const(value=1)
c = m.Const(value=1)
d = m.Const(value=1)
eta_RB = m.Const(value=0.01)
eta_HB = m.Const(value=0.01)
eta_RS = m.Const(value=0.01)
eta_RE = m.Const(value=0.01)
alpha = m.Const(value=0.01)
# ______Parameters______
T_z = m.Param([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
C_p_Gas = m.Param(np.ones([P]))
C_p_Elec = m.Param(np.ones([P]))
E_d = np.ones([Z,P])
H_d = np.ones([Z,P])
K_d = np.ones([Z,P])
# _______Variables______
E_purch = m.Array(m.Var, (Z,P), lb=0)
E_GT = m.Array(m.Var, (Z,P), lb=0)
F_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT_RB = m.Array(m.Var, (Z,P), lb=0)
Q_disp = m.Array(m.Var, (Z,P), lb=0)
Q_HB = m.Array(m.Var, (Z,P), lb=0)
K_RS = m.Array(m.Var, (Z,P), lb=0)
K_RE = m.Array(m.Var, (Z,P), lb=0)
delta_GT = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_HB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RS = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RE = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
# ____Intermediates_____
R = m.Intermediate((i*(1+i)**n)/((1+i)**n-1))
e_min_GT = m.Intermediate(E_min_GT/E_max_GT)
e_GT = [m.Intermediate(E_GT[z][p]/E_max_GT) for z in range(Z) for p in range(P)]
f_GT = [m.Intermediate(F_GT[z][p]/F_max_GT) for z in range(Z) for p in range(P)]
q_GT = [m.Intermediate(Q_GT[z][p]/Q_max_GT) for z in range(Z) for p in range(P)]
Q_RB = [m.Intermediate(eta_RB*Q_GT_RB[z][p]*delta_RB[z][p]) for z in range(Z) for p in range(P)]
F_HB = [m.Intermediate(eta_HB*Q_HB[z][p]*delta_HB[z][p]) for z in range(Z) for p in range(P)]
Q_RS = [m.Intermediate(eta_RS*K_RS[z][p]*delta_RS[z][p]) for z in range(Z) for p in range(P)]
E_RE = [m.Intermediate(eta_RE*K_RE[z][p]*delta_RE[z][p]) for z in range(Z) for p in range(P)]
F_Gas = [m.Intermediate(F_GT[z][p] + eta_HB*Q_HB[z][p]*delta_HB[z][p]) for z in range(Z) for p in range(P)]
Cc = m.Intermediate(R*(C_GT + C_RB + C_HB + C_RS + C_RE))
Cr_z = m.Intermediate((sum(C_p_Gas[p]*F_Gas[z][p] + C_p_Elec[p]*E_purch[z][p]) for p in range(P)) for z in range(Z))
Cr = m.Intermediate(sum(Cr_z[z]*T_z[z]) for z in range(Z))
# ______Equations_______
m.Equation(e_min_GT[z][p]*delta_GT[z][p] <= e_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(e_GT[z][p] <= 1*delta_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(f_GT [z][p]== a*delta_GT[z][p] + b*e_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(q_GT [z][p]== c*delta_GT[z][p] + d*e_GT[z][p] for z in range(Z) for p in range(P))
m.Equation(E_purch[z][p] + E_GT[z][p] == E_RE[z][p] + E_d[z][p] for z in range(Z) for p in range(P))
m.Equation(Q_GT[z][p] == Q_disp[z][p] + Q_GT_RB[z][p] for z in range(Z) for p in range(P))
m.Equation(Q_RB[z][p] + Q_HB[z][p] == Q_RS[z][p] + H_d[z][p] for z in range(Z) for p in range(P))
m.Equation(K_RS[z][p] + K_RE[z][p] == K_d[z][p] for z in range(Z) for p in range(P))
m.Equation(Q_disp[z][p] <= alpha*Q_GT[z][p] for z in range(Z) for p in range(P))
# ______Objective_______
m.Obj(Cc + Cr)
#_____Solve Problem_____
m.solve()
为了诊断问题，我在解决命令之前添加了打开 运行 文件夹的调用。
#_____Solve Problem_____
m.open_folder()
m.solve()
我用文本编辑器打开了 gk_model0.apm 模型文件来查看模型的文本版本。最下面显示最后两个中间体和9个方程有问题
i2327=<generator object <genexpr> at 0x0E51BC30>
i2328=<generator object <genexpr> at 0x0E51BC30>
End Intermediates
Equations
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
<generator object <genexpr> at 0x0E51BC30>
minimize (i2326+i2328)
End Equations
End Model
当 Param 或 Var 的 .value 属性 是列表或 numpy 数组而不是标量值时，就会发生这种情况。
# ______Parameters______
#T_z = m.Param([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
T_z = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
还有一些其他问题，例如
e_min_GT在第一个等式中被引用为列表的元素而不是标量值- 中间列表理解中缺少方括号，使列表成为二维的，例如
[[2 for p in range(P)] for z in range(Z)] - 尺寸不匹配
for z in range(Z)] for p in range(P)]而不是for p in range(P)] for z in range(Z)] - 为了提高效率，请使用
m.sum而不是sum
我还进行了其他一些更改。差分程序应该显示它们。
"""
Created on Fri Nov 22 10:18:33 2019
@author: julia
"""
# __Get GEKKO & numpy___
from gekko import GEKKO
import numpy as np
# ___Initialize model___
m = GEKKO()
# ___Global options_____
m.options.SOLVER = 1 # APOPT is MINLP Solver
# ______Constants_______
i = m.Const(value=0.05)
n = m.Const(value=10)
C_GT = m.Const(value=100000)
C_RB = m.Const(value=10000)
C_HB = m.Const(value=10000)
C_RS = m.Const(value=10000)
C_RE = m.Const(value=10000)
Z = 12
P = 24
E_min_GT = m.Const(value=1000)
E_max_GT = m.Const(value=50000)
F_max_GT = m.Const(value=100000)
Q_max_GT = m.Const(value=100000)
a = m.Const(value=1)
b = m.Const(value=1)
c = m.Const(value=1)
d = m.Const(value=1)
eta_RB = m.Const(value=0.01)
eta_HB = m.Const(value=0.01)
eta_RS = m.Const(value=0.01)
eta_RE = m.Const(value=0.01)
alpha = m.Const(value=0.01)
# ______Parameters______
T_z = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
C_p_Gas = np.ones(P)
C_p_Elec = np.ones(P)
E_d = np.ones([Z,P])
H_d = np.ones([Z,P])
K_d = np.ones([Z,P])
# _______Variables______
E_purch = m.Array(m.Var, (Z,P), lb=0)
E_GT = m.Array(m.Var, (Z,P), lb=0)
F_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT = m.Array(m.Var, (Z,P), lb=0)
Q_GT_RB = m.Array(m.Var, (Z,P), lb=0)
Q_disp = m.Array(m.Var, (Z,P), lb=0)
Q_HB = m.Array(m.Var, (Z,P), lb=0)
K_RS = m.Array(m.Var, (Z,P), lb=0)
K_RE = m.Array(m.Var, (Z,P), lb=0)
delta_GT = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_HB = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RS = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
delta_RE = m.Array(m.Var, (Z,P), lb=0, ub=1, integer=True)
# ____Intermediates_____
R = m.Intermediate((i*(1+i)**n)/((1+i)**n-1))
e_min_GT = m.Intermediate(E_min_GT/E_max_GT)
e_GT = [[m.Intermediate(E_GT[z][p]/E_max_GT) for p in range(P)] for z in range(Z)]
f_GT = [[m.Intermediate(F_GT[z][p]/F_max_GT) for p in range(P)] for z in range(Z)]
q_GT = [[m.Intermediate(Q_GT[z][p]/Q_max_GT) for p in range(P)] for z in range(Z)]
Q_RB = [[m.Intermediate(eta_RB*Q_GT_RB[z][p]*delta_RB[z][p]) for p in range(P)] for z in range(Z)]
F_HB = [[m.Intermediate(eta_HB*Q_HB[z][p]*delta_HB[z][p]) for p in range(P)] for z in range(Z)]
Q_RS = [[m.Intermediate(eta_RS*K_RS[z][p]*delta_RS[z][p]) for p in range(P)] for z in range(Z)]
E_RE = [[m.Intermediate(eta_RE*K_RE[z][p]*delta_RE[z][p]) for p in range(P)] for z in range(Z)]
F_Gas = [[m.Intermediate(F_GT[z][p] + eta_HB*Q_HB[z][p]*delta_HB[z][p]) for p in range(P)] for z in range(Z)]
Cc = m.Intermediate(R*(C_GT + C_RB + C_HB + C_RS + C_RE))
Cr_z = [m.Intermediate(m.sum([C_p_Gas[p]*F_Gas[z][p] + C_p_Elec[p]*E_purch[z][p] for p in range(P)])) for z in range(Z)]
Cr = m.Intermediate(m.sum([Cr_z[z]*T_z[z] for z in range(Z)]))
# ______Equations_______
m.Equation([e_min_GT*delta_GT[z][p] <= e_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([e_GT[z][p] <= 1*delta_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([f_GT [z][p]== a*delta_GT[z][p] + b*e_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([q_GT [z][p]== c*delta_GT[z][p] + d*e_GT[z][p] for z in range(Z) for p in range(P)])
m.Equation([E_purch[z][p] + E_GT[z][p] == E_RE[z][p] + E_d[z][p] for z in range(Z) for p in range(P)])
m.Equation([Q_GT[z][p] == Q_disp[z][p] + Q_GT_RB[z][p] for z in range(Z) for p in range(P)])
m.Equation([Q_RB[z][p] + Q_HB[z][p] == Q_RS[z][p] + H_d[z][p] for z in range(Z) for p in range(P)])
m.Equation([K_RS[z][p] + K_RE[z][p] == K_d[z][p] for z in range(Z) for p in range(P)])
m.Equation([Q_disp[z][p] <= alpha*Q_GT[z][p] for z in range(Z) for p in range(P)])
# ______Objective_______
m.Obj(Cc + Cr)
#_____Solve Problem_____
#m.open_folder()
m.solve()
问题求解时间为 1.6 sec。
--------- APM Model Size ------------
Each time step contains
Objects : 13
Constants : 20
Variables : 5209
Intermediates: 2320
Connections : 313
Equations : 5213
Residuals : 2893
Number of state variables: 5209
Number of total equations: - 2905
Number of slack variables: - 864
---------------------------------------
Degrees of freedom : 1440
----------------------------------------------
Steady State Optimization with APOPT Solver
----------------------------------------------
Iter: 1 I: 0 Tm: 1.53 NLPi: 3 Dpth: 0 Lvs: 0 Obj: 2.69E+04 Gap: 0.00E+00
Successful solution
---------------------------------------------------
Solver : APOPT (v1.0)
Solution time : 1.59730000000854 sec
Objective : 26890.6404951639
Successful solution
---------------------------------------------------
二维列表定义需要额外的方括号。这给出了一个 3 行和 4 列的二维列表。
[[p+10*z for p in range(3)] for z in range(4)]
# Result: [[0, 1, 2], [10, 11, 12], [20, 21, 22], [30, 31, 32]]
如果省略内部括号，它是一个长度为 12 的一维列表。
[p+10*z for p in range(3) for z in range(4)]
# Result: [0, 10, 20, 30, 1, 11, 21, 31, 2, 12, 22, 32]
当列表的每个元素都是 Gekko 时它也有效 Intermediate。
[[m.Intermediate(p+10*z) for p in range(3)] for z in range(4)]