I am trying to model the battery charging/discharging process with optimization program in Gekko. I used the following energy balance equation for the model:
Renewable power [t] - (charge[t] + surplus[t])+(discharge[t] + unmet[t])= Demand[t]
The state of charge/discharge (SOC
) is updated via the following equation:
SOC[t]=SOC[t-1]+charge[t]-discharge[t]
I formulated the Equations above in the Gekko
, but came across a problem in the discharging phase. For example, I do not receive a solution for the second row (yellow part in the image) of the data set shown in the image below. I expect that the SOC[second row]=196400
, but I do not receive a solution in the Gekko
model. I appreciate your advice in the case of this problem.
from gekko import GEKKO
import sys
Battery_capacity=589200 #kW-h
Battery_Initial_capacity=589200
max_discharge=Battery_capacity/3
max_charge=Battery_capacity
Big_M=1000000000000
epsilon=sys.float_info.epsilon
m = GEKKO(remote=False)
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1',\
'minlp_integer_tol 1e-100', \
'minlp_integer_leaves 2']
Demand_total=[5629940,5327262,4669939]
RE_total=[5204098,4410247,3694715]
Price_total=[15,17,17]
row=1
finish=0
counter=0
for r in range(3):
start=finish
finish=(r+1)*row
Demand = [Demand_total[r]]
RE = [RE_total[r]]
Price = [Price_total[r]]
# Objective function
c = [None] * row
for t in range(row):
c[t] = (Demand[t] * Price[t])
B_ch_over = [m.Var(lb=0) for t in range(row)]
B_disch_over = [m.Var(lb=0) for t in range(row)]
B_ch = [m.Var(lb=0) for t in range(row)]
B_disch = [m.Var(lb=0) for t in range(row)]
Unmet = [m.Var(lb=0) for t in range(row)]
Surplus = [m.Var(lb=0) for t in range(row)]
B_SOC = [m.Var() for t in range(row)]
if r==0:
starting_point = 0
Last_BOS=Battery_Initial_capacity
m.Equation(B_SOC[0]==Last_BOS) # initial charge at start
else:
m.Equation(B_SOC[0] == Last_BOS[0]) # initial charge at start
for t in range(starting_point,row):
m.Equation(RE[t] - B_ch_over[t] + B_disch_over[t] == Demand[t])
m.Equation(B_ch_over[t] * B_disch_over[t] <= 0)
B_ch[t] = m.if3(max_charge - m.if3(B_ch_over[t] - (Battery_capacity - B_SOC[t - 1]), B_ch_over[t],Battery_capacity - B_SOC[t - 1]), max_charge,m.if3(B_ch_over[t] - (Battery_capacity - B_SOC[t - 1]), B_ch_over[t],Battery_capacity - B_SOC[t - 1]))
B_disch[t]=m.if3(max_discharge-m.if3(B_disch_over[t]-B_SOC[t-1],B_disch_over[t],B_SOC[t-1]),max_discharge,m.if3(B_disch_over[t]-B_SOC[t-1],B_disch_over[t],B_SOC[t-1]))
B_SOC[t]=m.if3(Battery_capacity-(B_SOC[t - 1] + B_ch[t] - B_disch[t]),Battery_capacity,B_SOC[t - 1] + B_ch[t] - B_disch[t])
m.Equation(B_ch[t] + Surplus[t] == B_ch_over[t])
m.Equation(B_disch[t] + Unmet[t] == B_disch_over[t])
counter+=1
m.Maximize(m.sum(c))
m.options.SOLVER = 1
# m.open_folder()
m.solve(disp=False)
Last_BOS=B_SOC[t].value
if r==0:
print('Results')
for t in range(starting_point, row):
print('B_ch{}={} B_disch{}={} B_SOC{}={} B_ch_over{}={} B_disch_over{}={} Unmet{}={} Surplus{}={} RE{}={} Demand{}={}'.format(counter-row+t, B_ch[t].value, counter-row+t, B_disch[t].value, counter-row+t, B_SOC[t].value, counter-row+t, B_ch_over[t].value, counter-row+t,B_disch_over[t].value, counter-row+t, Unmet[t].value, counter-row+t, Surplus[t].value, counter-row+t, RE[t], counter-row+t,Demand[t]))
print('Objective: ' + str(m.options.objfcnval))
I used two additional variables (B_ch_if2 & B_disch_if2) and decomposed the nested m.if3 of the charging and discharging sections and received the results. However, I did not notice why the nested m.if3 did not work in the initial code. I appreciate your comments about this issue. Thanks
from gekko import GEKKO
import sys
Battery_capacity=589200 #kW-h
Battery_Initial_capacity=589200
max_discharge=Battery_capacity/3
max_charge=Battery_capacity
Big_M=1000000000000
epsilon=sys.float_info.epsilon
m = GEKKO(remote=False)
m.solver_options = ['minlp_gap_tol 1.0e-2',\
'minlp_maximum_iterations 10000',\
'minlp_max_iter_with_int_sol 500',\
'minlp_branch_method 1',\
'minlp_integer_tol 1e-100', \
'minlp_integer_leaves 2']
Demand_total=[5629940,5327262,4669939]
RE_total=[5204098,4410247,3694715]
Price_total=[15,17,17]
row=1
finish=0
counter=0
for r in range(3):
start=finish
finish=(r+1)*row
Demand = [Demand_total[r]]
RE = [RE_total[r]]
Price = [Price_total[r]]
# Objective function
c = [None] * row
for t in range(row):
c[t] = (Demand[t] * Price[t])
B_ch_over = [m.Var(lb=0) for t in range(row)]
B_disch_over = [m.Var(lb=0) for t in range(row)]
B_ch = [m.Var(lb=0) for t in range(row)]
B_disch = [m.Var(lb=0) for t in range(row)]
Unmet = [m.Var(lb=0) for t in range(row)]
Surplus = [m.Var(lb=0) for t in range(row)]
B_ch_if2 = [m.Var(lb=0) for t in range(row)]
B_disch_if2 = [m.Var(lb=0) for t in range(row)]
B_SOC = [m.Var() for t in range(row)]
if r==0:
starting_point = 0
Last_BOS=Battery_Initial_capacity
m.Equation(B_SOC[0]==Last_BOS) # initial charge at start
else:
m.Equation(B_SOC[0] == Last_BOS[0]) # initial charge at start
for t in range(starting_point,row):
m.Equation(RE[t] - B_ch_over[t] + B_disch_over[t] == Demand[t])
m.Equation(B_ch_over[t] * B_disch_over[t] <= 0)
B_ch_if2[t] = m.if3(B_ch_over[t] - (Battery_capacity - B_SOC[t - 1]), B_ch_over[t],Battery_capacity - B_SOC[t - 1])
B_ch[t] = m.if3(max_charge - B_ch_if2[t], max_charge, B_ch_if2[t])
B_disch_if2[t] = m.if3(B_disch_over[t] - B_SOC[t - 1], B_disch_over[t], B_SOC[t - 1])
B_disch[t] = m.if3(max_discharge - B_disch_if2[t], max_discharge, B_disch_if2[t])
B_SOC[t]=m.if3(Battery_capacity-(B_SOC[t - 1] + B_ch[t] - B_disch[t]),Battery_capacity,B_SOC[t - 1] + B_ch[t] - B_disch[t])
m.Equation(B_ch[t] + Surplus[t] == B_ch_over[t])
m.Equation(B_disch[t] + Unmet[t] == B_disch_over[t])
counter+=1
m.Maximize(m.sum(c))
m.options.SOLVER = 1
# m.open_folder()
m.solve(disp=False)
Last_BOS=B_SOC[t].value
if r==0:
print('Results')
for t in range(starting_point, row):
print('B_ch{}={} B_disch{}={} B_SOC{}={} B_ch_over{}={} B_disch_over{}={} Unmet{}={} Surplus{}={} RE{}={} Demand{}={}'.format(counter-row+t, B_ch[t].value, counter-row+t, B_disch[t].value, counter-row+t, B_SOC[t].value, counter-row+t, B_ch_over[t].value, counter-row+t,B_disch_over[t].value, counter-row+t, Unmet[t].value, counter-row+t, Surplus[t].value, counter-row+t, RE[t], counter-row+t,Demand[t]))
print('Objective: ' + str(m.options.objfcnval))