数据中心电池储能优化建模求助：SOC负值异常与约束报错

我正在搭建一个以成本最小化为目标的小时级购电优化模型，为数据中心提供风电、光伏、电网及电池储能四类供电方式，核心要求是达成指定的**无碳能源比率（CFE）**目标——即购得的低碳电量占总负荷的比例。

模型在仅包含风电、光伏、电网三类电源时能正常运行，但加入储能变量与约束后持续报错，且优化结果中荷电状态（SOC）出现负值，尽管已经添加了SOC非负约束。

模型设计逻辑：当风电、光伏发电量超过负荷时为电池充电；使用电网供电时放电，以此降低碳排放。

原代码 #

H = 1000

import numpy as np
import pandas as pd
from gurobipy import *
import matplotlib.pyplot as plt
import seaborn as sns
import math 

TEST = pd.read_excel('/Users/charles-henryduprez/Desktop/project/Datasets/Test2031.xlsx')
TEST_matrix = TEST.to_numpy()
Load = TEST_matrix[:H, 0]
Wind = TEST_matrix[:H, 1]
Solar = TEST_matrix[:H, 2]
GridCFE = TEST_matrix[:H, 3]
GridPrice = TEST_matrix[:H, 4]
OnlyGridCosts = [x * y for x, y in zip(Load, GridPrice)]
AverageOnlyGridCosts = sum(OnlyGridCosts)/len(OnlyGridCosts)
COSTWIND = 70
COSTSOLAR = 70
target = 0.90  # Define CFE1 target

# Define storage
MW_Storage = 50
Hours_MW = 2
MWh_Storage_Max = MW_Storage * Hours_MW
Capex_MW_Storage = 600000
Opex_MW_Storage = 30000
Period = H / 8760
Battery_eff = 0.87
Efficiency = math.sqrt(Battery_eff)
Storage_lifetime = 12 * Period

# Define the SOC 
SOC = {}  # Dictionary to store the SOC values
SOC[0] = 0.5 # Initialize the initial storage value
SOC_values =[]

# Define the objective
def objfunc(x1, x2, y, s):
    return quicksum(
        x1 * Wind[i] * COSTWIND + x2 * Solar[i] * COSTSOLAR + GridPrice[i] * y[i] + ((MW_Storage * Capex_MW_Storage) / Storage_lifetime)
        for i in range(H))

# Define the constraints
def c1(x1, x2, y, z, s):
    return [
        x1 * Wind[i] + x2 * Solar[i] + y[i] - z[i] + MWh_Storage_Max * s[i] * Efficiency - Load[i]
        for i in range(H)]

def c2(x1, x2, y, z):
    return (
        quicksum(x1 * Wind[i] + x2 * Solar[i] + y[i] * GridCFE[i] - z[i] for i in range(H))
        - target * quicksum(Load))

def c3(s, model):
    SOC[1] = 0  # Initialize SOC[1] to 0

    constraints = []
    for i in range(H):
        SOC[i + 1] = SOC[i] + s[i]  # Update the SOC[i+1] value
        constraints.append(model.addConstr(SOC[i + 1] >= 0, name="c3_lb_" + str(i)))
        constraints.append(model.addConstr(SOC[i + 1] <= 1, name="c3_ub_" + str(i)))
    return constraints


# Create Gurobi model
model = Model()

# Create variables
x1 = model.addVar(lb=0, ub=1000, name="x1")  # Number of MWs of Wind
x2 = model.addVar(lb=0, ub=1000, name="x2")  # Number of MWs of Solar
y = model.addVars(H, lb=0, ub=1000, name="y")  # Grid MWhs
z = model.addVars(H, lb=0, ub=1000, name="z")  # Excess WindSolar MWhs
s = model.addVars(H, lb=-1, ub=1, name="s")  # Storage/Dispatch

# Set objective function
obj = objfunc(x1, x2, y, s)
model.setObjective(obj, GRB.MINIMIZE)

# Add constraints
for i in range(H):
    model.addConstr(c1(x1, x2, y, z, s)[i] == 0, name="c1_" + str(i))

