Need help with gurobi code implementation for battery in microgrid

I am trying to implement the battery storge system within the microgrid environment using gurobi interfaced with python. 

As you have mentioned that, you used gurobi iterfaced with python in your paper, so I need your help for my code. 

Actually, it works fine without storage, i.e. Power is sell to the grid when:

total_load>=total power from all generators + wind+solar

and power is bought from the grid when:

  total_load<=total power from all generators + wind+solar  

but with storage the SOC status is not updated, it either goes to min of SOC or max of SOC.

Kindly help me around this.

Best regards,

Misbah Rani,

student of MSEE

 

 

 

 

CODE

# -*- coding: utf-8 -*-
"""
Created on Sun Sep 15 15:17:16 2019

@author: Misbah
"""

from gurobipy import *

#constants of the model
#revenue from serving load

#Load
PL=[15,15,15,15,15,15,15,15,15,15,15,15,15,15,15,15]
#FuelCell

#MicroTurbine
C_g= 0.09# natural gas price(currency/kWh)
C_d= 1.2# diesel price(currency/l)
#MT_1
a_DG1=2e-3
a_DG2=1.115e-3
a_FC=9.864e-3
a_MT1=1.924e-3
a_MT2=0
b_DG1=0.194
b_DG2=0.2232
b_FC=2.061
b_MT1=2.232
b_MT2=2.965
c_DG1=1.85
c_DG2=1.75
c_FC=1.152
c_MT1=82.33
c_MT2=32.2
#maintenance cost
CM_FC_c=0.004
CM_MT1_c=0.006
CM_MT2_c=0.015
CM_DG1_c=0.005
CM_DG2_c=0.02
#effiiencies
eta_FC_min=0.60
eta_DG1_min=0.50
eta_DG2_min=0.50
eta_MT1_min=0.32
eta_MT2_min=0.32
eta_c=0.85
eta_d=0.95

#generation limits
P_DG1_max=20
P_DG2_max=70
P_FC_max=50
P_MT1_max=200
P_MT2_max=65
P_DG1_min=5
P_DG2_min=7
P_FC_min=5
P_MT1_min=20
P_MT2_min=6.5
P_PV=0
P_wind=15
#cost of grid
C_grid_buy=20
C_grid_sell=15
#battery cost
CS_c=25
CS_d=18

R=5/100
P_grid_max=100
P_grid_min=0
Battery_rating=60 #A
SOC_min=0.3*Battery_rating
SOC_max=.9*Battery_rating
PS_min=0.10*Battery_rating
PS_max=0.20*Battery_rating
SOC_0=0.8*Battery_rating
PS_c_0=40
PS_d_0=10
C_L=10
# Model
model = Model('microgird_grid_connected')

# Create decision variables for the foods to buy
P_DG1 = model.addVar(vtype=GRB.INTEGER,name="P_DG1")
P_DG2 = model.addVar(vtype=GRB.INTEGER,name="P_DG2")
P_FC = model.addVar(vtype=GRB.INTEGER,name="P_FC")
P_MT1 = model.addVar(vtype=GRB.INTEGER,name="P_MT1")
P_MT2 = model.addVar(vtype=GRB.INTEGER,name="P_MT2")

P_buy = model.addVar( name="P_buy")
P_sell = model.addVar( name="P_sell")
PS_c = model.addVar(vtype=GRB.INTEGER,name="PS_c")
PS_d = model.addVar(vtype=GRB.INTEGER,name="PS_d")
SOC = model.addVar(vtype=GRB.INTEGER,name="SOC")
U_DG1 = model.addVar(vtype=GRB.BINARY,name="U_DG1")
U_DG2 = model.addVar(vtype=GRB.BINARY,name="U_DG2")
U_MT1 = model.addVar(vtype=GRB.BINARY,name="U_MT1")
U_MT2 = model.addVar(vtype=GRB.BINARY,name="U_MT2")
U_FC = model.addVar(vtype=GRB.BINARY,name="U_FC")
U_c= model.addVar(vtype=GRB.BINARY,name="U_c")
U_d= model.addVar(vtype=GRB.BINARY,name="U_d")

U_buy= model.addVar(vtype=GRB.BINARY,name="U_buy")
U_sell= model.addVar(vtype=GRB.BINARY,name="U_sell")

