Gurobi giving wrong results (not optimal solution)

here is my code for electricity biding strategy:
from gurobipy import *
# Define your functions

def f(x):

    return x / P_batt_max  # Adjust P_batt_max as needed

def f_inv(x):

    return x * P_batt_max  # Adjust P_batt_max as needed

def sing(x):

    if x >= 0:

        return 1

    else:

        return 0

# Create a new model

model = Model("energy_optimization")

# Decision variables

T = 6

I = 1

P_batt_max = 4

P_pv = {}

P_FlexibleEnergyDemand = {}

P_FlexibleEnergyExcess = {}

P_INFlexibleEnergyDemand = {}

P_INFlexibleEnergyExcess = {}

battCharge ={}

battDischarge ={}

x = {}

y = {}

a = {}

b = {}

SoC = {}

SoC_min = 0.2

SoC_max = 1

# Inflexible energy demand values

D_InflexibleDemand_values = [3, 6, 4, 10, 6, 18]

D_InflexibleDemand = {(t, i): D_InflexibleDemand_values[t] for t in range(T) for i in range(I)}

# Solar PV generation values

PV_generation_values = [5, 6, 4, 6, 6, 8]

PV_generation = {(t, i): PV_generation_values[t] for t in range(T) for i in range(I)}

# Cost parameters

C_lowBidImport = 5

C_HighBidImport = 18

C_lowOfferPrice = 3

C_HighOfferPrice = 4

for t in range(T):

    for i in range(I):

        P_pv[t, i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"P_pv_{t}_{i}")

        P_FlexibleEnergyDemand[t, i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"P_FlexibleEnergyDemand_{t}_{i}")

        P_FlexibleEnergyExcess[t, i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"P_FlexibleEnergyExcess_{t}_{i}")

        P_INFlexibleEnergyDemand[t, i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"P_INFlexibleEnergyDemand_{t}_{i}")

        P_INFlexibleEnergyExcess[t, i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"P_INFlexibleEnergyExcess_{t}_{i}")

       

        battCharge[t,i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"battCharge_{t}_{i}")

        battDischarge[t,i] = model.addVar(lb=0, vtype=GRB.CONTINUOUS, name=f"battDischarge_{t}_{i}")

        x[t, i] = model.addVar(vtype=GRB.BINARY, name=f"x_{t}_{i}")

        y[t, i] = model.addVar(vtype=GRB.BINARY, name=f"y_{t}_{i}")

        a[t, i] = model.addVar(vtype=GRB.BINARY, name=f"a_{t}_{i}")

        b[t, i] = model.addVar(vtype=GRB.BINARY, name=f"b_{t}_{i}")

        # Constraint: Set PV generation to given values

       

for t in range(T+1):

    for i in range(I):        

        SoC[t, i] = model.addVar(lb=0.2, ub=1, vtype=GRB.CONTINUOUS, name=f"SoC_{t}_{i}")

# Constraint: SoC at t = 0 is 0.2

for i in range(I):

    model.addConstr(SoC[0, i] == SoC_min)
   

for t in range(T):

    for i in range(I):

        if PV_generation[t,i]-D_InflexibleDemand[t,i]>=0:

            model.addConstr(b[t, i] == 1)

        else:

            model.addConstr(a[t, i] == 1)

# Constraints

for t in range(T):

    for i in range(I):

        model.addConstr(P_pv[t, i] == PV_generation[t, i])

        model.addConstr(P_FlexibleEnergyDemand[t, i] <= f_inv(SoC_max - SoC[t, i]-a[t,i]*f(battDischarge[t,i])+b[t,i]*f(battCharge[t,i])))

        model.addConstr(P_FlexibleEnergyExcess[t, i] <= f_inv(SoC[t, i] - SoC_min-a[t,i]*f(battDischarge[t,i])+b[t,i]*f(battCharge[t,i])))

        model.addConstr(P_FlexibleEnergyExcess[t, i] >= 0)

        model.addConstr(P_FlexibleEnergyDemand[t, i] >= 0)

        model.addConstr(battCharge[t, i] <= abs(PV_generation[t,i]-D_InflexibleDemand[t,i]))

        model.addConstr(battDischarge[t, i] <= abs(PV_generation[t,i]-D_InflexibleDemand[t,i]))

        model.addConstr(0 <= P_INFlexibleEnergyDemand[t, i])

        model.addConstr(0 <= P_INFlexibleEnergyExcess[t, i])


for t in range(T+1):

    for i in range(I):        

        model.addConstr(SoC[t, i] >= SoC_min)

        model.addConstr(SoC[t, i] <= SoC_max)

