batterygekko

Battery State of Charge (SOC) updating in Gekko


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.

enter image description here

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))

Solution

  • 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))