GEKKO+Python抽水蓄能建模：充放电互斥二进制变量实现报错

错误根源分析 #

1. 方程传递方式错误：m.Equation()不支持直接传入生成器表达式，需用循环逐个添加或使用m.Equations()（复数形式）处理批量方程。
2. 储能状态方程逻辑错误：原方程用不等式描述储能变化，不符合物理规律，应改为等式，并补充能量转换效率的影响。
3. 二进制约束不完整：仅限制B[i]+T[i]<=1无法约束实际抽水/放水动作，需关联Tur、Bom与二进制变量，确保二进制变量为0时对应动作也为0。
4. 目标函数写法错误：GEKKO数组变量与NumPy数组无法直接做元素级运算，需通过循环或m.sum()实现累加。

修正后的代码 #

import numpy as np
from gekko import GEKKO

horas = 4  # 最大储能（4小时放水时长）
Eff = 0.75  # 能量转换效率
PMD_calculo = np.array([1, 1.5, 3, 0.5, 0.2])  # 时段电价
Num_horas = len(PMD_calculo)

# 初始化GEKKO模型
m = GEKKO(remote=False)
m.options.SOLVER = 1  # 使用APOPT求解MINLP
m.solver_options = ['minlp_integer_tol 0.0']

# 定义变量
Carga = m.Array(m.Var, Num_horas, lb=0, ub=horas)  # 储能状态（0~4小时发电量）
Tur = m.Array(m.Var, Num_horas, lb=0, ub=1)  # 放水发电量（每时段最大1单位）
Bom = m.Array(m.Var, Num_horas, lb=0, ub=1/Eff)  # 抽水耗电量（每时段最大1/Eff单位）
T = m.Array(m.Var, Num_horas, lb=0, ub=1, integer=True)  # 放水二进制变量（1=放水，0=不放水）
B = m.Array(m.Var, Num_horas, lb=0, ub=1, integer=True)  # 抽水二进制变量（1=抽水，0=不抽水）

# 1. 初始储能约束（假设初始储能为0）
m.Equation(Carga[0] == 0)

# 2. 储能状态方程：每时段储能变化 = 抽水增加量 - 放水减少量
for i in range(1, Num_horas):
    # 抽水时，消耗Bom[i-1]电能，转化为Bom[i-1]*Eff的储能；放水时，消耗Tur[i-1]的储能，产生Tur[i-1]电能
    m.Equation(Carga[i] == Carga[i-1] + Bom[i-1]*Eff - Tur[i-1])

# 3. 充放电互斥约束
for i in range(Num_horas):
    m.Equation(B[i] + T[i] <= 1)  # 同一时段不能同时抽水和放水
    # 关联二进制变量与实际动作：二进制为0时，对应动作必须为0
    m.Equation(Tur[i] <= T[i] * 1)  # T=0时Tur必须为0
    m.Equation(Bom[i] <= B[i] * (1/Eff))  # B=0时Bom必须为0

# 4. 目标函数：最大化总收益（发电收益 - 抽水成本）
total_profit = m.sum(Tur[i]*T[i]*PMD_calculo[i] - Bom[i]*PMD_calculo[i]*B[i] for i in range(Num_horas))
m.Maximize(total_profit)

# 求解模型
m.solve(disp=True)

# 输出结果
print("储能状态：", [c.value[0] for c in Carga])
print("放水发电量：", [t.value[0] for t in Tur])
print("抽水耗电量：", [b.value[0] for b in Bom])
print("放水二进制变量：", [t_bin.value[0] for t_bin in T])
print("抽水二进制变量：", [b_bin.value[0] for b_bin in B])
print("总收益：", total_profit.value[0])

关键修正说明 #

• 方程传递：用循环逐个添加储能状态方程和互斥约束，避免生成器表达式导致的解析错误。
• 储能逻辑：将原不等式改为等式，明确储能的连续变化关系，同时加入效率系数Eff，体现抽水时电能到储能的转换损耗。
• 二进制约束强化：新增Tur[i] <= T[i]*1和Bom[i] <= B[i]*(1/Eff)，确保二进制变量与实际动作严格绑定，避免出现二进制为0但仍有充放电的情况。
• 目标函数：用m.sum()遍历每个时段的收益项，实现GEKKO变量与NumPy电价数组的正确运算。

内容的提问来源于stack exchange，提问作者lucassculp