I am writing an optimization model in python, using Gurobi, with the objective to maximize the revenue of a solar power plant by using a battery storage. My goal is to optimize when I will sell the produced energy immediately to the power market or store it in a battery storage and sell it at a later hour when the market price is higher. However, after running my model and analyzing the results I have noticed that the model's solution is unconventional.
For example, in hour 7 I have a high price of 104.54 EUR and my battery is filled with 1 MWh. Instead of selling the whole 1 MWh in that hour, the model decided to sell only 0.5 MWh for that price. Then in the next hour when the prices is lower (103.39 EUR) sells additional 0.25. Then at next hour sells 0.125 MWh at price 95 EUR, the next hour 0.0625 MWh at price 80.15 EUR. I can't figure out what constraint causes a result like this unfortunately. So I was wondering if someone could maybe pinpoint where the problem is?
Here is part of my model:
num_hours = len(solar_production)
# Create a new Gurobi model
model = gp.Model("SolarPowerPlantWithBattery")
# Decision variables
energy_sold = model.addVars(num_hours, lb=0, ub=5, vtype=GRB.CONTINUOUS, name="EnergySold")
battery_status = model.addVars(num_hours, lb=0, ub=battery_capacity, vtype=GRB.CONTINUOUS, name="BatteryStatus")
battery_discharge = model.addVars(num_hours, lb=0, ub=battery_power, vtype=GRB.CONTINUOUS, name="BatteryDischarge")
battery_charge = model.addVars(num_hours, lb=0, ub=battery_power, vtype=GRB.CONTINUOUS, name="BatteryCharge")
binary_variable = model.addVars(num_hours, vtype=GRB.BINARY, name="BinaryVariable")
# Objective function: Maximize total revenue
model.setObjective(gp.quicksum(energy_sold[i] * market_price[i] for i in range(num_hours)), sense=GRB.MAXIMIZE)
# Battery charge and discharge constraints
model.addConstr(battery_charge[0] == 0)
model.addConstr(battery_discharge[0] == 0)
model.addConstr(battery_status[0] == battery_capacity/2) # Initial condition for the battery charge
model.addConstrs(battery_status[i] == battery_status[i - 1] - battery_discharge[i] + battery_charge[i] for i in range(1, num_hours))
model.addConstrs(battery_discharge[i] <= battery_status[i] for i in range(0,num_hours))
# Battery discharge constraint
model.addConstr(gp.quicksum(battery_discharge[i] for i in range(0, num_hours)) <= discharge_limit*num_hours/24)
# Mutual exclusion constraint
model.addConstrs(battery_charge[i] <= battery_power * (1 - binary_variable[i]) for i in range(num_hours))
model.addConstrs(battery_discharge[i] <= battery_power * binary_variable[i] for i in range(num_hours))
# Connection between solar production, energy sold, and battery charge
model.addConstrs(solar_production[i] - energy_sold[i] - battery_charge[i] + battery_discharge[i] == 0 for i in range(num_hours))
# Optimize the model
model.optimize()
I have timeseries for the solar_production and market_price on hourly frequency. The battery_capacity is set at 2, the battery_power at 1 and the discharge_limit at 2.
I went through all my constraints multiple times and I can't figure out why is this happening. Also tried to force my model to discharge the battery less often by subtracting the binary value (1 when the battery is discharging) from the optimization criteria, it just resulted in everlasting calculation, which I had to interrupt.
I figured out what the problem was. It was this line:
model.addConstrs(battery_discharge[i] <= battery_status[i] for i in range(0,num_hours))
I am always looking the values at the end of the hour, so naturally at the end of hour 7, if the battery was discharged for 1MWh, my battery_status would have been 0MWh, which clearly breaks the constraint above.