#expressions
#cost of DG1
C_DG1=C_d*(a_DG1*P_DG1*P_DG1+b_DG1*P_DG1 +c_DG1)
#cost of DG2
C_DG2=C_d*(a_DG2*P_DG2*P_DG2+b_DG2*P_DG2 +c_DG2)
#cost of fuel cell
C_FC=C_g*(a_FC*P_FC*P_FC+b_FC*P_FC+c_FC)
#cost of MT1
C_MT1=C_d*(a_MT1*P_MT1*P_MT1+b_MT1*P_MT1 +c_MT1)
#cost of MT2
C_MT2=C_g*(a_MT2*P_MT2*P_MT2+b_MT2*P_MT2 +c_MT2)
#C&M cost of DG1
CM_DG1=CM_DG1_c*P_DG1
#C&M cost of DG2
CM_DG2=CM_DG2_c*P_DG2
#C&M cost of FC
CM_FC=CM_FC_c*P_FC
#C&M cost of MT1
CM_MT1=CM_MT1_c*P_MT1
#C&M cost of MT2
CM_MT2=CM_MT2_c*P_MT2
#Startup cost of DG1
C_up_DG1 = C_d*P_DG1_min/eta_DG1_min
#Startup cost of DG2
C_up_DG2 = C_d*P_DG2_min/eta_DG2_min
#Startup cost of MT1
C_up_MT1 = C_g*P_MT1_min/eta_MT1_min
#Startup cost of MT2
C_up_MT2 = C_g*P_MT2_min/eta_MT2_min
#Startup cost of FC
C_up_FC = C_g*P_FC_min/eta_FC_min
#Cost of seving load
C_load=C_L*sum(PL)
#Cost of op
C_G_op=C_DG1+C_DG2+C_FC+C_MT1+C_MT2
#Cost of C&M
C_G_CM=CM_DG1+CM_DG2+CM_FC+CM_MT1+CM_MT2
#cost of startup
C_G_up=C_up_DG1+C_up_DG2+C_up_FC+C_up_MT1+C_up_MT2
#total cost of generation
C_G_T=C_G_op+C_G_CM+C_G_up

#cost of grid
C_grid=C_grid_buy*P_buy - C_grid_sell*P_sell
P_G_T= P_DG1+ P_DG2+ P_FC+ P_MT1+ P_MT2
PL_T=sum(PL)
P_grid= U_buy*P_buy-U_sell*P_sell
C_S = CS_d*PS_d*U_d-CS_c*PS_c*U_c
if SOC_0 >= SOC_min and SOC_0 <= SOC_max-1:
U_c=1
elif SOC_0 >= SOC_max:
U_d=1

P_grid= U_buy*P_buy-U_sell*P_sell
#SOC = SOC_0+eta_c*PS_c*U_c-PS_d*U_d/eta_d
PS=PS_d*U_d-PS_c*U_c

#constraints
model.addConstr(U_c+U_d<=1)
model.addConstr(U_c*U_d<=0)
model.addConstr(PS_c<=PS_max*U_c)
model.addConstr(PS_c>=PS_min*U_c)

model.addConstr(PS_d<=PS_max*U_d)
model.addConstr(PS_d>=PS_min*U_d)
model.addConstr(SOC <= SOC_0 +eta_c*PS_c-PS_d/eta_d)
model.addConstr(SOC>=SOC_min)
model.addConstr(SOC<=SOC_max)

 

model.addConstr(U_sell+U_buy<=1)
model.addConstr(P_DG1<=P_DG1_max*U_DG1)
model.addConstr(P_DG1>=P_DG1_min*U_DG1)

model.addConstr(P_buy<=P_grid_max*U_buy)
model.addConstr(P_buy>=P_grid_min*U_buy)

model.addConstr(P_sell<=P_grid_max*U_sell)
model.addConstr(P_sell>=P_grid_min*U_sell)
model.addConstr(P_DG2<=P_DG2_max*U_DG2)
model.addConstr(P_DG2>=P_DG2_min*U_DG2)

model.addConstr(P_FC<=P_FC_max*U_FC)
model.addConstr(P_FC>=P_FC_min*U_FC)

model.addConstr(P_MT1<=P_MT1_max*U_MT1)
model.addConstr(P_MT1>=P_MT1_min*U_MT1)

model.addConstr(P_MT2<=P_MT2_max*U_MT2)
model.addConstr(P_MT2>=P_MT2_min*U_MT2)

model.addConstr( P_grid +PS+P_G_T+P_PV+P_wind>=sum(PL))
obj=C_load-C_G_T-C_grid-C_S

model.setObjective((obj), GRB.MAXIMIZE)
#model.setObjective((obj), GRB.MINIMIZE)
model.optimize()