# Sum of decision variables

for t in range(T):

    model.addConstr(quicksum(x[t, i] + y[t, i] for i in range(I)) == 1)

    model.addConstr(quicksum(a[t, i] + b[t, i] for i in range(I)) == 1)

# Constraints related to inflexible energy demand and SoC update equation

for i in range(I):

    for t in range(T):

        model.addConstr(SoC[t+1, i] == SoC[t, i] - f(-x[t,i]*P_FlexibleEnergyDemand[t, i] + y[t,i]*P_FlexibleEnergyExcess[t, i]+a[t,i]*battDischarge[t,i]-b[t,i]*battCharge[t,i]))

        model.addConstr((SoC[t, i] - f(a[t,i]*battDischarge[t,i]-b[t,i]*battCharge[t,i]))<=SoC_max)

        model.addConstr((SoC[t, i] - f(a[t,i]*battDischarge[t,i]-b[t,i]*battCharge[t,i]))>=SoC_min)

        model.addConstr(SoC[t+1, i] >= 0.2)

        model.addConstr(SoC[t+1, i] <= 1)

        model.addConstr(P_pv[t, i] + a[t,i]* battDischarge[t,i] + x[t,i]* P_INFlexibleEnergyDemand[t, i]  == D_InflexibleDemand[t, i] + y[t,i]* P_INFlexibleEnergyExcess[t, i]+ b[t,i]* battCharge[t,i])

model.addConstr(

    quicksum(P_pv[t, i] + a[t,i]* battDischarge[t,i] + x[t,i]* P_INFlexibleEnergyDemand[t, i]   for t in range(T) for i in range(I)) ==

    quicksum(D_InflexibleDemand[t, i] + y[t,i]* P_INFlexibleEnergyExcess[t, i]+ b[t,i]* battCharge[t,i] for t in range(T) for i in range(I))

)

       

model.setObjective(

    quicksum(

        0*(a[t,i]*battDischarge[t,i])-20*(b[t,i]*battCharge[t,i])+

        x[t, i] * (C_lowBidImport * P_FlexibleEnergyDemand[t, i] + C_HighBidImport * P_INFlexibleEnergyDemand[t,i]) +

        y[t, i] * (-C_HighOfferPrice * P_FlexibleEnergyExcess[t, i] - C_lowOfferPrice * P_INFlexibleEnergyExcess[t,i])

       

        for t in range(T) for i in range(I)

    ),

    GRB.MINIMIZE

)

# Solve the optimization problem

model.optimize()

# Access the solution

if model.status == GRB.OPTIMAL:

    for t in range(T):

        for i in range(I):

            print(f"********Results for: ", t,"**************")

            print(f"P_pv_{t}_{i}:", P_pv[t, i].x)

            print(f"D_InflexibleDemand_{t}_{i}:", D_InflexibleDemand[t, i])

            print(f"P_FlexibleEnergyDemand_{t}_{i}:", P_FlexibleEnergyDemand[t, i].x)

            print(f"P_FlexibleEnergyExcess_{t}_{i}:", P_FlexibleEnergyExcess[t, i].x)

            print(f"battDisch{t}_{i}:", battDischarge[t, i].x)

            print(f"battCharge{t}_{i}:", battCharge[t, i].x)

            print(f"P_INFlexibleEnergyDemand_{t}_{i}:", P_INFlexibleEnergyDemand[t, i].x)

            print(f"P_INFlexibleEnergyExcess_{t}_{i}:", P_INFlexibleEnergyExcess[t, i].x)

            print(f"x_{t}_{i}:", x[t, i].x)

            print(f"y_{t}_{i}:", y[t, i].x)

            print(f"a_{t}_{i}:", x[t, i].x)

            print(f"b_{t}_{i}:", y[t, i].x)

            print(f"SoC_{t}_{i}:", SoC[t, i].x)

            obj = model.getObjective()

            print(obj.getValue())

else:

    print("Optimization problem is infeasible or unbounded.")
However i could get manually results that give better solution that the one proposed by solver. Can any one help me. Thank you
David Torres Sanchez
Gurobi Staff
February 13, 2024 08:11
If you can obtain a better solution you could try to fix the variables to these values (set LB and UB values) and see what happens.
Cheers,
David
Gwyneth Butera
February 13, 2024 17:00
On addition, you may want to refer to the Knowledge Base article : How do I diagnose a wrong result?
Gurobi giving wrong results (not optimal solution)

Comments

Didn't find what you were looking for?
