Battery Storage modelling for a data centre

Question

I am trying to build an optimisation model (minimise costs) to procure electricity on an hourly basis to a data centre using wind,solar,grid and battery storage with an objective to achieve a certain carbon-free energy ratio (procured low-carbon electricity as a % of total data centre load) . The model works with wind,solar and grid but I keep having error messages when I introduce storage as a variable/constraint

The idea is to store the batteries when I have excess wind & solar compared to load and then to discharge when grid is used in order to reduce carbon emissions.

I have SOC that is shown negative in the results while I have included in the constraints to be positive

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)