model.addConstr(c2(x1, x2, y, z) == 0, name="c2")

# Call the c3 and c4 functions and pass the model object
c3(s, model)

# Optimize the model
model.optimize()

# Print the results
if model.status == GRB.OPTIMAL:
    print("Optimal solution found")
    print("Objective value:", model.objVal)
    print("Solution:")
    print("x1:", x1.X)
    print("x2:", x2.X)
    for i in range(H):
        print(f"y[{i}]:", y[i].X)
        print(f"z[{i}]:", z[i].X)
        print(f"s[{i}]:", s[i].X)

    # Store the values of x1, x2, y, z, s in lists
    Wind_values = [x1.X * Wind[i] for i in range(H)]
    Solar_values = [x2.X * Solar[i] for i in range(H)]
    Load_values = [Load[i] for i in range(H)]
    y_values = [y[i].X for i in range(H)]
    z_values = [z[i].X for i in range(H)]
    s_values = [s[i].X for i in range(H)]


    # Store the calculation values in lists
    GridCosts_values = [y[i].X * GridPrice[i] for i in range(H)]
    Storage_Discharge_MWhs = MWh_Storage_Max * np.array(s_values)  # Calculate Storage_Discharge_MWhs
    WindCosts = [x1.X * Wind[i] * COSTWIND for i in range(H)]
    SolarCosts = [x2.X * Solar[i] * COSTSOLAR for i in range(H)]
    CFECosts = [WindCosts[i] + SolarCosts[i] + GridCosts_values[i] for i in range(H)]
    ConsumedGridCFE = [y[i].X * GridCFE[i] for i in range(H)]
    CFERatio = [
        (w + s + g) / l if w + s + g < l else 1 for w, s, g, l in
        zip(Wind_values, Solar_values, ConsumedGridCFE, Load)]
    CFEAverage = sum(CFERatio) / len(CFERatio)
    CFECostsAverage = sum(CFECosts) / len(CFECosts)
    Storage_Discharge_MWhs = MWh_Storage_Max * np.array(s_values)
    Discharge_values = [max(0, val) for val in Storage_Discharge_MWhs]      # Separate Storage and Discharge values
    Storage_values = [-min(0, val) for val in Storage_Discharge_MWhs]
    SOC_values = [s[i].X for i in range(H)] # Calculate SOC values

else:
    print("Optimization failed. Status:", model.status)

问题根源与修复方案 #

核心问题 #

1. SOC未定义为优化变量：原代码中SOC是普通Python字典，而非Gurobi变量，导致约束无法被优化器识别，SOC上下限约束形同虚设。
2. 储能成本计算错误：目标函数中把储能 capex 分摊成本重复累加了H次，实际应仅计算一次周期内的分摊值。
3. 能量平衡约束逻辑混乱：充放电效率的应用未区分充电（输入损耗）和放电（输出损耗）场景，导致能量守恒逻辑错误。
4. SOC状态转移逻辑错误：未将充放电量与SOC变化建立正确的约束关联，优化器无法管控SOC的动态变化。

修复后代码 #

H = 1000

import numpy as np
import pandas as pd
from gurobipy import *
import math 

# 数据加载
TEST = pd.read_excel('/Users/charles-henryduprez/Desktop/project/Datasets/Test2031.xlsx')
TEST_matrix = TEST.to_numpy()
Load = TEST_matrix[:H, 0]
Wind = TEST_matrix[:H, 1]
Solar = TEST_matrix[:H, 2]
GridCFE = TEST_matrix[:H, 3]
GridPrice = TEST_matrix[:H, 4]

# 参数定义
COSTWIND = 70
COSTSOLAR = 70
target_cfe = 0.90

# 储能参数
MW_Storage = 50
Hours_MW = 2
MWh_Storage_Max = MW_Storage * Hours_MW  # 最大储能容量
Capex_MW_Storage = 600000
Opex_MW_Storage = 30000
Period = H / 8760  # 模拟周期占全年比例
Battery_charge_eff = 0.87  # 充电效率
Battery_discharge_eff = 0.87  # 放电效率
Storage_lifetime = 12  # 寿命（年）
storage_annual_cost = (MW_Storage * Capex_MW_Storage) / Storage_lifetime + MW_Storage * Opex_MW_Storage
storage_period_cost = storage_annual_cost * Period

# 创建模型
model = Model("DataCenter_Power_Opt")

# 定义变量
x1 = model.addVar(lb=0, ub=1000, name="wind_capacity")  # 风电装机容量（MW）
x2 = model.addVar(lb=0, ub=1000, name="solar_capacity")  # 光伏装机容量（MW）
y = model.addVars(H, lb=0, name="grid_power")  # 电网购电量（MWh）
z = model.addVars(H, lb=0, name="curtailment")  # 弃电量（MWh）
charge = model.addVars(H, lb=0, ub=MW_Storage, name="charge_power")  # 充电功率（MW）
discharge = model.addVars(H, lb=0, ub=MW_Storage, name="discharge_power")  # 放电功率（MW）
SOC = model.addVars(H+1, lb=0, ub=1, name="soc")  # SOC（0-1）

# 初始SOC约束
model.addConstr(SOC[0] == 0.5, name="initial_soc")

# 目标函数：总成本最小化
obj = quicksum(x1 * Wind[i] * COSTWIND + x2 * Solar[i] * COSTSOLAR + GridPrice[i] * y[i] for i in range(H)) + storage_period_cost
model.setObjective(obj, GRB.MINIMIZE)

# 能量平衡约束
for i in range(H):
    model.addConstr(
        x1 * Wind[i] + x2 * Solar[i] + y[i] + discharge[i] * Battery_discharge_eff 
        == Load[i] + charge[i] / Battery_charge_eff + z[i],
        name=f"energy_balance_{i}"
    )

# SOC状态转移约束
for i in range(H):
    model.addConstr(
        SOC[i+1] == SOC[i] + (charge[i] / Battery_charge_eff - discharge[i] * Battery_discharge_eff) / MWh_Storage_Max,
        name=f"soc_transition_{i}"
    )

# CFE目标约束
total_low_carbon_energy = quicksum(x1 * Wind[i] + x2 * Solar[i] + y[i] * GridCFE[i] - z[i] for i in range(H))
model.addConstr(total_low_carbon_energy == target_cfe * quicksum(Load), name="cfe_target")

# 优化求解
model.optimize()

# 结果输出与处理
if model.status == GRB.OPTIMAL:
    print("最优解已找到")
    print(f"总成本: {model.objVal:.2f}")
    print(f"风电装机容量: {x1.X:.2f} MW")
    print(f"光伏装机容量: {x2.X:.2f} MW")

    # 提取结果数据
    wind_gen = [x1.X * Wind[i] for i in range(H)]
    solar_gen = [x2.X * Solar[i] for i in range(H)]
    grid_purch = [y[i].X for i in range(H)]
    charge_vals = [charge[i].X for i in range(H)]
    discharge_vals = [discharge[i].X for i in range(H)]
    soc_vals = [SOC[i].X for i in range(H+1)]
    curtail_vals = [z[i].X for i in range(H)]

    # 计算CFE指标
    consumed_low_carbon = [wind_gen[i] + solar_gen[i] + grid_purch[i] * GridCFE[i] - curtail_vals[i] for i in range(H)]
    cfe_ratio = [consumed_low_carbon[i]/Load[i] if Load[i]>0 else 1 for i in range(H)]
    avg_cfe = sum(cfe_ratio)/H
    print(f"平均CFE比率: {avg_cfe:.2%}")

else:
    print(f"优化失败，状态码: {model.status}")

关键修改说明 #

1. 将SOC定义为Gurobi变量数组，让优化器直接管控SOC的动态变化
2. 拆分充放电变量，分别定义充电/放电功率，明确区分充放电效率的应用场景
3. 修正储能成本计算逻辑，仅计算一次周期内的分摊成本
4. 重新定义SOC状态转移约束，建立充放电量与SOC变化的数学关联
5. 优化能量平衡约束，确保能量守恒逻辑正确
6. 修正CFE计算逻辑，保证无碳电量统计准确

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